Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
123,512 characters · 40 sections · 149 citation commands
Recent Developments on Factor Models and its Applications in Econometric Learning
\singlespacing
{ Key words: factor models, spiked low rank matrix, matrix completion, unbalanced panel, multiple testing, high-dimensional
}
\eject
\onehalfspacing
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 fan2020statistical. 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 stock2016dynamic for dynamic factor models with applications on macroeconomics, to bai2016econometric for time series and panel data models, and to gagliardini2019estimation 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 giannone2008nowcasting and doz2011two 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 brillinger1964frequency. forni, forni2005generalized 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: bai2002determining, AH, onatski2010determining, li2017determining, bai2012statistical, bai2016maximum, onatski2012asymptotics, chudik2011weak, cheng2016shrinkage, massacci2017least, gagliardini2016time, gonccalves2018bootstrapping, baltagi2017identification,barigozzi2018simultaneous, ait2017using, chen2019five, liao2018uniform,li2019jump, pelger2019large, su2017time.
We use the following notation. For a matrix $\bfm A$, let $\lambda_i(\bfm A)$ denote the $i$ th largest singular value of $\bfm A$ and use $\lambda_{\min}(\bfm A)$ and $\lambda_{\max}(\bfm A)$ to denote its smallest and largest eigenvalues. We define the Frobenius norm $\|\bfm A\|_F=\sqrt{\operatorname{tr}(\bfm A'\bfm A)}$, the operator norm $\|\bfm A\|=\sqrt{\lambda_{\max}(\bfm A'\bfm A)}$, the element-wise norm $\|\bfm A\|_\infty=\max_{ij}|A_{ij}|$, and the matrix $\ell_1$-norm $\|\bfm A\|_{\ell_1}:=\max_{i\leq N}\sum_{j=1}^N|A_{ij}|$. In addition, define projection matrices $\bfm P_\bfm A=\bfm A(\bfm A'\bfm A)^{-1}\bfm A$ and $\bfm M_\bfm A=\bfm I-\bfm P_\bfm A$ when $\bfm A'\bfm A$ is invertible. Finally, for two (random) sequences $a_T$ and $b_T$, we write $a_T\ll b_T$ (or $b_T\gg a_T$) if $a_T=o_P(b_T)$.
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 $\bfsym \Sigma$ (either observable or not), having the following decomposition:
Such decomposition satisfies the following three properties:
The low-rank structure achieves dimension reductions: suppose the matrix $\bfsym \Sigma$ is of $N\times N_1$ dimensions, while the rank of $\bfm 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 $\bfm L$. Meanwhile, the spikedness helps seperate $\bfm L$ from $\bfm S$ approximately, and ensures that the large “signals" concentrate on $\bfm 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 $\bfm f_t$ is a $r$-dimensional vector of factors; $\bfm b_i$ is the loading vector and $u_{it}$ is the idiosyncratic noise. Specifically, ((ref)) applies to two decompositions of this model.
Factor Decomposition. The matrix form of the factor model gives $$ \bfm Y=\bfm M+\bfm U,\quad \bfm M:=\bfm B\bfm F', $$ where $\bfm Y$ and $\bfm U$ are $N\times T$ matrices of $y_{it}$ and $u_{it}$; $\bfm B$ is the $N\times r$ matrix of $\bfm b_i$ while $\bfm F$ is the $T\times r$ matrix of $\bfm f_t.$ Then corresponding to the notation ((ref)), $\bfsym \Sigma=\bfm Y$, $\bfm L=\bfm M$ and $\bfm S=\bfm U.$ In this decomposition, $\bfsym \Sigma$ is observable. Apparently, $\bfm M$ is a low-rank matrix with rank $r$. The nonzero singular values of $\bfm M$, under the strong factor assumption, grows much faster than those of $\bfm U$, which gives rise to the spikedness property. Now let $\ensuremath{\boldsymbol{\xi}}$ be the $N\times r$ matrix whose columns are the left singular vectors of $\bfm M$, and let $\ensuremath{\boldsymbol{\xi}}_i'$ denote its $i$ th row. Then under the assumption that the nonzero eigenvalues of $\bfm B'\bfm 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 ((ref)) that the covariance matrix of $\bfm y_t=(y_{1t},\cdots,y_{Nt})'$, denoted by $\bfsym \Sigma_y$, can be decomposed as follows:
where $\bfsym \Sigma_u$ denotes the covariance matrix of $\bfm u_t$. The above decomposition is well known for portfolio allocations and risk managements, where the total volatility is decomposed into the systematic risk $\bfm L$, plus the (sparse) idiosyncratic risk $\bfsym \Sigma_u$. It also leads to the spiked incoherent low-rank model, but $\bfsym \Sigma_y$ is unknown and needs to be estimated.
There are two general approaches to estimating model ((ref)): (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 (ref). We shall assume rank$(\bfm L)=r$ to be known.
For any matrix $\bfm A$, let $\bfm A=\bfm U_{A}\bfm D_{A}\bfm V'_{A}$ denote the singular value decomposition (SVD) of $\bfm A$. Define the singular value hard thresholding operator as
where $\bar \bfm D_R$ is a diagonal matrix that keeps the top $R$ diagonal elements of $\bfm D_A$ and replaces the remaining elements by zeros. So $H_R(\bfm A)$ is the best rank $R$ matrix approximation to $\bfm A$.
Suppose an estimator of $\bfsym \Sigma$, denoted by $\widehat\bfsym \Sigma$, is available, satisfying
for some sequences $\eta_N$ and $ c_N$. We use $\widehat\bfsym \Sigma$ as the input matrix, which can be the sample covariance matrix or its robustfied versions fan2019robust. The goal is to estimate $\bfm L$ in ((ref)) and its $N\times r$ matrix of the left singular vectors, denoted by $\ensuremath{\boldsymbol{\xi}}$ (also let $\ensuremath{\boldsymbol{\zeta}}$ denote its right singular vectors). We use respectively $ \widehat\bfm L:=H_R(\widehat\bfsym \Sigma) $ with $R=r$, which is the rank $r$ projection of $\widehat{\bfsym \Sigma}$, and the $N\times r$ matrix $\widehat\ensuremath{\boldsymbol{\xi}}$ whose columns are the left singular vectors of $\widehat\bfsym \Sigma$. The following theorem, adapted from fan2018eigenvector, provides deviation bounds of the estimators. To make the paper self-contained, we also provide a simpler proof with slightly different conditions.
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.
We observe an $N\times T$ data matrix $\bfm Y$, which can be decomposed as \[\bfm Y=\bfm M+\bfm U=\bfm B\bfm F'+\bfm U\] where $\bfm B$ is $N\times r$ factor loadings matrix, $\bfm F$ is $T\times r$ factors matrix and $\bfm U$ is $N\times T$ idiosyncratic errors, which are uncorrelated with $\bfm M:=\bfm B \bfm F'$. All the three parts $\bfm B$, $\bfm F$ and $\bfm U$ are unobserved. The $t$ th column of this expression can be written as
Under the model's specification, we have the covariance structure ((ref)). One of the most widely used estimation methods for the factor model is principal components analysis (PCA). Define the sample covariance $\bfm S_y=\frac1T\sum_{t=1}^T \bfm y_t\bfm y_t'=\frac1T\bfm Y\bfm Y'$. Let $\widehat\ensuremath{\boldsymbol{\xi}}_j$ be the $j$th eigenvector corresponding to the largest $j$ th eigenvalues of $\bfm S_y$. The PCA estimates $\bfm B$ by taking $\widehat\bfm B=\sqrt{N}(\widehat\ensuremath{\boldsymbol{\xi}}_1,\cdots,\widehat\ensuremath{\boldsymbol{\xi}}_R)$, which estimates $\bfm B$ up to a diagonal transformation. Given $\widehat\bfm B$, the factors can be estimated via the least squares: $$ \widehat\bfm F=\bfm Y'\widehat\bfm B (\widehat\bfm B'\widehat\bfm B)^{-1} =\frac{1}{N}\bfm Y'\widehat\bfm B. $$ This also leads to the estimated low-rank component $\frac{1}{T}\widehat\bfm B\widehat\bfm F'\widehat\bfm F\widehat\bfm B'$ for $\bfm B\operatorname{cov}(\bfm f_t)\bfm B' $.
PCA is equivalent to the singular value hard thresholding by taking the input matrix $\widehat\bfsym \Sigma=\bfm S_y$. Then $\frac{1}{T}\widehat\bfm B\widehat\bfm F'\widehat\bfm F\widehat\bfm B'=H_R(\bfm S_y)$. One can then apply Theorem (ref) to infer the rates of convergence of the PCA estimators, which were obtained by SW02. bai03 proved the asymptotic normality of PCA estimators for the factors and loadings. Results with general input $\widehat{\bfsym \Sigma}$ can be found in Chapter 10 of fan2020statistical.
Another popular method to estimate a factor model is the maximum likelihood (ML) method (see, e.g., Lawley, bai2012statistical, doz). Under the independence and normality assumptions, the log-likelihood function based on $\bfm y_t$ is, for some constant $C$, \[ \log L_{\mathrm{ML}}(\bfm B, \operatorname{cov}(\bfm f_t), \operatorname{diag}(\bfsym \Sigma_u))=C-\frac T2\ln|\bfsym \Sigma_y|-\frac12\sum_{t=1}^T\bfm y_t'\bfsym \Sigma_y^{-1} \bfm y_t. \] The log-likelihood function is then maximized with respect to the matrix parameters $(\bfm B, \operatorname{cov}(\bfm f_t), \operatorname{diag}(\bfsym \Sigma_u))$ under additional restrictions that $\bfsym \Sigma_u$ is diagonal bai2012statistical, bai2016maximum or sparse with regularizations bai2016efficient, wang2019penalized. Recently barigozzi2019quasi 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 $\bfm f_t$ and $\bfm u_t$, we have \[\mathbb E(\bfm f_t|\bfm y_t)=\bfm B'(\bfm B\bfm B'+\bfsym \Sigma_u)^{-1}\bfm y_t=(\bfm I_r+\bfm B'\bfsym \Sigma_u^{-1}\bfm B)^{-1}\bfm B'\bfsym \Sigma_u^{-1}\bfm y_t.\] This provides the basis of estimating factors. The other approach is the generalized least squares: for given $\bfm B$ and $\bfsym \Sigma_u^{-1}$, the GLS estimator for $\bfm f_t$ is \[\widehat\bfm f_t= (\bfm B'\bfsym \Sigma_u^{-1}\bfm B)^{-1}\bfm B'\bfsym \Sigma_u^{-1}\bfm y_t.\] 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.
Alternative to PCA, one can estimate $\bfm M$ directly taking advantage of its low-rank structure, based on the nuclear-norm regularization, the $\ell_1$-norm of singular values, that encourages the sparseness in singular values and hence low-rankness. For an $n\times m$ matrix $\bfm A$, let $\|\bfm A\|_n:=\sum_{i=1}^{\min\{m,n\}}\psi_i(\bfm A)$ be its nuclear-norm, where $\psi_i(\bfm A)$ is the $i$ th largest singular value of $\bfm A$.
Given the low-rank structure of $\bfm M$ (sparsity in singular value of $\bfm M$), we can estimate the model via solving the following penalized optimization:
for some tuning parameter $\nu>0$. The solution is $\widehat \bfm M= S_{\nu}(\bfm Y)$, where $S_{\nu}(\cdot)$ is the singular value thresholding operator ma2011fixed, defined as follows. Let $\bfm Y=\bfm U_y\bfm D\bfm V_y'$ be its SVD. Then $S_{\nu}(\bfm Y):= \bfm U_y\bfm D_{\nu}\bfm V_y',$ where $\bfm D_{\nu}= \operatorname{diag}(\{D_{ii}-\nu\}_+)$ with $D_{ii}$ being the diagonal entries of $\bfm D$. So $S_{\nu}(\bfm Y)$ applies “soft-thresholding" on the singular values of $\bfm 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 $\bfm M$ with $R$ factors is given by (see Section (ref)): $$ \widehat \bfm M_{\text{PC}} = H_R(\bfm Y),\quad H_R(\bfm Y):=\bfm U_y\bar \bfm D_{R} \bfm V_y'. $$ This estimator is the solution to the penalized least squares problem ((ref)) except that the nuclear norm is replaced by $\sum_{i=1}^{\min\{N,T\}} p_\nu(\psi_i(\bfm M))$, where $p_\nu(\theta) = \nu^2 - (\nu - |\theta|)_+^2$ is the harding thresholding penalty and $\psi_i(\bfm M)$ is the $i^{th}$ singular value of $\bfm M$.
Therefore the difference between ((ref)) 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 $\bfm Y.$ An important implication is that the factor estimator building on $\widehat\bfm 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 fan2019distributed.
Recall that $\bfsym \Sigma_y$ and $\bfsym \Sigma_u$ denote the $N\times N$ covariance matrices of $\bfm y_t$ and $\bfm u_t$ in model ((ref)), 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, $\bfsym \Sigma_u$ is sparse. While the definition of sparsity may differ in different contexts, here we mean $$ J:=\sum_{i\neq j} 1\{\operatorname{\mathbb E} u_{it} u_{jt}\} $$ should not grow too fast as $N\to\infty.$ This requirement can be weakened to approximate sparsity. In addition, $\bfm L$ is a low-rank matrix. Thus we can directly estimate the above covariance decomposition via solving the following penalized optimization:
where $\nu_1$ and $\nu_2$ are tuning parameters. Note that here we use the notation $\|\bfm A\|_1=\sum_{i,j}|A_{ij}|$ as the matrix 1-norm, distinguished from the usual matrix $\ell_1$-norm $\|\bfm A\|_{\ell_1}:=\max_{i\leq N}\sum_{j=1}^N|A_{ij}|$. The above optimization has been employed by many authors to study the low rank plus sparse decomposition, while some authors exclude the diagonal elements of $\bfsym \Sigma_u$ from the penalization, and additionally impose positive-definite and other constraints on $\bfm L$ and $\bfsym \Sigma_u$ klopp2017robust, agarwal2012noisy. Finally, given $\widehat\bfm L$, we can estimate the factors and loadings by extracting its eigenvectors.
The above optimization can be solved by alternating the estimation of $\bfm L$ and $\bfsym \Sigma_u$, and closed form solutions are available in both iterations. Given $\bfsym \Sigma_u$, solving for $\bfm L$ leads to the singular value soft-thresholding: $\widehat{\bfm L} = S_{\nu_1}(\bfm S_y- \bfsym \Sigma_u)$, and given $\bfm L$, solving for $\bfsym \Sigma_u$ leads to the element-wise soft-thresholding: $\widehat{\bfsym \Sigma}_u = \widetilde S_{\nu_2}(\bfm S_y-\bfm L)$. While both iterations solve convex problems, standard convergence analysis can be applied to show that the iterative algorithm converges in polynomial time.
agarwal2012noisy and klopp2017robust studied the statistical convergence properties of ((ref)). Let columns of $\bfm U_{L,2}$ be the singular vectors of the true $\bfm L$ corresponding to the zero singular values. Define projections $\mathcal P(\bfm A):=\bfm U_{L,2}\bfm U_{L,2}'\bfm A\bfm U_{L,2}\bfm U_{L,2}'$ and $\mathcal M(\bfm A):=\bfm A-\mathcal P(\bfm A)$. In addition, let $(\bfm A)_J$ and $(\bfm A)_{J^c}$ be the submatrices of $\bfm A$, whose elements respectively correspond to $\operatorname{\mathbb E} u_{it}u_{jt}\neq0$ and $\operatorname{\mathbb E} u_{it}u_{jt}=0.$ Additionally define $$ \mathcal C(\nu_1, \nu_2):=\{(\bfm A_1,\bfm A_2): \nu_1\|\mathcal P(\bfm A_1)\|_n+\nu_2\|(\bfm A_2)_{J^c}\|_1\leq 3\nu_1\|\mathcal M(\bfm A_1)\|_n+3\nu_2\|(\bfm A_2)_{J}\|_1 \}. $$ A key quantity is the restricted strong convexity (RSC) constant, which is defined as follows: $$ \kappa(\nu_1,\nu_2):=\sup\{c>0:\|\bfm A_1+\bfm A_2\|_F^2\geq c\|\bfm A_1\|_F^2+c\|\bfm A_2\|_F^2 \text{ for all } (\bfm A_1, \bfm A_2)\in \mathcal C(\nu_1, \nu_2) \}. $$
We then have the following theorem, adapted from agarwal2012noisy. To make the paper self-contained, we also provide a proof with slightly different conditions. See the online supplement.
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: $$ \|\bfm S_y-\bfsym \Sigma_y\|\leq \nu_1,\quad \|\bfm S_y-\bfsym \Sigma_y\|_\infty\leq \nu_2; $$ 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 (ref) 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 agarwal2012noisy for more refined lower bound of $\kappa(\nu_1,\nu_2)$.
The above problem is also called “robust PCA” candes2011robust. For recent advance and references, see chen2020bridging where factorization methods are also discussed.
POET proposed a nonparametric estimator of $\bfsym \Sigma_y$, named POET (Principal Orthogonal complEment Thresholding), when the factors are unobservable. It is basically an one-step solution to optimization ((ref)) with initialization $\bfsym \Sigma_u=0$. To motivate the estimator, suppose $r=R$. Then, heuristically $$ \bfm L\approx H_R(\bfsym \Sigma_y),\quad \bfsym \Sigma_u\approx \bfsym \Sigma_y -H_R(\bfsym \Sigma_y), $$ Thus, one estimates $\bfm L$ by $H_R(\bfm S_y)$ and sets $\bfm S_u:=\bfm S_y-H_R(\bfm S_y)$. To account for the sparsity assumption on $\bfsym \Sigma_u$, POET estimates $\bfsym \Sigma_y$ and $\bfsym \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 $\bfm S_u$POET; we can also take other form such as the adaptive thresholding in Cai11b. 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 AF for inference applications. One such example is known as SCAD. As noted in powerenhancement, the unbiased thresholding is required to avoid size distortions in a large class of high-dimensional testing problems involving a “plug-in" estimator of $\bfsym \Sigma_u$. In particular, this rules out the popular {soft-thresholding} function, which does not satisfy (iii) due to its first-order shrinkage bias.
In empirical asset pricing, factor loadings are known to depend on individual-specific observables $\bfm 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, CL07 and CMO model explicitly the loading matrix as a function of covariates $\bfm X_i$. fan2016projected extended the model to allowing components in factor loadings that are not explainable by characteristics:
Here $\bfm 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 $\{\bfm X_{i}\}_{i=1}^N$ (cross-sectional smoothing), and apply PCA to the smoothed data (fitted values). Let $\{\phi_j(\bfm 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 fan2020statistical. Set $\phi(\bfm X_i)'=(\phi_1(\bfm X_{i}),\cdots.,\phi_J(\bfm X_{i}))$ and $\Phi(\bfm X)=(\phi(\bfm X_1),\cdots,\phi(\bfm X_N))'$, an $N\times J$ matrix. Then the projection matrix on characteristics can be taken as $ \bfm P=\Phi(\bfm X)(\Phi(\bfm X)'\Phi(\bfm X))^{-1}\Phi(\bfm X)'. $ The projected data $\bfm P \bfm Y$ is the fitted value of regressing $\bfm Y$ on to the basis functions.
We make the following key assumptions:
The above conditions require that the strengths of the loading matrix should remain strong after the projection. Condition (ii) implies that if we apply $\bfm P$ to both sides of $\bfm Y=\bfm B\bfm F'+\bfm U$, then $$ \bfm P\bfm Y\approx\bfm P\bfm B\bfm F' = \bfm G \bfm F' $$ where $\bfm G=\bfm P\bfm B$ is the $N\times r$ matrix, which $\approx (\bfm g(\bfm X_i))_{N\times r}$ under additional assumption $ \operatorname{\mathbb E}(\bfsym \gamma_i|\bfm X_i)=0$ for all $i\leq N$. In other words, the noise $\bfm U$ is suppressed, while signals remain. Hence, the scaled sample covariance $(\bfm P\bfm Y)'\bfm P\bfm Y = \bfm Y'\bfm P\bfm Y\approx\bfm F\bfm G'\bfm G\bfm F'. $ For identification purpose, let us assume $\ensuremath{\boldsymbol{\Xi}}:=\bfm G'\bfm G$ is a diagonal matrix and $\bfm F'\bfm F/T=\bfm I$. Then from $$ \frac{1}{T} \bfm Y'\bfm P\bfm Y\bfm F\approx \bfm F\ensuremath{\boldsymbol{\Xi}}, $$ we infer that the columns of $\bfm F$ are approximately the eigenvectors of the $ \bfm Y'\bfm P\bfm Y$, scaled by a factor $\sqrt{T}$. This motivates estimating factors by using the top $R$ eigenvectors of $ \bfm Y'\bfm P\bfm Y$.
fan2016projected 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 $\bfm f_t$, we do know that factors are related to some proxy variables $\bfm W_t$. For example, the latent factors are unknown for equity markets, but they are related to Fama-French factors fama2015five; latent factors for disaggregated macroeconomics time series are unknown, but they are related to aggregated ones mccracken2016fred. Switching the roles of rows and columns, longitudinal regression of each series $\{y_{it}\}_{t=1}^T$ on $\{\bfm W_t\}_{t=1}^T$ yields the projected data matrix, from which latent factors and loadings can be extracted similarly. See fan2020augmented for details on how latent factor learning is augmented by instruments $\{\bfm W_t\}_{t=1}^T$.
In this section, we continue denoting by $R$ as the number of factors we use, and by $r$ as the true number of factors. fan2019learning proposed a simpler factor estimator that does not rely on eigenvectors, by using cross-sectional diversified projections (DP). Let $ \bfm W=(\bfm w_1,\cdots,\bfm w_R) $ be a given exogenous (or deterministic) $N\times R$ matrix, where each of its $R$ columns $\bfm w_k$ is an $N\times 1$ vector of “diversified weights", whose definition is to be clear below. We estimate $\bfm f_t$ by simply taking $$ \widehat\bfm f_t=\frac{1}{N}\bfm W'\bfm y_t. $$ By substituting $\bfm y_t=\bfm B\bfm f_t+\bfm u_t$ into the definition, immediately we have
Thus $\widehat\bfm f_t$ (consistently) estimates $\bfm f_t$ up to an $R\times r$ affine transform $\bfm H$, with the estimation error $\bfm e_t:=\frac{1}{N}\bfm W'\bfm u_t$. The assumption that $\bfm W$ should be diversified ensures that as $N\to\infty$, $\bfm e_t$ is “diversified away" (converging to zero in probability). More specifically, we impose the following assumption.
Conditions (i) and (ii) define the “diversified weights" $\bfm W$. When $(u_{1t},\cdots,u_{Nt})$ are cross-sectionally weakly dependent, they ensure that $ \bfm e_t$ is diversified away. Condition (iii) of Assumption (ref) is a key condition, which requires that $\bfm W$ should not diversify away the factor components in the time series. Several choices of $\bfm W$ can be recommended to satisfy this condition. For instance, if factor loadings satisfy ((ref)), then fix $R$ components of sieve basis functions: $(\phi_1(\cdot),\cdots,\phi_R(\cdot))$, we can define $$\bfm W:= (w_{i,k})_{N\times R},\quad \text{ where } w_{i,k}=\phi_k(\bfm X_i).$$ Alternatively, we can also use transformations of the initial observation $\bfm x_{t}$ for $t=0$, which was considered by juodis2020linear. If $\bfm y_0$ is independent of $\{\bfm u_t: t\geq 1\}$, we can apply $ w_{i,k} = \phi_k(y_{i,0})$ . These weights are correlated with $\bfm B$ through $\bfm y_0=\bfm B\bfm f_0+\bfm 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, barigozzi2018consistent 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$.
fan2019learning 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, karabiyik2019cce applied DP to the context of panel data models in the presence of common factors.
To apply either the PCA or the MLE to estimate the model, we need an initial covariance estimator $\bfm S_y$, whose application requires elements of $\bfm 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 ludvigson2010factor 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 $$ \widetilde y_{it}:=\mbox{sgn}(y_{it})\min\{|y_{it}|, \tau_i\} $$ with predetermined $\tau_i>0$. We then estimate $\mathbb Ey_{it}$ using the truncated-mean $\widetilde y_i:=\frac{1}{T}\sum_{t=1}^T\widetilde y_{it}$. Theorem (ref) shows that the high-dimensional means can be estimated uniformly well if $\mathbb E|y_{it}|^q<M$ for some $q\geq 2$.
catoni2012challenging constructed a robust M-estimator that shares the same Gaussian concentration. fan2017estimation, fan2019robust used the adaptive Huber's loss to define the mean estimator: $$ \widehat y_i=\arg\min_{\mu} \sum_{t=1}^T\psi_{\tau_i}(y_{it}-\mu) $$ where $\tau_i$ is a growing sequence, and $$ \psi_{\tau}(z)=
$$ The following theorem shows that $ \widehat y_i$ also estimates $\mathbb E y_{it}$ well provided that $\max_i\mathbb E y_{it}^2$ is bounded.
The robust mean estimation also applies to estimating covariance as its $(i,j)$ element is of form $\mathbb E y_{it} y_{jt}$. When the high-dimensional data have heavy-tailed components, we can replace the sample covariance by its robust version $\widehat \bfm S_y$ before estimating the factors. By the Gaussian concentration inequality, the robustly estimated covariance $\widehat \bfm S_y$ satisfies $$ \| \widehat \bfm S_{y}- \bfsym \Sigma_y\|_\infty=O_P\left(\sqrt{\frac{\log N}{T}}\right), $$ so long as $\operatorname{\mathbb E} y_{it}^2y_{jt}^2$ is uniformly bounded (and serial independence is assumed).
Based on the above robust covariance inputs, we can create factor estimators and derive their theoretical properties following the guidance of Section (ref). See Chapter 10 of fan2020statistical for further generalizations.
When factors are highly persistent but $\mathbb E \bfm u_t \bfm u_{t-h}^T = 0$, then the cross-covariance $$ \bfsym \Sigma_h = \mathbb E \bfm y_t \bfm y_{t-h}' = \bfm B (\mathbb E \bfm f_t \bfm f_{t-h}) \bfm B', \qquad h \geq 1 $$ contains valuable information about $\bfm B$. This motivates to estimate loadings by applying PCA to aggregated $\{\bfsym \Sigma_h: h=1,\cdots\}$, and we studied by lam2012factor. A related idea has been extended to matrix-variate PCA wang2019factor,chen2020constrained. fan2018optimal also provided a procedure to efficiently aggregate the cross-covariance information with the covariance information when $h=0$.
Many references have documented the comparisons among various estimation methods. westerlund2013estimation 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.
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.
SW02, bai2006confidence considered factor-augmented regression model for $h$-step ahead forecast:
Here $\bfm w_t$ in ((ref)) is the observed predictors, which may include lagged dependent variables. Equation ((ref)) is a high-dimensional factor model that includes a vector of latent factors $\bfm f_t$. The forecast can be implemented by regressing $y_{t+h}$ onto $\bfm w_t$ and estimated factors. The factor model ((ref)) serves as an important dimension reduction tool.
fan2015sufficient generalized ((ref)) to the nonlinear model with multi-indices. Consider the following forecasting model:
where $h(\cdot)$ is an unknown link function, and $\bfsym \varepsilon_{t+1}$ is the error independent of $\bfm f_t$ and $\bfm u_{t}$. Vectors $\bfsym \phi_1, \dots, \bfsym \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 ((ref)) is given by li1991sliced, which shows that under some regularity conditions such as $\bfm f_t$ is elliptically symmetric, we have
for a $R$-dimensional vector $\bfm a( y_{t+1})$, where $\bfsym \Phi=[\bfsym \phi_1, \bfsym \phi_2, \dots, \bfsym \phi_R]$ is an $r\times R$ matrix. In other words, the “inverse regression vector" $ \mathbb E(\bfm f_t| y_{t+1})$ falls in the column space spanned by $\bfsym \Phi$, which can be extracted by PCA. Indeed, since $ \mathbb E( \mathbb E(\bfm f_t| y_{t+1}))= \mathbb E(\bfm f_t)=0$, \[\operatorname{cov}( \mathbb E(\bfm f_t| y_{t+1}))=\bfsym \Phi \mathbb E[\bfm a( y_{t+1})\bfm a( y_{t+1})']\bfsym \Phi'\] The above matrix has $R$ nonvanishing eigenvalues if $ \mathbb E[\bfm a( y_{t+1})\bfm a( y_{t+1})']$ is non-degenerate. Their corresponding eigenvectors have the same linear span as $\bfsym \phi_1, \dots, \bfsym \phi_R$ do. If one can consistently estimate $\operatorname{cov}( \mathbb E(\bfm f_t| y_{t+1}))$, then the subspace spanned by $\bfsym \phi_1, \dots, \bfsym \phi_R$, which is of our primary interests, can be obtained by extracting the top $R$ eigenvectors of the estimated covariance matrix that correspond to the $R$ largest eigenvalues.
However, it is not an easy task to directly estimate the covariance of $ \mathbb E(\bfm f_t| y_{t+1})$. li1991sliced suggested the sliced covariance estimate, a widely used technique for dimension reductions: The sliced covariance matrix also satisfies the fundamental property ((ref)), namely $E(\bfm f_t | y_{t+1} \in \bfm I_k)$ falls in the column space spanned by $\bfsym \Phi$ for any given partition of the range of $ y_{t+1}$ into $H$ “slices" $\bfm I_1, \bfm I_2, \dots, \bfm I_H$. Correspondingly, let
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.
Implementing the above algorithm requires the number of slices $H$, the number of predictive indices $R$, and the number of factors $r$. In practice, $H$ has little influence on the estimated directions, as pointed out in li1991sliced and explained above that property ((ref)) holds. As regard to the choice of $R$, the first $R$ eigenvalues of $\operatorname{cov}(\mathbb E(\bfm f_t| y_{t+1}))$ must be significantly different from zero compared to the estimation error. Several methods such as li1991sliced and schott1994 have been proposed to determine $R$. For instance, the average of the smallest $r-L$ eigenvalues would follow $\chi^2$ distribution if the underlying factors are normally distributed. The number of factors can be determined by a number of methods.
Consider a high-dimensional regression model
where $\bfm g_t$ is a treatment variable whose effect $\bfsym \beta$ is of the main interest. The model contains high-dimensional exogenous control variables $\bfm 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 $(\ensuremath{\boldsymbol{\nu}}, \ensuremath{\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 $\bfm x_t$. Instead, if we substitute ((ref)) to ((ref)), we reach a factor-adjusted regression model:
where $\bfsym \alpha_g'=\ensuremath{\boldsymbol{\theta}}'\bfm B$, $\bfsym \alpha_y'=\bfsym \beta \bfsym \alpha_g'+\ensuremath{\boldsymbol{\nu}}'\bfm B$, and $\bfsym \gamma'=\bfsym \beta \ensuremath{\boldsymbol{\theta}}'+ \ensuremath{\boldsymbol{\nu}}'$. Here $(\bfsym \alpha_y, \bfsym \alpha_g,\bfsym \beta)$ are low -dimensional coefficient vectors while $(\bfsym \gamma, \ensuremath{\boldsymbol{\theta}})$ are high-dimensional sparse vectors. Importantly, the model contains high-dimensional latent controls $\bfm u_t$, which are weakly dependent due to the nature of idiosyncratic noises. The use of $\bfm u_t$ instead of $\bfm x_t$ validates conditions for many high-dimensional variable selection methods.
fan2020factor and hansen2018fac showed that the penalized regression can be successfully applied to ((ref)) to select components in $\bfm u_t$, which are cross-sectionally weakly correlated. Motivated by belloni2014inference, the algorithm can be summarized as follows. For notational simplicity, we focus on the univariate case $\dim(\bfsym \beta)=1$.
Note that $\bfsym \gamma:\to P_\tau(\bfsym \gamma)$ is a sparse-induced penalty function with a tuning parameter $\tau$. When $\ensuremath{\boldsymbol{\theta}}$ and $\bfsym \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 $\ensuremath{\boldsymbol{\varepsilon}}_{g,t}$ and $\eta_t\ensuremath{\boldsymbol{\varepsilon}}_{g,t}$.
More recently, fan2019learning 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 ((ref)) 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 $\bfm x_t=\bfm 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 $\bfm x_t$. This setting is empirically relevant as it allows to avoid pre-testing the presence of common factors for inference.
Figure (ref), taken from fan2019learning, plots the histograms of the t-statistics based on estimated $\bfsym \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 $\bfm W$ is $ (x_{it}, x_{it}^2,\cdots,x_{it}^R)$ at $t=0$. The “double selection" is the algorithm used in belloni2014inference that directly selecting among $\bfm 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.
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 fan2019farmtest and references therein, and Barras:2010fh,Harvey2016, harveyliu18, giglio2019thousands for applications in empirical asset pricing.
Suppose we observe realizations of a random vector $\{\bfm y_t=(y_{1t},\cdots,y_{Nt})'\}_{t=1}^T$. Let $\bfsym \alpha=(\alpha_1,\cdots,\alpha_N)'$ denote its mean vector. We are interested in testing individual hypotheses: $$ H_{0}^i: \alpha_i=0,\quad i=1,\cdots, N. $$
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 \[ \text{FDP}(x)= \frac{\mathcal{F}(x)}{\max\{\mathcal{V}(x),1\}},\quad \text{FDR}(x)=\mathbb E \{\text{FDP}(x)\}. \] 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.
One of the most popular procedures, proposed by benjamini1995controlling, 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 $$x^*= \max\{p_{(i)}: p_{(i)}\leq \tau i/N \}. $$ This method fits into Algorithm (ref) 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 storey2002direct, fan2012estimating, 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 (ref), romano2007control,romano2008control provided alternative procedures for FDR control.
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 $\bfm y_t $ is generated from a latent factor model:
where $\operatorname{\mathbb E} \bfm f_t=0$, and $\bfsym \alpha $ is the mean vector. In empirical asset pricing, the model can be used to identify nonzero alphas out of a large number of assets, and has been studied to identify skilled mutual fund managers, e.g., Barras:2010fh and Harvey2016. The presence of latent factors, however, leads to strong dependence among the t-statistics based on the naive sample means of $\bfm y_t$, which invalidates the weak dependence assumptions. As well documented in the literature, strong dependence creates fundamental challenges to multiple testing, including large standard errors among the estimated $\alpha_i$, unstable FDP's, and conservativeness of the test procedure. Learning dependence $\bfm B \bfm f_t$ and removing it from the model ((ref)) make the data not only weakly dependent but also less noisy (from $\bfm B \bfm f_t+\bfm u_t$ to $\bfm u_t$). This is the basic idea in factor-adjusted robust multiple tests (FarmTest) by using factor-adjusted data $\{\bfm y_{t}- \widehat{\bfm B}\widehat{\bfm f_t}\}_{t=1}^T$; see ((ref)). Furthermore, fan2019farmtest makes adjustments so that it is also robust to heavy tailed data.
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 $\bfm u_t$, $\bfm f_t$ and $\bfm 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 $\bfm y_t$, without using factor adjustments; 2) the factor-adjusted estimator based on PCA. We apply the method of benjamini1995controlling for multiple testing, setting $\tau=0.05$.
The top panels of Figure (ref) 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 (ref).
The middle two panels of Figure (ref) 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 (ref) 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.
giglio2019thousands 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: $$ y_{it}=\alpha_i+\bfm b_i'\ensuremath{\boldsymbol{\lambda}}+\bfm b_i'(\bfm f_t-\mathbb E\bfm f_t)+u_{it}. $$ In the model $\bfm f_t$ contains both observable and latent factors. The model allows nontradable observable factors and $\ensuremath{\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:
Because of many negative alphas from unskilled hund managers, the multiple testing problem should be properly formulated as one-sided hypotheses: $$ H_{0}^i: \alpha_i\leq 0,\quad i=1,\cdots, N. $$ 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, giglio2019thousands proposed to first screen off very bad funds, identified as: $$ \mathcal I=\{i\leq N,\widehat\alpha_i/\text{se}(\widehat\alpha_i)<-c_{NT}\} $$ where $c_{NT}>0$ is a slowly growing sequence to ensure sure screening fan2008sure: $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.
The issue of endogeneity is often encountered in real data applications. Consider the following instrumental variable (IV) regression model \[y_t=\bfm w_t'\bfsym \beta^0+\varepsilon_t=\bfm w_{1t}'\bfsym \beta_1+\bfm w_{2t}'\bfsym \beta_2+\varepsilon_t,\] where $\bfm w_{1t}$ is a $k_1$-dimensional vector of exogenous regressors and $\bfm w_{2t}$ is a $k_2$-dimensional vector of endogenous regressors. Meanwhile, we have an $N$-dimensional IV $\bfm x_t$ which admit a factor structure: $$ \bfm x_t=\bfm B\bfm f_t+\bfm u_t. $$
Below we introduce four estimators for $\bfsym \beta^0$, which differ on their choices of the instruments.
Use $\bfm f_t$ as the instruments. Project $\bfm w_{2t}$ on $\bfm f_t$: \[\bfm w_{2t}=\bfsym \phi'\bfm f_t+\bfm v_{t},\quad \mathbb E(\bfm v_t|\bfm f_t)=0\] where $\bfsym \phi$ is a $k_2\times r$ matrix. We need $r\ge k_2$ for identification. Let $\bfm z_t=(\bfm w_{1t}', \bfm f_t')'$ be the set of instruments. As $\bfm f_t$ is unobservable, we replace it with some factor estimator and apply the two stage least squares estimator $\widehat\bfsym \beta_\bfm f$ with the feasible instruments.
bai2010 studied this estimator, and showed that the estimation errors of $\widehat\bfsym \beta_\bfm f$ associated with the generated instruments (factor estimations) have no effect on the limiting variance. When $\bfm u_t$ and $\varepsilon_t$ are uncorrelated, this only requires $(N,T)\to\infty$ regardless of the relative growth rates. When some weak correlations are present but $\|\mathbb E \varepsilon_t \bfm u_{t}\|_1=O(1)$, we would require $\sqrt{T}=o(N)$ to offset the effect of estimating factors.
Use $\bfm x_t$ as the instruments. Project $\bfm w_{2t}$ on $\bfm x_t$:
where $\ensuremath{\boldsymbol{\theta}}$ is a $k_2\times N$ coefficient matrix. This projection motivates the use of $\bfm 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 \[\widehat\bfsym \beta_{\bfm x}=\Big(\bfm W'\bfm X\widehat \bfsym \Sigma_x^{-1}\bfm X'\bfm W\Big)^{-1}\bfm W'\bfm X\widehat \bfsym \Sigma_x^{-1}\bfm X'\bfm Y\] where $\bfm X$ is $T\times N$ matrix of $\bfm x_t$; $\bfm W$ and $\bfm Y$ are matrices of $\bfm w_t$ and $y_t$. Note that $\widehat \bfsym \Sigma_x$ is the estimated covariance of $\bfm 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\bfsym \beta_\bfm f$ with $\widehat\bfsym \beta_{\bfm x}$. bai2010 showed when $\bfm u_t$ and $\bfm w_t$ are uncorrelated, they have the same asymptotic variance, but $\widehat\bfsym \beta_{\bfm x}$ has a $O(\frac NT)$ bias term. So $\widehat\bfsym \beta_{\bfm x}$ is consistent only if $N=o(T)$.
Use selected $\bfm x_t$ as the instruments. We still consider the projection ((ref)), but assume that rows of $\ensuremath{\boldsymbol{\theta}}$ are sparse vectors so that we can apply penalized regression to select among the components of $\bfm x_t$:
where $P_\tau(\ensuremath{\boldsymbol{\theta}}_j)$ is a sparse-induced penalty with tuning $\tau$. Let $ \bfm x_{t, \text{selec}}$ be the vector of selected components corresponding to nonzero components of $\{\widehat\ensuremath{\boldsymbol{\theta}}_j: j\leq \dim(\bfm w_{2t})\}$. belloni2012sparse used $(\bfm w_{1t}, \bfm x_{t, \text{selec}})$ as the instruments to compute $\widehat\bfsym \beta_{\bfm 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 $\bfm x_t$ invalidates the variable selection procedure ((ref)).
Use $\bfm f_t$ and selected $\bfm 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 ((ref)), we obtain $$ \bfm w_{2t}= \ensuremath{\boldsymbol{\delta}}\bfm f_t+\ensuremath{\boldsymbol{\theta}}\bfm u_t+\bfm e_t $$ where $\ensuremath{\boldsymbol{\delta}}=\ensuremath{\boldsymbol{\theta}}\bfm B$. Hence we can carry out variable selections among $\bfm u_t$: $$ (\widehat\ensuremath{\boldsymbol{\delta}}_j, \widehat\ensuremath{\boldsymbol{\theta}}_j)=\arg\min_{\ensuremath{\boldsymbol{\delta}}_j\in\mathbb R^r, \ensuremath{\boldsymbol{\theta}}\in\mathbb R^N} \frac{1}{T}\sum_{t=1}^T(w_{2t,j}- \ensuremath{\boldsymbol{\delta}}_j'\widehat\bfm f_t-\widehat\bfm u_t'\ensuremath{\boldsymbol{\theta}}_j)^2+P_\tau(\ensuremath{\boldsymbol{\theta}}_j),\quad j\leq \dim(\bfm w_{2t}). $$ Let $ \widehat\bfm u_{t, \text{selec}}$ be the vector of selected components corresponding to nonzero components of $\{\widehat\ensuremath{\boldsymbol{\theta}}_j: j\leq \dim(\bfm w_{2t})\}$. We then use $(\bfm w_{1t}, \widehat\bfm f_t, \widehat \bfm u_{t, \text{selec}})$ as the instruments to compute $\widehat\bfsym \beta_{\bfm f, \bfm u}$, the two stage least squares estimator. This method is expected to work well because it marginalizes out the strong factors in $\bfm x_t$, leaving remaining components $\bfm 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=\bfm w_{2t}'\bfsym \beta^0+\varepsilon_t$, with a single endogenous regressor $\bfm w_{2t}$ generated from ((ref)) with $\bfm e_t= \varepsilon_t/2$. Here $\ensuremath{\boldsymbol{\theta}}=(2,1,-1,0...,0)$ and $\bfm x_t$ admits a two-factor structure. Variables $(\varepsilon_t, \bfm f_t, \bfm B,\bfm 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 ((ref)) we set $P_\tau(\ensuremath{\boldsymbol{\theta}})=\tau\|\ensuremath{\boldsymbol{\theta}}\|_1$ with $\tau= 2.2\|\frac{1}{T}\sum_t\bfm x_t\bfm e_t\|_\infty$.
Table (ref) reports the bias and standard deviation of each estimator calculated from 1000 replications. First, using only estimated factors as the instruments ($\widehat\bfsym \beta_\bfm f$) leads to the largest standard error. This is not surprising because it excludes the relevant information from $\bfm u_t$ while the latter is correlated with $\bfm w_{2t}$, so this method is less efficient. Secondly, using $\bfm x_t$ as instruments without variable selection ($\widehat\bfsym \beta_\bfm x$) has the smallest standard deviation, but is severely biased. Finally, the two instrumental selection based estimators ($\widehat\bfsym \beta_{\bfm x,\text{selec}}$ and $\widehat\bfsym \beta_{\bfm f,\bfm u}$) perform favorably and similarly. But $\widehat\bfsym \beta_{\bfm x,\text{selec}}$ is not as stable, as it occasionally selects none of the instruments in our numerical experiments.
Consider the following factor-augmented regression
where $\bfsym \alpha(\bfm L)=\bfsym \alpha_0+\bfsym \alpha_1\bfm L+\dots+\bfsym \alpha_p\bfm L^p$, $\bfsym \gamma(\bfm L)=\bfsym \gamma_0+\bfsym \gamma_1\bfm L+\dots+\bfsym \gamma_q\bfm L^q$ and $\bfsym \beta(\bfm L)=\bfsym \beta_0+\bfsym \beta_1\bfm L+\dots+\bfsym \beta_l\bfm L^l$, all are lag operator polynomials. Suppose that $\bfm w_t$ is a $k$-dimensional vector and $\bfm 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 bai2009boosting in this context.
Boosting is an ensemble meta-algorithm, which sequentially finds a “committee” of base learners and then makes a collective decisions by using a weighted linear combination of all base learners. The first successful and popular boosting algorithm is AdaBoost freund1997. friedman2001 proposes a generic functional gradient descent (FGD) algorithm, which views the boosting as a method for function estimation. If the squared loss function is specified, the FGD algorithm reduces to the $L_2$-Boosting, which is studied in friedman2001 and buhlmann2003. Suppose that $( y_t, \bfm z_t)_{t=1}^T$ are the observed target and predictive regressors over the sample period. The $L_2$-Boosting algorithm for estimating the conditional mean $\mathbb E (y_t|\bfm z_t)$ is given as follows.
One can apply the above $L_2$-Boosting to the factor-augmented predictive regression ((ref)). As seen in Algorithm (ref), one needs to specify the base procedure in step 3. bai2009boosting 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.
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 bai2009boosting for details.
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: $$ y_{t}= \bfm w_t'\bfsym \beta+\bfm w_t'\ensuremath{\boldsymbol{\delta}}1\{q_t>\gamma\}+\varepsilon_t, $$ where $\bfm w_{t}$ and $q_{t}$ are adapted to the filtration $\mathcal{F}_{t-1}$; $(\bfsym \beta, \ensuremath{\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 $\bfm w_{t}^{\prime }(\bfsym \beta +\ensuremath{\boldsymbol{\delta}})$; when $q_t \leq \gamma$, it reduces to $\bfm w_{t}^{\prime }\bfsym \beta$ chan1993consistency, hansen2000sample. 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-Linton and leemyung extended the model to multivariate threshold: $$ y_{t}= \bfm w_t'\bfsym \beta+\bfm w_t'\ensuremath{\boldsymbol{\delta}}1\{\bfsym \gamma'\bfm f_t>0\}+\varepsilon_t, $$ where $\bfm f_t$ is a vector of “factors" and $\bfsym \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: $ \bfm x_t=\bfm B\bfm f_t+\bfm u_t, $ where $\dim(\bfm x_t)=N\to\infty$. So $\bfm f_t$ can be unobserved factors that can be learned from $\bfm x_t.$ For the identification purpose, suppose $\frac{1}{T}\sum_t\bfm f_t\bfm f_t'=\bfm I$ and $\bfm B'\bfm B$ is diagonal, then $\bfsym \gamma$ and $\bfm 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: $$ \min_{\bfsym \beta,\ensuremath{\boldsymbol{\delta}},\bfsym \gamma}\sum_{t=1}^T(y_t-\bfm w_t'\bfsym \beta-\bfm w_t'\ensuremath{\boldsymbol{\delta}}1\{\bfsym \gamma'\widehat\bfm f_t>0\})^2, $$ where $\widehat\bfm f_t$ is the plugged-in PC-estimator of factors. Because the least squares problem is neither convex nor smooth in $\bfsym \gamma$, the computational task is demanding. leemyung recommended using algorithms based on mixed integer optimization (MIO). Introduce integers $d_t:=1\{\bfsym \gamma'\widehat\bfm 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 $\ensuremath{\boldsymbol{\delta}}$. Define $ M_t \equiv \max_{\bfsym \gamma\in\Gamma} | \bfsym \gamma'\widehat\bfm f_t| $, where $\Gamma$ is the parameter space for $\bfsym \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(\bfm w_t)$,
Then, we can apply modern MIO packages (e.g., Gurobi) to solve for the optimal $(\bfsym \beta,\ensuremath{\boldsymbol{\delta}},\bfsym \gamma)$.
Finally, leemyung 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 $(\bfsym \beta,\ensuremath{\boldsymbol{\delta}})$, but would affect both the rate of convergence and the limiting distribution of the estimated $\bfsym \gamma$.
The stochastic block model has been a popular approach to modeling networks (see abbe2017community for a recent review). We observe a graph of $N$ nodes. Let $\bfm A=(a_{ij})\in\mathbb R^{N\times N}$ be the adjancy matrix of edges so that $a_{ij}=1$ if nodes $i$ and $j$ are connected, and $a_{ij}=0$ otherwise. Suppose each node belongs to one of $r$ communities, and the community that node $i$ belongs to is denoted by an unknown $\pi_i\in\{1,\cdots,r\}$. In addition, elements of $\bfm A$ are random variables. Then stochastic block model assumes that $$ P(a_{ij}=1|\pi_i=k, \pi_j = l)= w_{k,l}, $$ where $w_{k,l}$ is an unknown probability. We observe the matrix $\bfm A$ and aim to recover the membership $\pi_i$ and the probabilities $w_{k, l}$ for all $k, l= 1, \cdots, r$.
Let $\bfm e_1,\cdots,\bfm e_r$ denote the canonical basis in $\mathbb R^r$, and $\bfm b_i=\bfm e_k$ where $\theta_i=k$. Then, $\bfm b_i$ indicates the community membership of node $i$, and the membership matrix is $$ \bfm B=(\bfm b_1,\cdots,\bfm b_N)',\quad N\times r, $$ whose rows represent nodes and columns represent communities. Let $\bfm W$ denote the $r\times r$ matrix of $(w_{k,l})$ and let $\bfm L:=\operatorname{\mathbb E}\bfm A$. It can easily be seen that $\bfm L=\bfm B \bfm W \bfm B' $ is a low-rank matrix, whose rank equals $r$, leading to the following low-rank decomposition: $$ \bfm A= \bfm L+ \bfm S,\quad \bfm S= \bfm A-\operatorname{\mathbb E} \bfm A. $$ Therefore, $\bfm A$ has the familiar decomposition ((ref)), with $\bfm L$ being similar to the systematic risk and $\bfm B$ as a low-rank loading matrix. Since the elements in $\bfm S$ are independent with mean-zero (Wigner matrix), the operator norm $\|\bfm S\|$ does not grow too fast, compared to that of $\bfm L$. We can then apply PCA on $\bfm A$ to estimate $\bfm B$. Suppose $r$ is known, then the estimator $\widehat\bfm B$ is defined as $\sqrt{N}$ times the eigenvectors of $\bfm A$, corresponding to the first $r$ eigenvalues.
Theorem (ref) 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 (\bfm W^{1/2}\bfm B'\bfm B\bfm W^{1/2}) = c_i g_N (1+o_P(1))$ for all $i\leq r$, then there is an $r\times r$ matrix $\bfm H$, so that $$ \|\widehat\bfm B-\bfm B\bfm H\|_\infty =O_P( g_N^{-2} N \|\bfm S\| + g_N^{-1} \sqrt{N\log N}). $$
Therefore, elements of a rotated $\bfm B$ can be estimated uniformly well. Moreover, because each community has many nodes belong to, $\bfm B\bfm 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 ke2015homogeneity on the rows of $\widehat\bfm B$ to consistently identify the communities.
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 $$ y_{it} = \bfm b_{i,t}'\bfm f_t+ u_{it} $$ where $\bfm 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 $\bfm b_{i,t}$ using a function of observed characteristics $\bfm z_{i,t-1}$: $$ \bfm b_{i,t}= \bfm b_i(\bfm z_{i,t-1}) $$ where $\bfm 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 $\bfm z_{i,t-1}$ is correctly specified and indeed can fully capture the degree of time-varyingness of the model, then $\bfm 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 $\bfm z_{i,t-1}$ and omitted variable problems. Above all, we refer to gagliardini2019estimation for an excellent review on conditional factor models using this approach, and their applications in empirical asset pricing.
The second approach assumes that factor loadings change smoothly over time. Suppose $\bfm b_i(\cdot)$ is an unknown smooth function, we assume $$ \bfm b_{i,t}= \bfm b_i\left(\frac{t}{T}\right),\quad \forall t\leq T. $$ Then locally, $ \bfm b_{i,t} \approx \bfm 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: $$ y_{it}\approx \bfm b_{i,r}'\bfm f_t+ u_{it},\quad t\in\mathcal B(r). $$ Motivated by this assumption, ang2012testing and ma2020testing tested the market mean-variance efficiency assumption in the case of known factor case. In the unknown factor case, su2017time 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).
Consider a continuous-time factor model $$ d\bfm y_t=\bfsym \alpha_tdt+\bfm B_td\bfm f_t+d\bfm u_t $$ where $\bfm y_t$, $\bfm f_t,\bfm u_t$ are vectors of asset prices, factors and idiosyncratic risks; $\bfsym \alpha_t$ is a drift term. The time-varying loading matrix $\bfm B_t$ is an $N\times r$ matrix that is assumed to be continuous and locally bounded It\^{o} semimartingale of the form: $$ \bfm B_{t}=\bfm B_{0} +\int_0^t\widetilde\bfsym \alpha_sds+\int_0^t\bfsym \sigma_{s}d\bfm W_s, $$ where $\widetilde\bfsym \alpha_s$ and $\bfsym \sigma_s$ are optional processes and locally bounded; $\bfm W_s$ is a Brownian motion. Roughly speaking, by the Burkholder-Davis-Grundy inequality (cf. chapter 2 of jacod2011discretization), $\bfm 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: ait2017using, chen2019five, liao2018uniform,li2019jump, pelger2019large.
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 cai2016structured,bai2019matrix. Besides, the EM algorithm is also a classical approach to dealing with unbalanced panels. We refer to stock2002macroeconomic, su2019factor, zhu2019high for detailed discussions on related issues.
Recall that the covariance matrix of $\bfm y_t$, under the factor model ((ref)), has the following decomposition, $ \bfsym \Sigma_y= \bfm B\operatorname{cov}(\bfm f_t)\bfm B'+\bfsym \Sigma_u, $ where columns of $\bfm B$ are approximately equal to the eigenvectors of $\bfsym \Sigma_y$ corresponding to the first $r$ eigenvalues. As such, let $\widehat\bfsym \Sigma_y$ be an input matrix, serving as an estimator for $\bfsym \Sigma_y$. Then as described in Section (ref), we can estimate the space spanned by $\bfm B$ using the leading eigenvectors of $\widehat\bfsym \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\bfsym \Sigma_y=(\widehat\sigma_{ij})$ with $$ \widehat\sigma_{ij}:=\frac{1}{w_iw_jT}\sum_{t=1}^Ty_{it}y_{jt}x_{it}x_{jt}. $$ In the matrix form, let $\bfm Y$ and $\bfm X$ be the $N\times T$ matrices of $y_{it}$ and $x_{it}$. So we only observe $\bfm Y\circ \bfm X$, where $\circ$ represents the element-wise matrix product, the Hadamard product. Also let $\bfm W $ be the diagonal matrix with $w_i$ being its $i$ th diagonal entry. Then $$\widehat\bfsym \Sigma_y= \frac{1}{T} \bfm Z\bfm Z',\quad \bfm Z:=\bfm W^{-1}\bfm Y\circ \bfm X. $$ Therefore, columns of the loading matrix estimator $\widehat\bfm B$ equal to $\sqrt{N}$ times the top right singular vectors of $\bfm Z$. This method simply replaces the missing entries of $\bfm 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., imbens2015causal). 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}^Tx_{it}$. But in the case of homogeneous missing, that is, $w_1=\cdots=w_N$, the IPW is not needed, because $\bfm W$ equals the identity matrix up to a constant, which does not affect the PCA on $\bfm Y\circ \bfm 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 abbe2017entrywise, su2019factor under the assumption of homogenous missing. su2019factor used this estimator as their initial value for the EM algorithm. xiong2019large allowed heterogenous missing and proved that the estimators are also asymptotically normal (they estimated $w_iw_j$ directly by $\frac{1}{T}\sum_{t=1}^Tx_{it}x_{jt}$). We can also quickly derive the rate of convergence by applying Theorem (ref). 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 (ref).
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., athey2018matrix, bai2017principal,moon2018nuclear,giglio2019thousands.
In the matrix form $\bfm Y=\bfm M+\bfm U$, the goal is to recover the factor component $\bfm M=\bfm B\bfm F'$ when $\bfm 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\bfm M$. negahban2011estimation and koltchinskii2011nuclear 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 $$ \frac{1}{NT}\|\widehat\bfm M-\bfm M\|_F^2=O_P\left(\frac{1}{T}+\frac{1}{N}\right). $$ chen2020noisy certifies further that the convex optimization ((ref)) 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.
Several recent progress in this literature focuses on debiasing the regularized regression in order to have valid confidence intervals, e.g., chen2019inference,xia2019statistical, chernozhukov2019inference. When the missing is homogeneous, $P(x_{it}=1)=p$ for all $(i,t)$, chen2019inference proposed the following simple debiased estimator
where $H_R(\cdot)$ is the best rank $R$ approximation in ((ref)), $\widehat{\bfm M}$ is given by ((ref)), $\widehat{p}$ is the sample proportion of missing data. The idea is very intuitive. Ignoring the weak-dependence between $\widehat{\bfm M}$ and $\bfm X$ and estimating error in $\widehat p$, we have $$ \operatorname{\mathbb E} (\widehat{\bfm M} + \widehat{p}^{-1}(\bfm Y-\widehat{\bfm M})\circ \bfm X) \approx \operatorname{\mathbb E} \widehat{\bfm M} + \operatorname{\mathbb E} (\bfm Y-\widehat{\bfm M}) = \bfm M, $$ which is approximately unbiased. However, the estimator $\widehat{\bfm M} + \widehat{p}^{-1}(\bfm Y-\widehat{\bfm M}\circ \bfm X)$ is no longer of rank $R$, which increases the variances. This leads to use the projection as in ((ref)), which is asymptotically efficient in terms of both rate and pre-constant.
Alternatively, the debiasing can be achieved through the iterative least squares chernozhukov2019inference. Suppose the true number of factors, $r$, is known.
A key technical argument is to ensure that the estimation error in $\widehat\bfm B$ (step 2) has no impact on the factor estimator (step 3); this is achieved by chen2019inference using an “auxiliary leave-one-out" argument.
When the missing probability $P(x_{it}=1)$ varies across $i$, there are two ways to revise the previous algorithm to achieve the asymptotic normality. One way is to replace ((ref)) with a weighted regularization:
where $\widehat\bfm W$ is a diagonal matrix, whose $i$ th diagonal entry equals $\widehat w_i:=\frac{1}{T}\sum_{t=1}^Tx_{it}$. This debiases the least squares part of the loss function, adopting the same idea of inverse probability weighting. The remaining steps of Algorithm (ref) are the same. Then the same “auxiliary leave-one-out" technical argument of chen2019inference still goes through. The other way is to apply “sample splitting", which evenly split the columns of $\bfm Y$ into two parts: on one part we run the penalized regression as in ((ref)) and obtain $\widehat\bfm B$, on the other part we run iterative least squares. Then exchange the two parts and re-do the estimations. The final estimator is taken as the average of the two. Suppose $u_{it}$ is serially independent, the sample splitting then artificially creates independences among various statistics from the splitting sample. See chernozhukov2019inference for detailed descriptions of this approach.
In an attempt to handle endogenous missing, bai2019matrix proposed a block-rearrangement method. At the cost of this generality, they require that the data matrix $\bfm Y$ should have a sufficiently large balanced sub-block after elementary rearrangements. See cai2016structured,fan2019structured for related ideas.
Specifically, a preliminary step of their estimation is to rearrange the data in a shape that all the factor loadings can be estimated in one sub-block and all the factors can be estimated in another sub-block. The following example is adapted from bai2019matrix, which gives a good illustration on this manipulation: example of the $N\times T$ matrix for $y_{it}$: \[
\Longrightarrow
.\] The left matrix is the originally collected data and the right is the rearranged one. The symbols with asterisk denote the missing data. From the column perspective, the 1st, 2nd and 4th columns have missing values and therefore are rearranged as the last three columns in the right panel; from the row perspective, the 2nd, 3rd and 4th rows have missing values and therefore are rearranged as the last three rows in the right panel. bai2019matrix name the black block “bal”, name the black plus the red blocks “tall”, and name the black plus the blue block “wide”.
Consider the missing value $\bfm y_{22}^*$. We want to replace it with its expected value $\mathbb E(\bfm y_{22}^*)=\bfm b_2'\bfm f_2$. Note that $\bfm y_{22}^*$ shares the same factor loadings $\bfm b_2$ with data points $\bfm y_{23}$ and $\bfm y_{25}$ in the wide block; and shares the same factors $\bfm f_2$ with data points $\bfm y_{12}$ and $\bfm y_{52}$ in the tall block. Meanwhile, $\bfm b_2$ can be estimated using data in the“tall" block; $\bfm f_2$ can be estimated using data in the “wide" block. As a result, one might expect to recover $\mathbb E(\bfm y_{22}^*)$ with these two estimators. However, we must take into account the rotational indeterminacy inherent with the factor models. For a generic missing value $\bfm y_{it}$, \[\widehat\bfm b_{\mathrm{tall},i}=\bfm H_{\mathrm{tall}}'\bfm b_i+o_P(1), \quad \widehat\bfm f_{\mathrm{wide},t}=\bfm H_{\mathrm{wide}}^{-1}\bfm f_t+o_P(1).\] Therefore \[\bfm b_i'\bfm f_t=\widehat\bfm b_{\mathrm{tall},i}' \bfm A\widehat\bfm f_{\mathrm{wide},t}+o_P(1),\quad \bfm A:= \bfm H_{\mathrm{tall}}^{-1}\bfm H_{\mathrm{wide}} .\] To estimate $\bfm A$, by $\widehat\bfm f_{\mathrm{wide},t}=\bfm H_{\mathrm{wide}}^{-1}\bfm f_t+o_P(1)$ and $\widehat\bfm f_{\mathrm{tall},t}=\bfm H_{\mathrm{tall}}^{-1}\bfm f_t+o_P(1)$, we have $$\widehat\bfm f_{\mathrm{tall},t}=\bfm A\widehat\bfm f_{\mathrm{wide},t}+o_P(1).$$ So one can run the regression of $\widehat\bfm f_{\mathrm{tall},t}$ on $\widehat\bfm f_{\mathrm{wide},t}$ to consistently estimate $\bfm A$. This leads to the following estimation procedure.
Once $\widetilde\bfm Y$ is obtained, we apply the PCA again to the imputed data $\widetilde\bfm Y$ to get more efficient estimates of $\bfm B$ and $\bfm F$. Suppose the size of the “tall” block is $N\times T_0$ and the size of the “wide” block is $N_0\times T$. So the size of the “bal” block is $N_0\times T_0$. The whole sample size (including missing data points) is $N\times T$. Bai and Ng require that $$ \max\{\sqrt{N}, \sqrt{T}\} =o(N_0), \quad \text{ and } \quad \max\{\sqrt{N}, \sqrt{T}\} =o(T_0). $$ An implication of the above condition is that the missing data points should not be too frequent in the sense that the balanced subblock is large enough. Though this condition rules out the case of random missing (e.g., missing occurs as outcomes of Bernoulli trials), it is not stringent given the nature of endogenous missing.
We conduct a simulation study to compare six matrix completion approaches, namely:
IPW. The inverse probability weighting.
ReUW. Unweighted regularization. The eigenvectors of the estimator ((ref)).
ReW. Weighted regularization. The eigenvectors of the estimator ((ref)).
ReDebias. The debiased regularized estimator from Algorithm (ref).
EM. The EM algorithm.
We generate a two-factor model where loadings, factors and $u_{it}$ are independent standard normal. Under the homogeneous missing we generate $x_{it}\sim$ Bernoulli$(0.5)$; under the heterogeneous missing we generate $x_{it}|w_i\sim$ Bernoulli$(w_i)$, and $w_i\sim $Uniform$[0.1,1]$. The three regularized methods require choosing $\lambda$, the tuning parameter. Write the penalized loss function to be $\|( \bfm W^{-1/2}\bfm Y- \bfm W^{-1/2}\bfm M)\circ\bfm X\|_F^2+\lambda\|\bfm M\|_n$ where $\bfm W$ is a diagonal weighting matrix. The theory requires that with a high probability, there is $c>0$, $$ (2+c)\|\bfm U\circ(\bfm W^{-1}\bfm X)\|<\lambda. $$ So we set $\lambda$ to be the 0.95 quantile of $ 2.2\|\bfm Z\circ(\bfm W^{-1}\bfm X)\|$ where $\bfm Z$ is an $N\times T$ matrix of standard normal variables. In practice, one can also simulate $\bfm Z$ using the estimated idiosyncratic covariance matrix.
We compare the performance of estimating the loading space, measured by $\bfm P_{\bfm B}=\bfm B(\bfm B'\bfm B)^{-1}\bfm B'$. Table (ref) reports $\|\bfm P_{\widehat\bfm B}-\bfm P_\bfm B\|$ averaged over 100 replications for each method. In all scenarios, the IPW performs the worst among all estimators. Under the homogeneous missing, all the other four methods perform similarly, but the difference is much more noticeable under the heterogeneous missing. The general ranking is that $$ \text{IPW} \prec \text{ReUW} \prec \text{ReW} \prec\text{ReDebias} \approx \text{EM}. $$ This ranking is as expected: IPW is the least efficient method among the five; ReUW uses the nuclear-norm regularized estimation that does not take into account the heterogeneous missing or debias; ReW accounts for the heterogeneous missing probabilities, and ReDebias further removes the regularization bias.
Finally, it is not surprising to see that ReDebias and EM perform similarly because both start with an initial low-rank estimator (ReDebias initializes from ReW while EM initializes from IPW), then proceed via iterative least squares. But we note that ReDebias operates much faster because it only iterates once, so is more attractive than EM in handling large scale problems. We also implemented the “early-stop-EM" (which only iterates twice), it performs only slightly better than IPW and is worse than all the other estimators. Therefore we conclude that the ReDebias is a recommended method for handling large scale low-rank matrix completion problems.
We have conducted a selective overview on the recent developments of the factor model and its application on statistical learning. We focus on the perspective of the low-rank structure of factor models, and particularly draws attentions to estimating the model from the low-rank recovery point of view. New estimation and inference methods, and matrix completion problems have been discussed.