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.
88,365 characters · 12 sections · 63 citation commands
Approximate Factor Models with Strongly Correlated Idiosyncratic Errors
\def\spacingset#1{ {#1}} \spacingset{1}
\if11 \fi
\if01 \fi \thispagestyle{empty}
{\it Keywords:} convex optimization; alternating minimization; convergence; high-probability error bounds.
\setcounter{page}{1} \spacingset{1}
Factor models are widely used in a number of scientific fields for reducing the dimension of data sets comprising of a large number of variables anderson2003introduction. A factor model assumes that each variable under consideration can be expressed as a linear combination of a small number of {\em latent} factors plus an idiosyncratic component (error term). Co-movements among the variables can be accounted for by these few factors thus aiding interpretation. When {\em exact} factor models are used in the analysis of cross-sectional data, it is assumed that the idiosyncratic components are mutually uncorrelated anderson2003introduction. However, for time series data such assumptions are often too restrictive, especially if a large panel of time series is considered where the common factors may not fully capture all relationships among the observed time series; in this case, it is of interest to examine {\em approximate} factor models that also allow for correlations amongst the idiosyncratic components.
Such an approximate factor model was introduced in chamberlain1983arbitrage for the analysis of portfolios comprising of a large number of assets. Since then, a number of papers have appeared in the literature investigating properties of such approximate factor models, under the assumption that the correlations between the common factors and the idiosyncratic component, as well as those amongst the idiosyncratic components are {\em weak}. Formally, the approximate factor model is defined as
where $X_t$ is a vector of $p$-dimensional time series, $F_t$ a $K$-dimensional latent factor process, $\Lambda$ a $p\times K$ matrix of {\em factor loadings} and $u_t$ the vector of idiosyncratic components. It is often further assumed that the factor process exhibits Vector Autoregressive dynamics, namely $F_t = \sum_{i=1}^{q} \Phi_i F_{t-i}+\eta_t$, where $\eta_t$ is a serially uncorrelated error process that is independent across its coordinates, and $\Phi_i$ are $K\times K$ transition matrices. The model in (ref) is typically estimated through principal component (PC) decomposition, which operates under the assumption that as the time series panel size $p\rightarrow\infty$, the leading $K$ eigenvalues of $\Sigma_X:=\mathbb{E}(X_tX_t^\top)$ diverge, whereas all eigenvalues of $\Sigma_u:=\mathbb{E}(u_tu_t^\top)$ are bounded, thus enabling the separation between the common factors and the idiosyncratic components. Some key theoretical results for this model are given in bai2002determining,bai2003inferential, where asymptotic normality of the estimated factors and factor loadings\footnote{up to some invertible transformation} obtained from PC analysis is established, under a $\sqrt{p}/T\rightarrow 0$ scaling for the former result and a $\sqrt{T}/p\rightarrow 0$ scaling for the latter. Further, if $T/p\rightarrow 0$, then the maximum time-indexed deviation of the estimated factors relative to their true values vanishes. In later work, stock2005implications consider the same factor model representation, but each coordinate of $u_t$ is allowed to exhibit serial correlation, and assumed to be uncorrelated with $F_t$ across all time leads and lags. By decorrelating the coordinates of the $u_t$ error process\footnote{With a slight abuse of notion, here we use $F_t$ to denote the term that collects the lags of the actual factors that enter the model after decorrelation.}, the model can written in the form of
where $D(L)=\text{diag}(\delta_1(L),\dots,\delta_p(L))$ is a {\em diagonal} matrix with each entry being the autoregressive polynomial corresponding to coordinates of $u_{t}$, while $\epsilon_t$ is a pure noise term that is neither cross-sectionally nor serially correlated.
The presence of strongly correlated idiosyncratic components in the model can lead to distorted estimation and inference, resulting in overestimation of the number of factors greenaway2012estimating and is detrimental for forecasting purposes anderson2007forecasting. Within the DFM framework, a common remedy entails the inclusion of lagged terms of $X_t$ anderson2007forecasting,carare2010spillovers,liu2013modelling,eichengreen2012subprime, which however augments the number of model parameters at the rate of $p^2$. To overcome the technical issues arising from jointly estimating a large number of parameters, these methods either implicitly assume a small panel size anderson2007forecasting so that their application does not suffer from the curse of dimensionality, or resort to estimating $p\choose2$ models based on pairwise univariate series from the $X_t$ components liu2013modelling.
In a related line of work, forni2000generalized introduced the generalized dynamic factor model (GDFM) framework that dictates the existence of two {\em mutually orthogonal} processes that capture the common and idiosyncratic components, respectively. The dynamic factors spanning the common space can be general $L^2$-integrable processes and estimated through principal components in the frequency domain forni2005generalized,forni2015dynamic.
Despite the generality of the GDFM framework, whose formulation ensures orthogonality between the independent and identically distributed noise process and the common space, its common space recovery relies on estimated spectral density matrices that can exhibit numerical instabilities when the dimensionality of $X_t$ becomes large fiecas2014data. On the other hand, the aforementioned shortcomings of the DFM framework can be largely mitigated, through the inclusion of lagged terms and a formulation that {\em jointly} estimates model parameters via a computationally stable procedure. To this end, following stock2005implications, we write the approximate factor model in the form given in (ref), but allow for $D(L)$ to exhibit cross-correlation structure; i.e. $D(L)$ is not restricted to be diagonal, but merely {\em sparse}. Hence, the dynamics of the $p$ time series in $X_t$ can be written in the form of a lag-adjusted static factor model, with the lag term impacting the current values through {\em sparse} transition matrices. Through cross-sectional de-correlation, $\epsilon_t$ becomes a strictly exogenous noise process comprising of independent and identically distributed shocks, and the model representation aligns with that under the GDFM framework, with the lagged term(s) and the factors collectively capturing the common space and accounting for all pervasive shocks. In the proposed model specification, two key quantities are the space spanned by the common factors (the “factor hyperplane" henceforth), and the sparse transition matrices of the time lags of the observable process. To obtain their estimates, we formulate a {\em penalized maximum likelihood} objective function, introduce a block-coordinate descent algorithm to solve the posited optimization problem and establish finite sample high-probability error bounds for the convergent solution estimates. Finally, note that the transition matrices of the lagged values can provide useful and interpretable information, as shown in our application study and noted in eichengreen2012subprime,liu2013modelling.
Together with the proposed model specification that allows for strongly correlated idiosyncratic components in approximate factor models, key contributions of this work entail the convex formulation that leads to joint estimation of the model parameters, as well as the technical developments that provide insights on appropriately handling the interaction between the {\em latent} factor space and the lagged space spanned by the past history of the observed process. In particular, the strategy used to establish error bounds is applicable to other high-dimensional statistical models involving simultaneously observed and latent components.
The remainder of this paper is organized as follows. In Section (ref), we introduce our model setup, estimation procedure for model parameters, as well as steps for performing forecast. Theoretical properties of the proposed estimators are established in Section (ref), including their high-probability statistical error bound, and convergence property. In Section (ref), we introduce an empirical implementation procedure and present the performance evaluation of the estimates based on synthetic data. In Section (ref), an application of our model to weekly stock return data of large US financial institutions for the period from 2001 to 2017 period is considered. Finally, Section (ref) concludes the paper.
{\em Notation.} Throughout this paper, for some generic matrix $A$ of dimension $m\times n$, we use ${\vert\kern-0.25ex\vert\kern-0.25ex\vert \cdot \vert\kern-0.25ex\vert\kern-0.25ex\vert}$ to denote its matrix norms, including the operator norm ${\vert\kern-0.25ex\vert\kern-0.25ex\vert A \vert\kern-0.25ex\vert\kern-0.25ex\vert}_{\text{op}}$, the Frobenius norm ${\vert\kern-0.25ex\vert\kern-0.25ex\vert A \vert\kern-0.25ex\vert\kern-0.25ex\vert}_{\text{F}}$, the nuclear norm ${\vert\kern-0.25ex\vert\kern-0.25ex\vert A \vert\kern-0.25ex\vert\kern-0.25ex\vert}_{*}$, ${\vert\kern-0.25ex\vert\kern-0.25ex\vert A \vert\kern-0.25ex\vert\kern-0.25ex\vert}_1 = \max_{1\leq j\leq n}\sum_{i=1}^m |a_{ij}|$, and ${\vert\kern-0.25ex\vert\kern-0.25ex\vert A \vert\kern-0.25ex\vert\kern-0.25ex\vert}_\infty = \max_{1\leq i\leq m}\sum_{j=1}^n|a_{ij}|$. We use $\|A\|_1=\sum_{i,j}|a_{ij}|$ and $\|A\|_\infty = \max_{i,j}|a_{ij}|$ to denote the elementwise $1$-norm and infinity norm. Additionally, we use $\varrho(A)$ to denote its spectral radius $(\max|\lambda(A)|)$. For two matrices $A$ and $B$ of commensurate dimensions, denote their inner product by $\savebox{\@brx}{\(\m@th{\langle}\)} \mathopen{\copy\@brx\mkern2mu\kern-0.9\wd\@brx\usebox{\@brx}}A, B\savebox{\@brx}{\(\m@th{\rangle}\)} \mathclose{\copy\@brx\mkern2mu\kern-0.9\wd\@brx\usebox{\@brx}}= \text{trace}(A^\top B)$. Finally, we write $A\succsim B$ if there exists some absolute constant $c$ that is independent of the model parameters such that $A\geq cB$.
We start by introducing the model assuming that the idiosyncratic component follows the aforementioned sparse $\mathrm{VAR}(d)$ model, which simultaneously incorporates the cross-sectional and serial structure among its coordinates. To convey the main arguments, we assume without loss of generality that $d=1$ for the ease of exposition, and present the extension to the general lag case in the Supplement.
The starting point is the dynamic factor representation of the observable process $X_t = \tilde{\lambda}(L)f_t + u_t$, where $f_t$ is the common latent factor; $u_t$ is the idiosyncratic component whose dynamics satisfy $\mathcal{B}(L)u_t = \epsilon_t$ with $\mathcal{B}(L)=\mathrm{I}_p-BL$ being the lagged matrix polynomial for some {\em weakly sparse} $B$. Multiplying $\mathcal{B}(L)$ on both sides leads to the dynamic factor model consisting of (ref) and (ref), where $F_t\in\mathbb{R}^K (K\ll p)$ collects the lags of $f_t$ so that it only enters the dynamics of $X_t$ contemporaneously, and is additionally assumed to follow some VAR model with lagged polynomial $\Phi(L)$:
$\epsilon_t$ is a mean zero noise process that is both serially and cross-sectionally uncorrelated. Moreover, it is strictly exogenous satisfying $\text{Cov}(X_{t-1},\epsilon_{t+h})=0$ and $\text{Cov}(F_t,\epsilon_{t+h})=0$, $\forall\,h\geq 0$. The parameters of interest are the factor hyperplane (to be specified later) and the sparse transition matrix $B$. Note that here we only require $B$ to be {\em weakly sparse}, the notion of which can be formalized through the definition of an $\ell_q$ ball with radius $R_q$ negahban2012unified:
The case of exact sparsity corresponds to $q=0$ where $B\in\mathbb{B}_q(R_0)$ has at most $R_0$ nonzero entries; whereas for $q\in(0,1]$, the $R_q$ ball imposes constraints on the decay rate of $|B_{ij}|$'s.
To ensure that $X_t$ is covariance stationary, we require that the spectral radius of $B$ satisfies $\varrho(B)<1$ without further restricting $\Lambda$. Additionally, note that under the assumption that the spectral density of $X_t$ exists, the spectral density of the filtered process $Z_t := \mathcal{B}(L) X_t = \Lambda F_t + \epsilon$ satisfies
Correspondingly, the spectral density of $X_t$ is given by
where $g_X(\omega)$ and $g_{X,Y}(\omega)$ respectively denote the spectrum and cross-spectrum of some generic process $\{X_t\}$ and $\{Y_t\}$:
with $\Gamma_X(h) :=\mathbb{E}(X_tX_{t-h}^\top)$ and $\Gamma_{X,Y}(h) :=\mathbb{E}(X_tY_{t-h}^\top)$.
Given a sample of the $p$-dimensional observable process $X_t$, denoted by $\{x_0,x_1,\dots,x_T\}$, let
where $\mathbf{X}_{T}\in\mathbb{R}^{T\times p}$ and $\mathbf{X}_{T-1}\in\mathbb{R}^{T\times p}$ respectively denote the contemporaneous response matrix and the lagged predictor matrix, and $\mathbf{F}\in\mathbb{R}^{T\times K}$ denotes the latent factor matrix with the latent factor $F_t$ at time point $t$ stacked in its rows. The noise matrix $\mathbf{E}$ is analogously defined. We additionally define the {\em factor hyperplane} associated with the latent factor $\mathbf{F}$ as $\Theta:=\mathbf{F}\Lambda^\top \in \mathbb{R}^{T\times p}$, and note that $\Theta$ has rank at most $K$. With the above notations, the model in (ref) for the observed samples can be written as $\mathbf{X}_{T} = \Theta + \mathbf{X}_{T-1}B^\top + \mathbf{E}$. With the transition matrix $B\in\mathbb{R}^{p\times p}$ assumed sparse and the factor hyperplane $\Theta=\mathbf{F}\Lambda^\top\in\mathbb{R}^{T\times p}$ being low rank, we formulate the following constrained optimization problem:
with the feasible region determined through a rank constraint imposed on $\Theta$ and a sparsity-inducing norm constraint imposed on $B$.
The rank constraint in (ref) leads to a {\em non-convex} feasible region, making it particularly hard to characterize the obtained solution analytically as it depends on the initial values provided to the algorithm. Thus, as commonly undertaken in the literature agarwal2012noisy, we consider a tight convex relaxation of the rank constraint, and the solution to the convexified program has convergence guarantees independent of the initializer. Formally, we consider obtaining the estimator through the convex program in (ref), which can be obtained from (ref) by alternatively considering the nuclear norm constraint for the factor hyperplane and the $\ell_1$ norm constraint for the sparse transition matrix $B$ in Lagrangian form: {
} where $\lambda_B$ and $\lambda_\Theta$ are tuning parameters. The solution $(\widehat{B},\widehat{\Theta})$ can be obtained by a block-coordinate descent algorithm which alternately minimizes with respect to $B$ and $\Theta$, as outlined in Algorithm (ref).
\paragraph{Reconstruction of the factors.} The solution to (ref) provides an estimate of the factor hyperplane, based on which realizations of the $K$-dimensional latent factors process can be reconstructed under certain identifiability restrictions. As mentioned in Section (ref), for any invertible matrix $R\in\mathbb{R}^{K\times K}$, the following equality holds
hence, given a factor hyperplane and the latency of the factors, to fully identify the factors and the corresponding loading matrix $(\mathbf{F},\Lambda)$ from their observationally equivalent counterpart $(\check{\mathbf{F}},\check{\Lambda})$, a total number of $K^2$ restrictions is required to address their indeterminacy. Various choices for the identification restrictions have been discussed in the literature bai2008large, including the most popular PC estimator stock2002forecasting which assumes orthogonality for both the factors and the loadings, as well as the ones that implicitly assume certain ordering of the factors and impose specific structural restrictions on the loading matrix bai2013principal. Under these restrictions, the factors and the loading matrix can always be uniquely identified\footnote{For the PC estimator or under the PC2 restriction, where $\mathbf{F}'\mathbf{F}/n=\mathrm{I}_K$ and $\Lambda$ is assumed lower-triangular, the identification is up to sign rotation; under the PC3 one, where the upper $K\times K$ upper sub-matrix of $\Lambda$ is assumed an identity matrix and $\mathbf{F}$ is left unrestricted, the identification is exact bai2013principal.} and obtained based on the SVD of the estimated hyperplane $\widehat{\Theta}$. It is worth noting that regardless of the identification restrictions that lead to different versions of the estimated factors, the space spanned by the estimated factors is invariant once $\widehat{\Theta}$ is obtained. Specifically, forecasting future values of $x_t$ does not require an exact recovery of $F_t$, as discussed next.
Given estimates $\widehat{B}$ of the transition matrix and $\widehat{\Theta}$ of the hyperplane, we consider the following procedure that first obtains forecasts of the filtered process $Z_t:=X_t - BX_{t-1}$ through projection onto the factor space, followed by a lag adjustment to obtain those of the $X_t$.
To this end, according to the model in (ref), the filtered process $Z_t$ can be represented as $Z_t = \Lambda F_t + \epsilon_t$, whose $h$-step-ahead best linear predictor based on $F_{T-k},k\geq 0$ is given by the projection $\text{Proj}(Z_{T+h}\,|\,\text{Span}(\mathbf{F},T))$, where $\text{Span}(\mathbf{F},T)$ denotes the linear space spanned by $\{F_t\}_{t=1}^T$ stock2002forecasting,forni2005generalized. In particular, based on estimate $\widehat{B}$, the filtered process $Z_t$ can be estimated through $\widehat{z}_t:= x_t - \widehat{B}x_{t-1}$, whose common space estimate corresponds to $\widehat{\Theta}$. Using the surrogate process $\{\widehat{z}_t\}$, let the sample covariance be $\widehat{\Gamma}_Z(h):=\tfrac{1}{T-h}\sum_{t=h+1}^T \widehat{z}_t \widehat{z}^\top_{t-h}$; the $h$-step-ahead forecast of $\{z_{t}\}$ is then given by
where columns of $\widehat{V}$ are the right singular vectors of $\widehat{\Theta}$ corresponding to nonzero singular values, and are effectively an orthonormal basis for the factor space. In the case where $h=1$, $\widehat{x}_{T+1} = Bx_T + \widehat{z}_{T+1|T}$; in the case where $h>1$, $\widehat{x}_{T+h}$ can be obtained inductively by sequentially estimating $x_{T+i}$, for all $0<i\leq h$. Algorithm (ref) outlines the forecasting procedure.
\paragraph{Connections to GDFM.} To conclude this section, we discuss similarities of the proposed formulation to the GDFM forni2000generalized,forni2005generalized. GDFM encompasses a broader class of factor models wherein the observed process admits a decomposition into two mutually orthogonal processes that respectively capture the common and the idiosyncratic component forni2000generalized, with the former not limited to a VAR representation. From a modeling perspective, as pointed out in lutkepohl2014structural, the distinction between GDFM and the state-space form of DFM stock2005implications is not that substantial, since stationary processes can be approximated arbitrarily well by VAR processes with unrestricted order of lags. From the estimation perspective, however, to accommodate the potentially more complex dynamics of the factor processes, GDFM recovers the common space leveraging the spectral domain features of the processes, and the dynamic factors are obtained by solving a generalized eigen-equation w.r.t. estimated covariance matrices of the common and the idiosyncratic components forni2005generalized. In our formulation, we posit an optimization that jointly estimates the factor space and the lagged space, which accounts for the lag information explicitly, under mild sparsity assumptions.
The empirical performance of the two procedures is considered and compared in Section (ref) under various data generating mechanisms.
To establish statistical properties of the estimators, a ball constraint on the feasible region of $\Theta$ is required to incur additional compactness on the low rank component that limits the spikiness of its entries, and this enables identification of the sparse component $B$. To this end, throughout this section, we consider estimators that are solutions to the following convex program:
where $\mathbb{B}_\infty(\phi,\mathbf{X}_{T-1})$ is a box constraint given by
$\phi$ is chosen such that the true value of the parameters $\Theta^\star$ is always feasible. We will provide further illustration on the interpretation of such a box constraint in Section (ref) and Remark (ref). $(\widehat{B},\widehat{\Theta})$ falls into the class of {\em regularized $M$-estimators}, whose theoretical properties have been extensively studied in the statistical literature for diverse settings agarwal2012fast,loh2012high.
A road map to establish properties of the estimators for $(\Theta, B)$ is given next: first in Section (ref) we derive non-asymptotic statistical error bounds of $\widehat{\Theta}$ and $\widehat{B}$ under certain regularity conditions, when the proposed estimation procedure is based on a {\em deterministic realization} of the observable process $\{X_t\}$. In particular, the required regularity conditions primarily entail the {\em restricted strong convexity} (RSC) condition agarwal2012noisy and that the choice of $\lambda_B$ and $\lambda_\Theta$ is in accordance with some {\em deviation condition} loh2012high. Subsequently, in Section (ref), we establish that the required conditions are satisfied with high probability, and provide probabilistic analogues of key model parameters' error bounds for {\em random realizations} drawn from the underlying observable Gaussian process $\{X_t\}$ and the latent process $\{F_t\}$. We also briefly discuss how the model identifiability issue is tackled through the constrained formulation adopted in (ref). Finally in Section (ref), from a numerical perspective, we establish the convergence of the proposed iterative algorithm to a stationary point. All proofs are deferred to Appendices (ref) and (ref). Throughout our exposition, we use superscript $\star$ to denote the true value of the parameters of interest, and denote the errors of the estimators by $\Delta_{\Theta}:=\widehat{\Theta}-\Theta^\star$ and $\Delta_B:=\widehat{B}-B^\star$, respectively.
We start by introducing some additional notation needed in the ensuing technical developments. Let $\ell_T(B,\Theta;\mathbf{X})$ denote the loss function, given by
The true number of latent factors is given by $K$ and thus $\text{rank}(\Theta^\star)=K$. Further, given some $\eta>0$ (to be specified later), let $S^\star_\eta$ denote the thresholded support set of $B^\star$, and we use $s_\eta$ to denote its cardinality, that is, $S^\star_\eta := \{(i,j)\,|\, |B^\star_{ij}|>\eta \}$ and $s_\eta := \|S^\star_\eta\|_0$. Finally, let $S_\mathbf{E} : = \mathbf{E}^\top_{T-1}\mathbf{E}_{n-1}/T$ denote the sample covariance matrix of the error process and let $\Lambda_{\max}(S_\mathbf{E})$ be its maximum eigenvalue. Formally, the {\em RSC condition} agarwal2012noisy,negahban2012unified is defined as follows.
Further, for high dimensional sparse VAR models ($\Theta=0$ in the current setup), the tuning parameter $\lambda_B$ needs to satisfy a {\em deviation condition} loh2012high,basu2015estimation, namely,
which can be simplified to $\lambda_B\geq \|\mathbf{X}_{T-1}^\top\mathbf{E}/T\|_\infty$. Under the current model setup, however, the deviation condition is significantly more involved and requires proper modifications to incorporate quantities associated with the factor hyperplane, as seen in Theorem (ref).
Before stating the main results, we provide a brief discussion on the box constraint on $\Theta$, which aims to “limit" the spikiness of the low rank component, and hence the interaction between the latent factor space and the observable lag space spanned by $X_{t-1}$ --- in particular, for $\Theta$ and $B$ to be properly recovered, such interaction can not be too large. Due to the basis vectors of the factor space being latent, a direct restriction on the interaction is impractical and conceptually unsatisfying, whereas the box constraint adopted effectively restricts the product of the signals from the two spaces and serves our objective, as shown in the proof of Theorem (ref) and Remark (ref). Note that this constraint is in the same spirit to similar ones in the literature agarwal2012noisy,negahban2012restricted, and the norm of $\mathbf{X}$ is necessary since the two spaces have distinct bases.
Next, we comment on the error bound in (ref) and the required conditions in (ref). The error bound encompasses three terms that are respectively associated with the transition matrix $B$, the low rank factor space $\Theta$, and the tolerance $\tau_T$ which measures the extent to which the log-likelihood function deviates from strong convexity (see Definition (ref)). Both $\mathcal{E}_B$ and $\mathcal{E}_\Theta$ depend on three components: (1) the overall curvature of the log-likelihood function as captured by $\alpha_{\text{RSC}}$, (2) the interaction structure between various components of the underlying process, as captured by the tuning parameters $\lambda_B$ and $\lambda_\Theta$, and (3) the inherent structure of the parameters as captured by $s_\eta$, $\sum_{(i,j)\notin S^{\star}_{\eta}}|B_{ij}|$ and $K$ --- in particular, due to the approximately sparse structure of $B^\star$, both the density level $s_\eta$ of its strong support set and the magnitude of its “weak" entries play a role, with the two respectively reflecting the {\em estimation error} and the {\em approximation error} agarwal2012noisy after proper scaling. The curvature as measured by $\alpha_{\text{RSC}}$ dictates the constraint to which the tolerance $\tau_T$ needs to conform (see Equation (ref)), and such a constraint is also interrelated to $s_\eta$ and $K$: for (ref) to be satisfied, neither $K$ nor $s_\eta$ can be too large. Moving to the tuning parameters, $\lambda_B$ can be sub-divided into two terms: the cross-product term $\|\mathbf{X}_{T-1}^\top \mathbf{E}/T\|_\infty$ measures the maximum interaction between the design matrix $\mathbf{X}_{T-1}$ and the noise $\mathbf{E}$, which according to model assumption (population level) should center around 0; where the term $\phi/\sqrt{Tp}$ corresponds to an upper bound on the interaction between the latent (factor) space and the observed one ($\mathbf{X}_{T-1}$). For $\lambda_\Theta$, we require that it dominates the maximum signal coming from the error process in the form of $\Lambda^{1/2}_{\max}(S_{\mathbf{E}})$. Thus, a smaller $\lambda_B$ is needed when interactions between associated terms are weaker and similarly a smaller $\lambda_\Theta$ is needed if the magnitude of the noise is weaker, thus leading to a tighter error bound for the estimates. Finally, it is worth noting that $\mathcal{E}_{\tau_T}$ is a result of the approximately sparse structure of $B$; in the special case where $B$ is exactly sparse, this term would be 0.
Corollary (ref) gives the bound of $\Delta_B$ and $\Delta_\Theta$ with specific choice of the thresholded level $\eta$ , when the true value $B^\star$ lies in the $\ell_q$ ball of radius $R_q$ (see definition in Equation (ref)):
Next, we provide high probability bounds/concentrations for key quantities associated with the derived error bound in Section (ref), for random Gaussian realizations of the underlying factor and error processes. Specifically, this involves the verification of the RSC condition, as well as the examination of quantities associated with the deviation condition to which the choice of $(\lambda_B,\lambda_\Theta)$ needs to conform, as shown in (ref).
We introduce additional notation for the subsequent technical developments. For some generic process $\{X_t\}$, in addition to the auto-covariance function $\Gamma_X(h)$ and its spectral density $g_X(\omega)$, we define its maximum and minimum eigenvalue associated with the spectral density $g_X(\omega)$ introduced in Section (ref) as follows basu2015estimation:
For two generic centered processes $\{X_t\}$ and $\{Y_t\}$ that are assumed jointly covariance stationary, whose spectral density is given by $g_{X,Y}(\omega):=\tfrac{1}{2\pi}\sum_{h=-\infty}^{\infty} \Gamma_{X,Y}(h)e^{i\omega h}$ where $\Gamma_{X,Y}(h)=\mathbb{E}(X_tY_{t+h}^\top)$, the upper extreme for $g_{X,Y}(\omega)$ is analogously defined as
In general $g_{X,Y}(\omega)\neq g_{Y,X}(\omega)$, but $\mathcal{M}(g_{X,Y})=\mathcal{M}(g_{Y,X})$.
For the processes involved in our proposed model, recall that $\{X_t\}$, $\{\epsilon_t\}$ and $\{F_t\}$ are mean zero Gaussian processes. In particular, $\{\epsilon_t\}$ is a noise process that does not exhibit temporal nor cross-sectional dependence, hence it is effectively a Gaussian random vector with covariance $\Sigma_\epsilon=\sigma_\epsilon^2 \mathrm{I}_p$, and its spectral density simplifies to $g_\epsilon(\omega)=\frac{\Sigma_\epsilon}{2\pi}$. Further, we define the shifted process $\{\widetilde{\epsilon}_t:=\epsilon_{t+1}\}$ for notation convenience.
The following lemma verifies that with high probability, for random realizations of the process $\{X_t\}$, the RSC condition is satisfied provided that the sample size is sufficiently large:
The next lemma establishes a high probability bound for the interaction term $\mathbf{X}_{T-1}^\top \mathbf{E}/T$ that influences the choice of $\lambda_B$ through its elementwise $\ell_\infty$ norm.
Note that with the definition of the shifted processes $\{\widetilde{\epsilon}_t\}$, we have $g_{X,\widetilde{\epsilon}}(\omega) = e^{-ih\omega} g_{X,\epsilon}(\omega)$, which implies $\mathcal{M}(g_{X,\widetilde{\epsilon}})=\mathcal{M}(g_{X,\epsilon})$. Hence, the term that measures the upper extreme of the cross-spectrum between $X_t$ and the shifted process in (ref) can be replaced by its unshifted counterpart. Moreover, since $g_\epsilon(\omega)=\tfrac{\sigma_\epsilon}{2\pi}$, its upper extreme is given by $\mathcal{M}(g_\epsilon) = \Lambda_{\max}(\Sigma_\epsilon)/(2\pi)$.
The next lemma provides an upper bound for the maximum eigenvalue of the sample covariance matrix.
Proofs for Lemmas (ref) to (ref) can be found in Appendix (ref).
Up to this stage, we have verified the RSC condition and obtained the high probability bounds for quantities that are associated with the choice of $(\lambda_B,\lambda_\Theta)$, for random realizations from the underlying processes. Theorem (ref) combines the results in Corollary (ref) and Lemmas (ref) to (ref), and provides a high probability error bound of the estimates when the data are random realizations from the underlying processes, as stated next.
The convergence property of Algorithm (ref) can be established using familiar arguments and exploiting its convex nature. Specifically, define the objective function is given by
and is {\em jointly convex} in $(B,\Theta)$, with a convex feasible region $\mathbb{B}_\infty(\phi,\mathbf{X}_{T-1})$. Thus, it directly follows from tseng2001convergence that the alternating minimization that generates the sequence $\{(\widehat{B}^{(k)},\widehat{\Theta}^{(k)})\}$ converges to a stationary point which is also a global optimum, though the global optimum is not necessarily unique.
To conclude this section, we remark that the theoretical formulation in (ref) can be solved in an analogous way to Algorithm (ref). Specifically, the update of $\Theta$ requires modification to satisfy the constraint on the feasible region of $\Theta$, and the partial minimization can be solved by employing the composite gradient descent algorithm of nesterov2007gradient that involves singular value thresholding steps. Nevertheless, the modified algorithm is also convergent, as the one in Algorithm (ref).
In this section, we present results for simulation studies under various settings to demonstrate the performance of our proposed model. As comparison, we also present the common space (to be defined later) recovery error and the one-step-ahead forecast error across our proposed method, the vanilla factor analysis using PC method a la stock2005implications, and the proposed method in forni2005generalized.
\paragraph{An empirical algorithmic relaxation.} The actual implementation of Algorithm (ref) requires $\lambda_B$, $\lambda_\Theta$ as inputs, which in practice are challenging to select. On the other hand, the computation procedure designed for solving the convex program in (ref) suggests that to obtain the estimates boils down to alternating between the following two steps: (1) a regularized regression (lasso) update on the rows of $B$; and (2) an SVT update on $\Theta$. This naturally motivates the following steps in the implemented version of the algorithm, outlined next in Algorithm (ref).
Algorithm (ref) outlines the algorithmic relaxation to obtaining $(\widehat{B},\widehat{\Theta})$ in (ref), and it can be viewed as an alternating minimization algorithm that solves
For each update, the partial minimization step with respect to $\Theta$ or $B$ ensures that the value of the objective function is always non-ascending, which together with the fact that the objective function is bounded below guarantees convergence of the objective function iterates. In practice, the algorithm is terminated when the descent magnitude of the objective function between successive iterations is smaller than some pre-specified tolerance level. This algorithm does not provide guarantees of convergence to a stationary point of the sequence of $(\bar{\Theta}^{(k)},\bar{B}^{(k)})$ iterates, which requires stronger assumptions --- either the convexity of the objective function and the constraint region, or the uniform compactness of the generated sequence of iterates.
\paragraph{Choice of the tuning parameter $\lambda_B$ and the rank constraint $r$.} The implementation of Algorithm (ref) requires a specific pair of $(\lambda_B,r)$ as input. We consider choosing the optimal pair of $(\lambda_B,r)$ based on the information criterion proposed in ando2015selecting, called the Panel Information Criterion (PIC) and defined as:
where $\widehat{\sigma}^2 = \tfrac{1}{np}{\vert\kern-0.25ex\vert\kern-0.25ex\vert \mathbf{X}_n-\widehat{\Theta}_{\text{emp}} - \mathbf{X}_{n-1}\widehat{B}_{\text{emp}}^\top \vert\kern-0.25ex\vert\kern-0.25ex\vert}_{\text{F}}^2$ and $(\widehat{B}_{\text{emp}},\widehat{\Theta}_{\text{emp}})$ are solutions to (ref) with the specific pair of plug-in $(\lambda_B,r)$. The optimal pair $(\lambda_B,r)$ is then selected in two steps: in step 1, we obtain $(\lambda_{B}^0,r^0)$ that gives the smallest PIC over a lattice $\mathcal{G}_{\lambda_B}\times \mathcal{G}_r := \{\lambda_B^{(1)},\dots,\lambda_B^{(j_1)}\}\times \{r^{(1)},\dots,r^{(j_2)} \}$; in step 2, we fix $r$ at $(d+1)\times r^0$ where $d$ is the number of lags corresponding to the sparse $\mathrm{VAR}(d)$ model, and seek for $\lambda_B^{\text{opt}}$ over a grid that minimizes PIC$(\lambda_B, (d+1)r^0)$. The optimal pair of tuning parameters is then given by $(\lambda_B^{\text{opt}},r^{\text{opt}}):=(\lambda_B^{\text{opt}},(d+1)\times r^0)$.
\paragraph{Data generating mechanism.} Synthetic data are generated according to the lag-adjusted factor model representation $X_t = \Lambda F_t + BX_{t-1} + \epsilon_t$. Starting from the standard approximate factor model representation $X_t = \widetilde{\Lambda} f_t + u_t$, $u_t$ is serially correlated and follows a $\text{sparse }\mathrm{VAR}(d)$ model\footnote{Throughout this section, we assume $d=1$; additional results for $d>1$ have been deferred to Supplement (ref).}, at each timestamp $t$, the $K$-dimensional factor is generated according to a $\mathrm{VAR}(q)$ model $f_t = \Phi_1 f_{t-1} +\cdots+\Phi_q f_{t-q} + \eta_t$ where $\eta_t\sim \mathcal{N}(0,\sigma_\eta^2\mathrm{I})$; decorrelating $u_t$ leads to following dynamic of $X_t$:
where $\Lambda = [\widetilde{\Lambda},B\widetilde{\Lambda}]$ and $F_t = (f_t^\top,f^\top_{t-1})^\top\in\mathbb{R}^{2K}$.
We consider several simulation settings as listed in Table (ref) to test various facets of the model, primarily encompassing the dimensionality of the system $p$ and the number of factors $K$, as well as the sparsity structure of $B$ and its spectral radius that captures the level of autocorrelation. In addition to settings S0 to S4 where $\epsilon_t$ is Gaussian, to test the robustness of the proposed model to the presence of heavy tails, we consider also cases where $\epsilon_t$ follows some multivariate $t$ distribution (S5 to S7). Throughout all numerical experiments presented in this section, the sample size $T$ is fixed at 200 and the spectral radius of the $\mathrm{VAR}(q)$ system is set randomly from $\mathsf{Unif}[0.6,0.8]$.
To generate the sparse transition matrix $B$, for each row that corresponds to the coefficients of each single time series regression, its (strong) support set is randomly generated to meet the specified density level (i.e., $2/p$ or $2/5$), and nonzero entries are then generated from $\pm\mathsf{Unif}[m_B-0.1,m_B+0.1]$. In the case of a weakly sparse $B$, entries in the weak support set are generated from $\mathsf{Unif}([-10\%m_B,10\%m_B])$. Finally, all entries are scaled to meet the specified $\varrho(B)$ level, to ensure that the system is stationary. For the dense factor loading matrix $\Lambda$, its entries are generated from $\pm\mathsf{Unif}[m_\Lambda-0.1,m_\Lambda+0.1]$. It is worth noting that the value of $m_B$ and $m_\Lambda$ are set so that the factor/lag space relative strength is satisfied, measured by the empirical relative signal-to-noise ratio for the $\Lambda F_t$ and the $BX_{t-1}$ component.
\paragraph{Performance evaluation.} To measure the accuracy of the obtained estimates and forecast, we focus on the following four components of the model:
As Table (ref) demonstrates, for all three components, estimates obtained from Algorithm (ref) exhibit good performance. In particular, (i) the proposed method is robust to the sparsity structure of $B$, as both exactly-sparse and weakly-sparse settings yield very satisfactory strong support recovery (see S1, S2 and S4). (ii) A larger panel size $p$ leads to improved factor hyperplane recovery, as manifested in the form of smaller relative error in magnitude estimation although it requires the sparsity of the transition matrix to decrease accordingly (recall that it is set to $2/p$); however, the performance deteriorates as the dynamics of $f_t$ become more complex (e.g., S4). (iii) A strong signal in the lag-space leads to improved recovery of $B$, despite the presence of stronger temporal dependence which empirically incurs the algorithm to take more iterations to converge (e.g., S3). For all settings, PIC correctly selects the number of factors, which translates into the correct identification of the rank constraint.
Next, we compare the performance of common space recovery and forecasting for the following three methods: the posited model, standard SW formulation and GDFM. For SW stock2005implications, the reported error is based on the minimum error among estimates obtained under different rank constraints ranging between $K$ and $2K$; for GDFM forni2005generalized, the reported error is based on the minimum error among estimates obtained under different combinations of $(\widetilde{q},\widetilde{r})$ that determines the number of common factors when loaded dynamically and contemporaneously. For all settings, the proposed method (lag-adjusted DFM) outperforms the other two by explicitly incorporating the lag space spanned by $X_{t-1}$; specifically, it outperforms its competitors by a wide margin when the lag space possesses a stronger signal, in which case SW becomes particularly susceptible (e.g., S3 and S6). However, as the dynamics of $f_t$ becomes more involved, its advantage becomes less pronounced (S4).
Finally, the proposed model is relative robust to the presence of heavy tails, although the performance deteriorates compared to the Gaussian case. Specifically, when the distribution shows significant deviation from Gaussian (e.g., S5), the degradation manifests itself through less satisfactory recovery in the support of $B$ and larger error of the estimated factor space; whereas the forecasting performance isn't affected. On the other hand, with lighter tails (e.g., S6), the performance becomes comparable to the Gaussian case.
To further compare the performance across all three methods, we consider settings where data generating processes deviate from the proposed model in (ref). Specifically, we adopt the data generating mechanism in forni2017dynamic, that is,
where coordinates of $u_t$ and $\xi_t$ are i.i.d standard Gaussian white noises processes, with $u_t$ capturing the structural shocks. Entries of $\Lambda$ and $K$ are drawn independently from $\mathsf{Unif}[-1,1]$; entries of $D$ are first drawn independently from $\mathsf{Unif}[-1,1]$ then scaled so that the spectral norm of $D$ satisfies some pre-specified target, with the latter drawn from $\mathsf{Unif}[0.4,0.9]$. We focus on the performance of common space recovery and the one-step-ahead forecast, under various combinations of the model parameters, as listed in Table (ref):
As the results show, with the data generating procedure deviating from the proposed model in (ref), with properly chosen tuning parameters, the performance of the proposed methodology matches that of SW by effectively having $\widehat{B}=0$. Meanwhile, it is worth noting for all three methods, the performance of common space recovery shows significantly less variability compared with that of forecasting; in particular, the forecasting performance of GDFM exhibits the highest variance across replications, among the three methods.
Factor models have been widely used in financial applications. In particular, they have been employed in analyzing the dynamics of asset returns, either for the purpose of identifying risk factors, or for estimating the covariance structure amongst assets for better portfolio diversification and asset allocation fan2012vast. We applied the proposed modeling framework to a set of stocks return data corresponding to 75 large US financial institutions, which also exhibit strong (serial) correlation in the error terms. Specifically, we analyze the risk-free returns\footnote{The risk-free return of Stock $i$ at time $t$ is calculated as $\widetilde{r}_{i,t} = r_{i,t} - r_{\text{rf},t} = \frac{p_{i,t}-p_{i,(t-1)}}{p_{i,(t-1)}}- r_{\text{rf},t}$, where $p_{i,t}$ is its stock price at time $t$ and $r_{\text{rf},t}$ is the risk-free rate.} of 25 banks, 25 insurance companies and 25 broker/dealer firms for the period of 2001-17. Note that this time period contains a number of significant events for the financial industry, including the growth of mortgage bank securities ashcraft2010MBS in the early 2000s, rapid changes in monetary policy in 2005-06, the great financial crisis eichengreen2010tale in 2008-09 and the European debt crisis in 2011-12 and their aftermath. Our analysis identifies a number of interesting patterns, especially around the period 2007-09 encompassing the beginning, height and immediate aftermath of the US financial crisis, both through changes in the factor structure and the partial autocorrelation one governed by the VAR model transition matrix of the log-returns of these financial assets.
\paragraph{Data.} The data consist of weekly stock return data corresponding to 75 large financial institutions in terms of market capitalization, for the period of January 2001 to December 2017 and were obtained from the Center for Research in Security Prices (CRSP) database. The 75 companies are categorized into three sectors: banks (SIC code 6000--6199), broker/dealers (SIC code 6200--6299) and insurance companies (SIC code 6300--6499), with 25 in each sector billio2012econometric. As we require that the data be available for the entire time span under consideration, 56 firms are kept for further analysis, since the remaining ones either went bankrupt or were forced to merge with financially healthier companies (e.g. Lehman Brothers and Merill Lynch in 2008, resp.). To get an overview of the correlation structure amongst the stocks after accounting for the first principal component that captures the weighted average return of the portfolio they constitute avellaneda2010statistical, we plot the correlation among the principal component regression residuals. Specifically, the entire time span is broken into three sub-periods that have been previously considered in the literature billio2012econometric: 2001--2006 (pre-crisis), 2007--2009 (crisis), 2010-2017 (post-crisis), and plot the correlation maps corresponding to samples in each period. As Figure (ref) demonstrates, overall, we observe positive correlation within each sector and negative correlation across them. Such a structural pattern is predominant in the pre-crisis period especially within the insurance sector, and becomes significantly weaker in the post crisis one; whereas during the crisis, stronger negative correlation across blocks is present as well as scattered positive correlations. This suggests that different factor and auto-regressive structures emerge during the crisis period. Further, note that similar results hold if we examine the residuals after removing a second principal component, so as to capture a larger percentage of variance of the stock returns.
The analysis is based on 104-week-long rolling windows to avoid issues with non-stationarity that potentially depends on length of the period under consideration. This strategy has also been used in billio2012econometric, lin2017regularized and allows monitoring change in the number of factors over time, as well as the sparsity level of $B$ which measures the connectivity of the partial autocorrelation network across these financial institutions. Note that 48% of the rolling samples fail to reject the null hypothesis that they are multivariate normality\footnote{we consider the Henze-Zirkler's multivariate normality test.}, with the exception of samples during the crisis period and those at the end of the sampling period post-crisis (around 2015--17). We fit the proposed lag-adjusted factor model in each time window, with tuning parameters selected according to a modified PIC criterion\footnote{The criterion is modified to $\text{PIC}^*(B,r) = \log\widehat{\sigma}^2 + \big(\frac{\log T}{T}\big)\|B\|_0 + r\cdot \big(\frac{T+p}{Tp}\big)\log(Tp)$.} that does not depend on the range within which the number of factors is being searched.
As Figure (ref) shows, sharp changes are observed in the temporal dependence structure of stock returns during the crisis period. In particular, two change points respectively correspond to the beginning of the 2007 sub-prime mortgage crisis and the ending of the 2008--2009 global financial crisis. Specifically, for the pre- and post-crisis periods, the density of the transition matrix stays at a level close to zero, suggesting that not much serial correlation exists in the idiosyncratic component after the common factor (proxy for the market portfolio) is accounted for. During the crisis period, however, the connectivity level of $\widehat{B}$ witnesses a sharp increase, reaching its maximum in the sampling window corresponding to Dec 2006--Dec 2008, during which period multiple major events of the financial crisis occurred. We also track the change in R-squared and the R-squared attributed to the factor over time, as a surrogate for the quality of the model fit, as shown in Table (ref).
In accordance with the connectivity level of the VAR component that captures the temporal dependence among the stocks, under normal market conditions, the majority fit of the model (R-squared) comes from the factor hyperplane; whereas during the crisis period, the gap between the total R-squared and the factor R-squared widens with the lag term explaining a non-negligible proportion of the R-squared, which indicates the presence of significant cross autocorrelations in the returns.
To further investigate the factor composition and the temporal dependence structure during the crisis period, we zoom in on the sampling frequency and focus on the year of 2008. Specifically, we consider daily data from January 2008 to December 2008 that cover 253 consecutive trading days and fit the proposed lag-adjusted factor model. Note that for this part of the analysis, the sample consists of 72 stocks. Using $\text{PIC}^*$, 2 factors are identified, with the dominant one capturing 55% of the R-squared followed by 11% from the second. For reconstruction purposes, we assume they are orthogonal so that the factor composition can be retrieved from the singular vectors of $\widehat{\Theta}$.
As depicted in the left panel of Figure (ref), all financial institutions contribute positively to the first factor, with dominating contributors spread in all sectors. The composition of the second factor shows an interesting pattern: two negative contributors are FRE (Freddie Mac) and FNM (Fannie Mae), and the positive ones are primarily in the insurance sector. However, AIG---unlike its peers---shows almost zero contribution to the second factor, albeit its strong contribution to the first one. The latter is consistent with other findings that it played a prominent role during the crisis.\footnote{According to an estimate as of January 2010, AIG accounted for 38% of the total losses incurred by insurance companies (98.2 out of 261.0 billions) since 2007. Source: Bloomberg, see also schich2010insurance.} In the right panel of Figure (ref), we plot the partial auto-correlation network of the firms during the crisis after properly thresholding the entries that have small magnitudes, with red edges denoting positive links and with blue negative ones. Nodes that belong to the same sector are colored identically. A careful examination of the node weighted in/out-degrees shows that the top emitters are relatively uniform, in the sense that their weighted out-degrees do not differ by much; whereas the top receivers are dominant, since the weighted in-degrees for top receivers are significantly higher compared with the rest. Further top emitters heavily concentrate in the insurance sector. Meanwhile, some of the top receivers are also major contributors to the factors' composition, e.g., AIG to the 1st factor, HIG to the 2nd, etc. This finding partially aligns with the role that many insurance companies played in magnifying the impact of the crisis on the overall stability of the financial system, due to their large insurance underwriting of Credit Default Swaps and subsequent exposure to accentuated risks eichengreen2010tale. However, this analysis points to the importance of insurance companies based on publicly available data and before their role in the crisis was fully revealed and understood. It is worth noting that with the same set of data, vanilla factor analysis using the information criterion proposed in bai2002determining only identifies 1 factor, which further substantiates the aforementioned point that classical factor analysis may lead to skewed inference when strong correlation among the idiosyncratic component is present.
To conclude this section, we compare and contrast our results with those obtained in billio2012econometric, in which the authors consider 100 financial institutions comprising of the largest 25 among each of the four categories: hedge funds, broker/dealers, banks and insurers; thus, that data set is enhanced by the inclusion of big hedge funds for which publicly available stock quotes are not accessible. From the systemic risk standpoint, the authors measure the connectedness of the system based on principal component analysis (PCA) and Granger-causality network analysis during the 1994--2008 period, and identify increased level of interconnectedness during the crisis period and the asymmetry in the degree of connectedness amongst different sectors. Our results are qualitatively similar to these results, and the conclusions broadly match. However, we would like to highlight some key differences in both modeling and in the empirical results obtained. From the modeling perspective, billio2012econometric consider two separate modeling strategies: (i) a Principal Components Analysis (akin to a static factor model) and (ii) a Granger-causality based analysis through fitting a VAR model for each pair of stocks returns. The PCA analysis examines a {\em fixed} number of principal components/factors and the authors argue that the increasing proportion of variation explained by them is an indication of the systematic response of the financial system to the crisis. Their pairwise based Granger-causal network also reveals increased connectivity during the crisis period. Our model considers latent factors and lead-lag relationships among stock returns {\em simultaneously}, thus gaining better and more informative insights. In addition, the lead-lag relationships are considered across all firms simultaneously rather than in a pairwise fashion. By incorporating the strong correlations present in the idiosyncratic component, our model is more parsimonious. Specifically, during the crisis, billio2012econometric uses 10 principal components to account for 85% of the returns variance, whereas only 5 suffice in our model; further, the leading factor in their analysis only accounts 37% of the variance, compared to 50% in our model. Finally, extending the analysis period to 2017 shows that after 2011 the influence of banks and insurance companies on stock returns waned, as the marker slowly returned to normalcy. However, we are in broad agreement with the billio2012econometric conclusion on the heightened role of banks and insurers up to 2009.
In this paper, we introduced a novel modeling framework that generalizes the classical approximate factor model to include lags of the observable process, so that stronger correlations among the idiosyncratic component can be accommodated. The autoregressive structure is assumed to be sparse, which enables its estimation for large time series panels. Estimation of the model parameters is based on a maximum likelihood formulation that leads to a convex optimization problem, and the resulting estimates come with high probability error bounds guarantees that can be expressed in terms of key structural parameters ($n, p, K, s$, etc.), and exhibit superior empirical performance in synthetic data.
In addition to generalizing the model in stock2005implications, our proposed model can also be perceived as a robust treatment of endogeneity. Specifically, as noted by anderson2007forecasting, in the presence of large values in $u_t$ and for a relatively small panel size $p$, the factor estimates will be distorted as a result of this endogeneity. In this work, by explicitly taking into consideration the lagged terms in the dynamics of $X_t$, the noise term $\epsilon_t$ becomes strictly exogenous. Our proposed model and estimation procedure has the capacity of handling much stronger correlation between $F_t$ and $u_t$, although ultimately we do require $\text{Cov}(F_t,u_t)$ to be indirectly bounded in some appropriate way.
\setcounter{page}{1} \pagenumbering{arabic}
\titleformat{\section}{\normalfont}{APPENDIX \arabic{section}.}{0.5em}{\uppercase} \numberwithin{lemma}{section} \numberwithin{remark}{section}