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.
85,181 characters · 15 sections · 29 citation commands
Optimal Estimation of Large-Dimensional Nonlinear Factor Models
Keywords: nonlinear factor model, latent variables, low-rank method, high-dimensional data, principal component analysis
\thispagestyle{empty}
\onehalfspacing \setcounter{page}{1}
\pagestyle{plain}
High-dimensional data have become increasingly available due to technological advances in data collection, which are typically characterized by a large number of cross-sectional units and a large number of features. Factor analysis is a useful tool for summarizing information in such big data sets and has wide applications in statistics, economics, and many other data science disciplines.
One crucial insight of factor analysis is that much of the data variation can be explained by the interaction between a few important, usually unobserved characteristics associated with two dimensions (“cross section” and “features”). For example, a classical factor model has the following representation \[ x_{il}=\bm{f}_l'\bm\alpha_i+u_{il}, \quad 1\leq i\leq n, \; 1\leq l\leq p, \] where $x_{il}$ is the $l$th feature of the $i$th cross-sectional unit; $\bm{f}_l\in\mathbb{R}^r$ is usually termed common factors, which only vary across $l$; $\bm\alpha_i\in\mathbb{R}^r$ is termed factor loadings, which only vary across $i$; and $u_{il}$ is some idiosyncratic error. Both the number of cross-sectional units $n$ and the number of features $p$ are large. Such specifications are common in many problems. For instance, in panel data analysis $x_{il}$ may be a variable collected over time, and $\boldsymbol{\lambda}_l$ and $\boldsymbol{\alpha}_i$ are time-specific and individual-specific effects respectively; in analysis of a recommender system $x_{il}$ may be a variable representing the preference of an individual $i$ for an item $l$, which is explained by item-specific features $\boldsymbol{\lambda}_l$ and individual-specific features $\boldsymbol{\alpha}_i$; and in causal inference and program evaluation $x_{il}$ could be repeated measurements of some underlying unobserved confounders $\boldsymbol{\alpha}_i$, and $\boldsymbol{\lambda}_l$ is a (linear) transformation specific to the $l$th measurement.
The linear structure above, however, is usually restrictive and may be unrealistic in many problems. For example, test scores in multiple subjects or time periods are often used to measure fundamental abilities of students in empirical research Cunha-Heckman-Schennach_2010_ECMA. It is difficult to justify a linear relationship between test scores and unobserved abilities in this context. A more appealing approach is to consider a possibly nonlinear factor model Yalcin-Amemiya_2001_SS
The observed variables $x_{il}$ (e.g., test scores) and the unobserved variables $\boldsymbol{\alpha}_i$ (e.g., latent abilities) are linked through a “production function” $\eta_l$, whose functional form is left unspecified. This setup encompasses the classical linear factor model Bai-Wang_2016_ARE as a special case, but also allows for other possibly nonlinear relationships between the observables and the unobservables.
In this paper we analyze the nonlinear factor model (ref) based on local principal subspace approximation. The procedure begins with $K$-nearest neighbors ($K$-NN) matching approximation for each unit on the observed features $x_{il}$'s. Within each local neighborhood formed by the $K$ matches, the underlying possibly nonlinear factor structure is then approximated by a linear factor structure which can be estimated using principal component analysis (PCA). Given the locality of this procedure, we term it local principal component analysis in this paper.
This methodology has an intuitive geometric interpretation. The set of latent functions $\{\eta_l:1\leq l\leq p\}$ generates a low-dimensional, possibly nonlinear “surface” embedded in a high-dimensional space, if the number of features $p$ is large but the number of latent variables $r$ is small. Suppose that different values of latent variables $\boldsymbol{\alpha}_i$ can induce non-negligible differences in many observed features $x_{il}$'s. In this case, the $K$ nearest neighbors of each unit as appropriately measured by the observed features should also be close in terms of the latent variables, thus forming a local neighborhood in the high-dimensional space. Furthermore, such units are approximately lying on a local principal surface, up to errors governed by the number of nearest neighbors and the number of principal components extracted in each local neighborhood. The availability of many observables as “measurements” of the latent variables $\boldsymbol{\alpha}_i$ is crucial for validity of this approximation, which affects the matching discrepancy of nearest neighbors and the estimation error of local principal components.
The idea of locally approximating a nonlinear latent surface embedded in a high-dimensional space has been widely used in the modern machine learning literature and is popular in applications such as face recognition, motion segmentation, and text classification. Typical examples include Zhang-Zha_2004_SIAM, Peng-Lu-Wang_2015_NN,Arias-Lerman-Zhang_2017_JMLR, among others. This paper provides a formal theoretical foundation for such methods that use the similar idea of local approximation.
We establish the statistical properties of our proposed local PCA in large-dimensional settings, which makes several contributions to the literature. First, we derive a sharp bound on the implicit discrepancy of nearest neighbors in terms of the latent variables (Theorem (ref)). This result is established under a generic choice of the distance function, encompassing and extending previous studies of matching techniques based on specific metrics Zhang-Zha_2004_SIAM, Zhang-Levina-Zhu_2017_BIMA. Crucially, we show that the closeness of nearest neighbors, indirectly obtained through matching on noisy measurements, relies on two conditions: (i) the selected distance can “denoise” the data to reveal the latent factor structure $\eta_l(\boldsymbol{\alpha}_i)$, and (ii) the distance of the noise-free structure $\eta_l(\boldsymbol{\alpha}_i)$ needs to be informative about that of the unobserved variables $\boldsymbol{\alpha}_i$. Therefore, our first contribution provides theoretical guidance for practitioners who have to rely on noisy measurements to match on some unobserved variables of interest.
Second, we derive the sup-norm error bounds for the estimated local factors and factor loadings obtained by applying PCA to nearest neighbors (Theorem (ref)). The target quantities characterize the latent functions $\eta_l$'s and the latent variables $\boldsymbol{\alpha}_i$'s respectively. Importantly, due to the potential nonlinearity, the strength of factors in the local approximation is possibly heterogeneous and needs to be properly taken into account. This result appears to be new in the literature, complementing the studies of linear factor models with weak or semi-strong factors onatski2012asymptotics, Wang-Fan_2017_AoS, Abbe-et-al_2020_AoS. Furthermore, we note that though the latent functions $\eta_l$'s and latent variables $\boldsymbol{\alpha}_i$'s cannot be separately recovered without additional restrictions, the factors and loadings from local PCA suffice for flexible out-of-sample forecasts based on nonparametric regression and can be used in general causal inference problems with mismeasured confounders miao2018identifying,nagasawa2018treatment. The details of this method are discussed in Feng_2023_wp.
Third, building on the first two results, we show that local PCA can consistently estimate the nonlinear factor structure $\eta_l(\boldsymbol{\alpha}_i)$, deriving a convergence rate that is uniform over both individuals and features (Corollary (ref)). Under rather general conditions, the local PCA estimator can attain a uniform convergence rate that coincides with the optimal one for the infeasible cross-sectional nonparametric estimation of the heterogeneous functions $\eta_l$'s stone1982optimal. To the best of our knowledge, this paper is the first to show this rate can be achieved in this general nonlinear factor model, contributing to the literature on low-rank approximation of data matrices udell2019big, fernandez2021low.
Fourth, in Section (ref) and the online Supplemental Appendix we extend the basic nonlinear factor model (ref) by including observable regressors that have high-rank variation in both dimensions. This provides a new tool for studying, for example, linear regression models with nonlinear fixed effects, complementing the vast literature on panel regression with interactive fixed effects Bai_2009_ECMA,Bai-Li_2014_AoS.
Finally, we apply local PCA to certain matrix completion problems with a few missing entries (Theorem (ref)). An empirical application studying the effect of tax cuts on economic growth is used to illustrate the potential usefulness of the proposed method in policy evaluation settings such as synthetic controls Abadie_2020_JEL
The rest of the paper is organized as follows. In Section (ref) we formally set up the nonlinear factor model. Section (ref) gives a detailed description of the estimation procedure. Section (ref) presents the main theoretical results. Section (ref) summarizes Monte Carlo results. An empirical application to synthetic control problems is given in Section (ref). Section (ref) concludes. The appendices collects several technical results that may be of independent interest, including properties of several distance functions, verification of the local approximation of latent factor structure (Appendix (ref)), and selected proofs for the main results (Appendix (ref)). The online Supplemental Appendix (SA hereafter) contains all omitted proofs and additional technical and numerical results. Replications of the simulation study and empirical illustration are available at \url{https://github.com/yingjieum/Replication_NonlinearFactorModel_2023}.
Let $\bm{x}_i=(x_{i1}, \cdots, x_{ip})'$ be a $p$-vector of observed variables for the $i$th unit in the sample. We can write the nonlinear factor model as
where $\boldsymbol{\alpha}_i\in\mathbb{R}^{r}$ is a vector of latent variables, $\boldsymbol{\eta}=(\eta_{1}, \cdots, \eta_{p})': \mathbb{R}^r\mapsto\mathbb{R}^p$ is a vector of latent functions, and $\bm{u}_{i}=(u_{i1},\cdots,u_{ip})'$ is the idiosyncratic error. In matrix notation, $$\bm{X}=\bm{H}+\bm{U}$$ where $\bm{X}=(\bm{x}_1,\cdots, \bm{x}_n)$, $\bm{H}=(\boldsymbol{\eta}(\boldsymbol{\alpha}_1),\cdots, \boldsymbol{\eta}(\boldsymbol{\alpha}_n))$ and $\bm{U}=(\bm{u}_1,\cdots, \bm{u}_n)$ are $p\times n$ matrices. Throughout the paper, $(\boldsymbol{\alpha}_i: 1\leq i\leq n)$ and $\boldsymbol{\eta}(\cdot)$ are understood as random elements, which generate the $\sigma$-field $\mathscr{F}$. Our main analysis below is conducted conditional on $\mathscr{F}$. In this sense, $\boldsymbol{\alpha}_i$ and $\boldsymbol{\eta}(\cdot)$ are akin to the “fixed effects” commonly incorporated in panel data models.
The usual linear factor model, also known as the interactive fixed-effect model, is covered as a special case by this setup, where the latent function is assumed to be linear in $\boldsymbol{\alpha}_i$, e.g., $\eta_l(\boldsymbol{\alpha}_i)=\bm{f}_l'\boldsymbol{\alpha}_i$ for some $\bm{f}_l\in\mathbb{R}^r$, $l=1, \cdots, p$. By construction, the latent mean structure $\bm{H}$ is exactly low-rank ($r\ll p\wedge n$). In the more general case, however, the potential nonlinearity of the latent function $\boldsymbol{\eta}(\cdot)$ can make $\bm{H}$ full-rank, and traditional methods based on assuming a linear factor structure become inappropriate.
The main insight in nonlinear factor analysis is that due to the low-dimensionality of $\boldsymbol{\alpha}_i$, the variation of the large-dimensional $\bm{x}_{i}$ can still be explained by a few low-dimensional components in a possibly nonlinear way, which implies that $\bm{H}$ is still approximately low-rank. To gain some intuition, consider the local linear approximation approach widely used in the manifold learning literature (e.g., Zhang-Zha_2004_SIAM):
where $\nabla \boldsymbol{\eta}_l(\boldsymbol{\alpha}_i)$ is the vector of first-order partial derivatives of $\eta_l$, evaluated at $\boldsymbol{\alpha}_i$. This representation amounts to an approximately linear factor model with $r+1$ “factors” (i.e., an intercept plus $r$ partial derivatives of $\eta_l$), which motivates applying the usual principal component analysis to the local neighborhood of $\boldsymbol{\alpha}_i$.
More generally, if $\eta_l$ is sufficiently smooth, a higher-order approximation can be employed, which amounts to extracting more local factors from the approximation error in equation (ref). Intuitively, as derivatives of $\boldsymbol{\eta}(\cdot)$, the “factors” in such representations signify the degree of nonlinearity of the underlying latent structure, and the “loadings” reflect the magnitude of different approximation terms.
In practice, since $\boldsymbol{\alpha}_i$ is never observed by the researcher, one has to first employ some indirect strategy to construct the local neighborhood for each unit based on observables, and then conduct principal component analysis locally. This estimation procedure is described in Section (ref) and then theoretically formalized in Section (ref).
Now, before we proceed to the estimation procedure, we summarize the regularity conditions on the latent variables $\boldsymbol{\alpha}_i$, latent functions $\boldsymbol{\eta}(\cdot)$ and the error terms $\bm{u}_i$ in the next assumption, which are imposed throughout our main analysis.
Parts (a) and (b) are commonly used in the nonparametric regression literature. For simplicity, we assume all latent functions are sufficiently smooth, i.e., they belong to a H\"{o}lder class of order $\bar{m}$. Part (c) is a standard condition on idiosyncratic errors in factor analysis and graphon estimation. The independence requirement is imposed to simplify some analysis and can be further relaxed to allow for some weak correlation in one or both dimensions.
Derivatives. For a generic sequence of functions $h_j(\cdot)$, $j=1, \cdots, M$, defined on a compact support, let $\nabla^{\ell}\bm{h}_j(\cdot)$ be a vector of $\ell$th-order partial derivatives of $h_j(\cdot)$. The derivatives on the boundary are understood as limits with the arguments ranging within the support. When $\ell=1$, $\nabla\bm{h}_j(\cdot):=\nabla^1\bm{h}_j(\cdot)$ is the gradient vector, and the Jacobian matrix is $\nabla\bm{h}(\cdot):=(\nabla\bm{h}_1(\cdot), \cdots, \nabla\bm{h}_{M}(\cdot))'$.
Matrices. For a vector $\bm{v}\in\mathbb{R}^\mathsf{d}$, $\|\bm{v}\|=\sqrt{\bm{v}'\bm{v}}$ is the Euclidean norm of $\bm{v}$, and for an $m\times n$ matrix $\bm{A}$, $\|\bm{A}\|_{\max}=\max_{1\leq i\leq m, 1\leq j\leq n}|a_{ij}|$ is the entrywise sup-norm of $\bm{A}$. $s_{\max}(\bm{A})$ and $s_{\min}(\bm{A})$ denote the largest and smallest singular values of $\bm{A}$ respectively. Moreover, $\bm{A}_{i\cdot}$ and $\bm{A}_{\cdot j}$ denote the $i$th row and the $j$th column of $\bm{A}$ respectively. $\bm{1}_{\mathsf{d}}$ denotes the $\mathsf{d}$-vector of ones.
Asymptotics. For sequences of numbers or random variables, $a_n\lesssim b_n$ or $a_n=O(b_n)$ denotes $\limsup_n|a_n/b_n|$ is finite, $a_n\lesssim_\mathbb{P} b_n$ denotes $\limsup_{\varepsilon\rightarrow\infty}\limsup_n\mathbb{P}[|a_n/b_n|\\ \geq\varepsilon]=0$, $a_n=o(b_n)$ implies $a_n/b_n\rightarrow 0$, and $a_n=o_\mathbb{P}(b_n)$ implies that $a_n/b_n\rightarrow_\mathbb{P} 0$, where $\rightarrow_\mathbb{P}$ denotes convergence in probability. $a_n\asymp b_n$ implies that $a_n\lesssim b_n$ and $b_n\lesssim a_n$.
Others. For two numbers $a$ and $b$, $a\vee b=\max\{a,b\}$ and $a\wedge b=\min\{a,b\}$. For a finite set $\mathcal{S}$, $|\mathcal{S}|$ denotes its cardinality. For a $\mathsf{d}$-tuple $\bm{q}=(q_1, \cdots, q_{\mathsf{d}})\in\mathbb{Z}_{+}^{\mathsf{d}}$ and $\mathsf{d}$-vector $\bm{v}=(v_1, \cdots, v_{\mathsf{d}})'$, define $\bm{v}^{\bm{q}}=v_1^{q_1}v_2^{q_2}\cdots v_{\mathsf{d}}^{q_{\mathsf{d}}}$. We use $[m]$ to denote the set $\{1, 2, \cdots, m\}$ for any positive integer $m$.
This section describes the main procedure for local PCA which typically consists of two steps. First, choose a proper function $\rho:\mathbb{R}^p\times\mathbb{R}^p\mapsto\mathbb{R}$ to define the “distance” between different units in the sample based on the observed variables. We use the notion of distance in a loose sense, that is, $\rho(\cdot, \cdot)$ does not have to satisfy all axioms in the standard definition of the distance function. Second, apply principal component analysis to $K$ nearest neighbors of each unit defined by the distance calculation in the first step. See Algorithm \hyperlink{t2}{1} for a short summary. The main tuning parameters in this procedure are the number of nearest neighbors $K$ and the number of extracted principal components $\mathsf{d}_i$ in the local neighborhood of each unit $i$.
We recommend users separate the $K$ nearest neighbors matching and principal component analysis by row-wise sample splitting, which guarantees desired theoretical properties of local PCA as will be explained below. Specifically, split the row index set $[p]$ of $\bm{X}$ into two non-overlapping subsets: $[p]=\mathcal{R}^\dagger\cup\mathcal{R}^\ddagger$ with $p^\dagger=|\mathcal{R}^\dagger|$, $p^\ddagger=|\mathcal{R}^\ddagger|$ and $p^\dagger\asymp p^\ddagger\asymp p$. Accordingly, the data matrix $\bm{X}$ is divided into two submatrices $\bm{X}^{\dagger}$ and $\bm{X}^{\ddagger}$ with row indices in $\mathcal{R}^\dagger$ and $\mathcal{R}^\ddagger$ respectively. $\bm{H}^\dagger$, $\bm{H}^\ddagger$, $\bm{U}^\dagger$ and $\bm{U}^\ddagger$ are defined similarly. By Assumption (ref), $\bm{X}^\dagger$ and $\bm{X}^\ddagger$ are independent conditionally on $\mathscr{F}$. In principle, the two portions of the data only need to be “approximately” independent conditional on $\mathscr{F}$, thus allowing for weakly dependent errors. In practice, if rows of $\bm{X}$ are conditionally independent, one could, for example, randomly split the index set into two portions. For unit-time panel data, one could use, for example, the first half of time periods for $K$-NN matching and then apply PCA to the second half, thus respecting the original time series structure.
This step makes use of the subsample labeled by $\dagger$, i.e., the submatrix of $\bm{X}$ with row indices in $\mathcal{R}^{\dagger}$. For a generic unit $i\in[n]$, search for a set of indices $\mathcal{N}_i$ for its $K$ nearest neighbors (including $i$ itself) in terms of the distance $\rho(\cdot, \cdot)$:
Our theory is established for generic choices of the distance $\rho(\cdot,\cdot)$ under high-level conditions, which covers many usual choices in practice such as (1) the usual Euclidean distance $\rho(\bm{X}^{\dagger}_{\cdot i}, \bm{X}^{\dagger}_{\cdot j})=\frac{1}{\sqrt{p^\dagger}}\|\bm{X}^{\dagger}_{\cdot i}-\bm{X}^{\dagger}_{\cdot j}\|$, (2) the pseudo-max distance $\rho(\bm{X}^{\dagger}_{\cdot i}, \bm{X}^{\dagger}_{\cdot j})=\max_{l\neq i, j} |\frac{1}{p^\dagger}(\bm{X}^{\dagger}_{\cdot i}-\bm{X}^{\dagger}_{\cdot j})'\bm{X}^{\dagger}_{\cdot l}|$ proposed in Zhang-Levina-Zhu_2017_BIMA, and (3) the distance of “average” $\rho(\bm{X}_{\cdot i}^\dagger, \bm{X}_{\cdot j}^\dagger) =\frac{1}{p^\dagger}|\bm{1}_{p^\dagger}'(\bm{X}_{\cdot i}^\dagger-\bm{X}_{\cdot j}^\dagger)|$. These distance functions have different properties and may affect the behavior of the resulting nearest neighbors. For example, the pseudo-max distance can reveal the differences of units in terms of the noise-free factor structure even when the error $\bm{u}_i$ is (condition-on-$\mathscr{F}$) heteroskedastic, while the Euclidean distance cannot in this case. More detailed discussion is available in Section (ref) and Appendix (ref).
Note that when the observed variables differ in scale or importance for revealing information on the latent variables, it may be desirable to rescale or reweight different features when searching for nearest neighbors. Such transformations can be viewed as particular choices of the distance. See Remark (ref) below for more discussion.
In general, the “distance” function $\rho(\cdot,\cdot)$ needs to fulfill two purposes: (i) the distance of observables can be translated into that of noise-free factor structure (“denoising”), and (ii) the distance of the noise-free structure can be translated into that of the unobservables $\boldsymbol{\alpha}_i$. To grasp some intuition, take the pseudo-max distance as an example. This “metric” is defined based on averaging information across different features. If the error $u_{il}$ is independent or weakly dependent across $l$, their impact on the distance becomes negligible as the dimensionality $p^\dagger$ grows large, in which sense we “denoise” the measurements $\bm{x}_i$ and recover the noise-free component $\boldsymbol{\eta}(\boldsymbol{\alpha}_i)$. On the other hand, if $\boldsymbol{\eta}(\boldsymbol{\alpha}_i)$ is informative about $\boldsymbol{\alpha}_i$ in the sense of Assumption (ref) below, any two points found close in terms of the factor structure should also be close in terms of the underlying latent variables. Therefore, the nearest neighbors obtained by matching on the observables are similar in terms of the unobservables, which is the key building block of subsequent analysis.
Figure (ref) gives a conceptual illustration of this idea. An artificial two-dimensional surface is embedded in a three-dimensional space. If the requirements outlined before are satisfied, $K$-NN matching for a particular unit $i$ (colored in red) would generate a local neighborhood (the region within the red “ball”).
\FloatBarrier
\FloatBarrier
This step makes use of the subsample labeled by $\ddagger$, i.e., the submatrix of $\bm{X}$ with row indices in $\mathcal{R}^{\ddagger}$. Given a set of nearest neighbors $\mathcal{N}_i$ from the previous step, define a $p^\ddagger\times K$ matrix $\bm{X}_{\langle i \rangle}=(\bm{X}^{\ddagger}_{\cdot j_1(i)}, \cdots, \bm{X}^{\ddagger}_{\cdot j_K(i)})$. The subscript $\langle i \rangle$ indicates that the data matrix is defined locally for unit $i$.
If the neighbors we find for each unit $i$ are truly close in the latent variables, i.e., $\boldsymbol{\alpha}_j\approx\boldsymbol{\alpha}_i$ for all $j\in\mathcal{N}_i$, then the possibly full-rank matrix $\bm{H}_{\langle i \rangle}=(\bm{H}^{\ddagger}_{\cdot j_1(i)}, \cdots, \bm{H}^{\ddagger}_{\cdot j_K(i)})$ can be locally approximated by a low-rank structure: \[ \bm{H}_{\langle i \rangle}\approx \bm{F}_{\langle i \rangle}\boldsymbol{\Lambda}_{\langle i \rangle}', \] where $\bm{F}_{\langle i \rangle}\in\mathbb{R}^{p^\ddagger\times \mathsf{d}_i}$ and $\boldsymbol{\Lambda}_{\langle i \rangle}\in\mathbb{R}^{K\times \mathsf{d}_i}$ are termed local factors and local factor loadings respectively, and $\mathsf{d}_i$ is a user-specified parameter that governs the number of approximation terms. In general, we can write an approximately linear factor model:
where $\boldsymbol{\Xi}_{\langle i \rangle}$ is a matrix of corresponding approximation/smoothing bias, and $\bm{U}_{\langle i \rangle}=(\bm{U}^{\ddagger}_{\cdot j_1(i)}, \cdots, \bm{U}^{\ddagger}_{\cdot j_K(i)})$ is the idiosyncratic error matrix. This decomposition motivates the application of PCA to $\bm{X}_{\langle i \rangle}$:
such that $\frac{1}{p^\ddagger}\tilde\bm{F}_{\langle i \rangle}'\tilde\bm{F}_{\langle i \rangle}=\bm{I}_{\mathsf{d}_i}$ and $\frac{1}{K}\tilde\boldsymbol{\Lambda}_{\langle i \rangle}'\tilde\boldsymbol{\Lambda}_{\langle i \rangle}$ is diagonal.
The idea underlying (ref) is similar to the step of learning local tangent spaces in Zhang-Zha_2004_SIAM. The main difference is that $K$-NN matching and PCA in my procedure are conducted on different rows of $\bm{X}$. This is motivated by the fact that searching for nearest neighbors has implicitly used the information on $\bm{u}_i$. Without sample splitting, for units within the same local neighborhood, the nonlinear factor components $\bm{F}_{\langle i \rangle}\boldsymbol{\Lambda}_{\langle i \rangle}'+\boldsymbol{\Xi}_{\langle i \rangle}$ would be correlated with the noise $\bm{U}_{\langle i \rangle}$, rendering the standard PCA technique inapplicable. Row-wise sample splitting is a simple remedy, when the noise $u_{il}$ is independent (or weakly dependent) across $l$.
The idea of local PCA is illustrated in Figure (ref). Units around the red dot are approximately lying on a (local) linear tangent plane (colored in purple). Intuitively, this approximation is analogous to the local linear regression in the nonparametrics literature, though conditioning variables in this context are unobserved. More generally, if more leading local factors can be differentiated from the noise, then a local nonlinear principal surface can be constructed for a higher-order approximation of the underlying surface. However, the extracted principal components from the noisy data may not always be informative about the latent structure $\bm{H}_{\langle i \rangle}$. Instead, they could be partially or completely determined by the noise matrix $\bm{U}_{\langle i \rangle}$, and the potential local degeneracy of the nonlinear structure may further complicates this issue. Thus, we recommend users only take a few leading principal components associated with large eigenvalues. See formal discussion in Section (ref).
\FloatBarrier
\FloatBarrier
To theoretically formalize the estimation procedure in Section (ref), we need to (i) ensure the closeness between $\boldsymbol{\alpha}_j$ and $\boldsymbol{\alpha}_i$ for $j\in\mathcal{N}_i$, and (ii) show that the linear factor structure locally extracted from the observables are “consistent” in some proper sense for the nonlinear factor structure (“signals”) of our interest. The first task is nontrivial since $\boldsymbol{\alpha}_i$ is not observed by the researcher, and as described above, some indirect strategy is usually adopted to construct the desired neighborhood. Then, the key challenge is translating the distance measured in observables into that in unobservables under appropriate conditions. On the other hand, the second task is complicated by the fact that the factors in approximation (ref) are of different strength and might be degenerate at some point(s) $\boldsymbol{\alpha}\in\mathcal{A}$. In the following we will discuss each task and provide formal results.
The local neighborhood described in Section (ref) is constructed indirectly based on observed noisy measurements $\bm{x}_i$ of the latent $\boldsymbol{\alpha}_i$. To show the closeness of the resulting nearest neighbors, we typically need to guarantee the noise $\bm{u}_i$ is approximately negligible and the noise-free structure $\boldsymbol{\eta}(\boldsymbol{\alpha}_i)$ is informative about $\boldsymbol{\alpha}_i$, which are formalized in Assumption (ref).
Condition (a) formalizes the idea of “denoising” the data by choosing a proper distance $\rho(\cdot, \cdot)$. It guarantees that the distance of the observables between any pair of units is approximately determined by that of the noise-free components, up to a fixed constant $\rho_0$. Then, the closeness in terms of the observables can be translated into that of the latent factor structure. This requirement is usually mild, if we have many features and “average” them in a proper way.
On the other hand, condition (b) concerns the noise-free structure only. The upper bound condition is usually mild and can be deduced from the smoothness of $\boldsymbol{\eta}$ imposed in Assumption (ref), given a particular choice of $\rho(\cdot,\cdot)$. By contrast, the lower bound is the key requirement for $\bm{x}_i$ to be informative about $\boldsymbol{\alpha}_i$. Intuitively, it says if $\bm{H}^\dagger_{\cdot i}$ is close to $\bm{H}^\dagger_{\cdot j}$ in terms of the distance $\rho(\cdot,\cdot)$, $\boldsymbol{\alpha}_i$ needs to be close to $\boldsymbol{\alpha}_j$. The parameters $\underline{\varsigma}$ and $\bar{\varsigma}$ govern how the distance of the unobservables and that of the observables are linked. Typically, these requirements depend not only on the nonlinear factor structure itself, but also the chosen distance function. The idea underlying such informativeness requirement is also related to the completeness condition widely used in econometric identification problems schennach2020mismeasured. Roughly speaking, for a family of distributions, completeness requires that the density of a variable sufficiently vary across different values of the conditioning variable. Analogously, the lower bound in (b) amounts to saying that there is enough variation observed on the latent surface for different values of the latent variables.
To gain more intuition about the two conditions, we discuss several specific choices of $\rho(\cdot,\cdot)$ that are common in the literature. Formal technical results are deferred to Appendix (ref).
Now, we present our first main result, which characterizes the indirect matching discrepancy of nearest neighbors.
As shown in the above theorem, the matching can be made up to errors consisting of two terms in an asymptotic sense. The first part $(K/n)^{\bar{\varsigma}/(\underline{\varsigma} r)}$ reflects the direct matching discrepancy for $\boldsymbol{\alpha}_i$. It grows quickly with the number of latent variables, which coincides with the results in the nearest neighbors matching literature (e.g., Gyorfi-et-al_2002_bookchapter). The second term $a_n^{1/\underline{\varsigma}}$ arises from the existence of the idiosyncratic error $\bm{u}_i$. In the three examples described above, the distance is defined based on certain averages across different features, and thus the impact of $\bm{u}_i$ vanishes as $p$ grows large.
Note that if $\boldsymbol{\alpha}_i$'s were observed, matching could be directly implemented on it with the number of matches $K$ fixed. In this paper, however, $\boldsymbol{\alpha}_i$ is unobservable, and matching can only be done on their noisy measurements, leading to the indirect matching discrepancy characterized by the second term above. Using a fixed (or small) number of nearest neighbors is unable to further reduce bias and thus is not recommended in this scenario.
Theorem (ref) generalizes the existing results in the literature that relies on specific metrics and provides a way to precisely quantify the indirect matching discrepancy. For example, Zhang-Levina-Zhu_2017_BIMA proposes the pseudo-max distance for neighborhood smoothing in the graphon estimation context, but they provide no results regarding the distance in the latent variables. Moreover, if the conditions specified in Theorem (ref) hold and $\eta_l(\boldsymbol{\alpha}_i)=\eta_l(\alpha_i, \varpi_l)$ for scalar latent variables $\alpha_i$ and $\varpi_l$, Theorem (ref) implies that the neighborhood smoothing (local average) estimator of the factor structure can achieve a sup-norm convergence rate of order $O(n^{-\frac{1}{3}})$, up to $\log n$ terms, which improves upon the $L^2$-type convergence rate of order $O(n^{-1/4})$, up to $\log n$ terms, given in Zhang-Levina-Zhu_2017_BIMA. See more detailed discussion about the uniform convergence rate in Section (ref).
Now, we proceed to discuss the properties of local principal component analysis. Throughout this subsection, we assume a set $\mathcal{N}_i$ of $K$ nearest neighbors for each $i\in[n]$ has been obtained. Recall that PCA is applied to the submatrix $\bm{X}_{\langle i \rangle}$ of $\bm{X}$ formed by a subset of observed features indexed by $\ddagger$ for the $K$ nearest neighbors of each unit $i$. Define the constants $\delta_n=(K\wedge p)^{1/2}/\sqrt{\log (n\vee p)}$ and $h_n=(K/ n)^{\bar{\varsigma}/(\underline{\varsigma} r)}$.
We need some regularity conditions on the local approximation of $\bm{H}_{\langle i \rangle}$. Formally, we consider the following $L^2$-approximation
where $\bm{F}_{ \langle i \rangle}= \mathbb{E}^\ddagger[\bm{H}_{\langle i \rangle}\boldsymbol{\Lambda}_{\langle i \rangle}] \mathbb{E}^\ddagger[\boldsymbol{\Lambda}_{\langle i \rangle}'\boldsymbol{\Lambda}_{\langle i \rangle}]^{-1}$ and $\mathbb{E}^\ddagger$ denotes the expectation operator conditional on $\bm{X}^\dagger$. Accordingly, $\boldsymbol{\Xi}_{\langle i \rangle}$ should be understood as a matrix of $L^2$-projection errors. We introduce a diagonal (scaling) matrix $\boldsymbol{\Upsilon}_{\langle i \rangle}=\operatorname*{diag}\{\upsilon_{1,\langle i \rangle}, \cdots, \upsilon_{\mathsf{d}_i,\langle i \rangle}\}$, denoting the possibly heterogeneous strength of local factors for the neighborhood of unit $i$. Without loss of generality, we assume $\upsilon_{1, \langle i \rangle}\geq \upsilon_{2,\langle i \rangle}\geq\cdots\geq \upsilon_{\mathsf{d}_i,\langle i \rangle}$.
Among the three conditions, (a) and (b) are usually mild and similar to those required for sieve approximation in the nonparametric regression literature. Typically, they can be verified by properly choosing an approximation basis. For example, let $\boldsymbol{\alpha}\in\mathcal{A}\mapsto \boldsymbol{\lambda}_{\langle i \rangle}(\boldsymbol{\alpha}):=(\lambda_1(\boldsymbol{\alpha}), \cdots, \lambda_{\mathsf{d}_{\langle i \rangle}}(\boldsymbol{\alpha}))'$ be an $r$-variate monomial basis of degree no greater than $m-1$ centered at $\boldsymbol{\alpha}_i$ (including the constant term), with a typical element given by $(\boldsymbol{\alpha}-\boldsymbol{\alpha}_i)^{\bm{q}}$ for some $\bm{q}=(q_1, \cdots, q_r)'$ such that $\sum_{j=1}^rq_j\leq m-1$. The length of this basis is $\mathsf{d}_i=\binom{r+m-1}{r}$. Define $\boldsymbol{\Lambda}_{\langle i \rangle}=(\boldsymbol{\lambda}_{\langle i \rangle}(\boldsymbol{\alpha}_{j_1(i)}), \cdots, \boldsymbol{\lambda}_{\langle i \rangle}(\boldsymbol{\alpha}_{j_K(i)}))'$. Accordingly, we can let $\boldsymbol{\Upsilon}_{\langle i \rangle}=\operatorname*{diag}(1, h_n\bm{I}_{\ell_1}, \cdots, h_n ^{m-1}\bm{I}_{\ell_{m-1}})$ with $\ell_j=\binom{r+j-1}{j}$ for $j\in[m-1]$, which properly normalizes the approximation basis functions of different order. Then, parts (a) and (b) can be verified under the regularity conditions on $\boldsymbol{\eta}(\cdot)$ specified in Assumption (ref).
The requirement on $\upsilon_{j,\langle i \rangle}/\upsilon_{j+1,\langle i \rangle}$ in part (a) merely formalizes the possibility that the local approximation terms may be of different strength in an asymptotic sense. Consider the example of the monomial basis discussed above. The first element corresponds to the constant term whose magnitude is (asymptotically) greater than the next $r$ elements that corresponds to local polynomials of degree one. The other elements correspond to local polynomials of the second or higher order whose magnitude is even smaller. On the other hand, the condition $\delta_n^{-1}/\upsilon_{\mathsf{d}_i,\langle i \rangle}=o(1)$ in part (b) is key for the local factors to be consistently estimable. Intuitively, $\upsilon_{\mathsf{d}_i, \langle i \rangle}$ determines the signal strength of the “weakest” factors one desires to extract, which has to be stronger than the strength $\delta_n^{-1}$ of the noise.
Part (c) should be taken with some extra care. Such conditions, usually termed non-degeneracy of factors, are common in linear factor analysis. In the nonlinear setting, however, local factors $\bm{F}_{\langle i \rangle}$ have a specific meaning: they are transformations of derivatives of $\boldsymbol{\eta}$. The non-degeneracy condition (c) indeed reflects to what extent the functions $\{\eta_l: l\in\mathcal{R}^\ddagger\}$ are nonlinear at each point in the support. In general, the degree of nonlinearity may be heterogeneous across the evaluation points, and condition (c) could fail with a universal choice of the number of local factors such as $\mathsf{d}_i=\binom{r+m-1}{r}$ discussed above. Therefore, it is important in this context to allow $\mathsf{d}_i$ to vary across $i$ (different local neighborhoods), which substantially weakens or rationalizes this nonlinearity requirement. To gain more intuition, we discuss two examples below.
In Appendix (ref) we provide a general result about the verification of Assumption (ref) under primitive conditions, and Section SA-3 in the SA gives further discussion with concrete examples.
Now, we present the uniform convergence properties of the estimated local factors and loadings, the second main result of this paper.
Recall that Theorem (ref) showed that the distance between $\boldsymbol{\alpha}_i$ and its nearest neighbors $\{\boldsymbol{\alpha}_{j_k(i)}: 1\leq k\leq K\}$ is diminishing as $n$ and $p$ diverge. Theorem (ref) above further shows that the loading matrix $\boldsymbol{\Lambda}_{\langle i \rangle}$, as the “approximation basis”, can be consistently estimated up to a rotation, which indeed characterizes the relationship of different units within each local neighborhood. Accordingly, the factor matrix $\bm{F}_{\langle i \rangle}$, as (transformations of) derivatives of $\boldsymbol{\eta}(\cdot)$, can also be consistently estimated, which characterizes the local nonlinearity of the latent surface.
The matrix convergence result above is in terms of the sup-norm, which also holds uniformly over the local neighborhoods indexed by $\langle i \rangle$. The estimation error of both $\widehat{\bm{F}}_{\langle i \rangle}$ and $\widehat{\boldsymbol{\Lambda}}_{\langle i \rangle}$ consists of two parts. The first term in each upper bound, i.e., $\delta_n^{-1}\upsilon_{\ell, \langle i \rangle}^{-1}$ and $\delta_{n}^{-1}$, reflects the estimation variance. Since latent variables $\boldsymbol{\alpha}_i$'s are not observed, the variability of the local factor and loading estimators relies on the sizes of both dimensions $K$ and $p$, which is akin to the results in classical linear factor analysis Bai_2003_ECMA. The second term in each upper bound, i.e., $h_n^m\upsilon_{\ell, \langle i \rangle}^{-1}$ and $h_n^m$, arises from the smoothing bias and is not present in linear factor analysis.
By Assumption (ref)(a), the local factors in $\bm{F}_{ \langle i \rangle}$ may be of heterogeneous strength, reflected by the magnitude of the associated loadings. Consequently, the estimation error of $\widehat{\bm{F}}_{\langle i \rangle}$ involves a penalty factor $\upsilon_{\ell, \langle i \rangle}^{-1}$ that is inversely related to the strength of different factors. This is in line with the intuition: higher-order approximation terms are weaker signals about the latent structure, which is less precisely estimated. As emphasized before, the rate condition $\delta_{n}^{-1}/\upsilon_{\mathsf{d}_i,\langle i \rangle}\rightarrow\infty$ is key to ensure the “weakest” local factors can still be differentiated from the remainder in Equation (ref).
To get some sense of the rate conditions required, consider the simple case where $p\asymp n$ and $\boldsymbol{\Upsilon}=\operatorname*{diag}(1, h\bm{I}_{\ell_1}, \cdots, h^{m-1}\bm{I}_{\ell_{m-1}})$. Assumption (ref)(b) needs $(n/K)^{\frac{2m-2}{r}}=o(K/\log n)$, and the additional rate restriction imposed in the theorem can be simplified to $n^{\frac{4}{\nu}}\lesssim K/\log n$. If we let $K=n^{A}$, $A>\max\{\frac{4}{\nu}, \frac{2m-2}{2m-2+r}\}$ suffices. In particular, if $\nu$ is sufficiently large, this restriction can be satisfied by setting, for example, $K\asymp n^{\frac{2m}{2m+r}}$, which coincides with the MSE-optimal choices of tuning parameters in nonparametric regression.
Finally, using Theorem (ref), we immediately have the following corollary.
This corollary shows that the latent nonlinear factor component can be consistently estimated, and this error bound holds uniformly over all units and features. Note that if $\boldsymbol{\alpha}_i$'s were observed, the natural alternative to estimate the heterogeneous functions $\eta_l$'s (and thus $\bm{H}$) would be the cross-sectional nonparametric regression of $x_{il}$ on $\boldsymbol{\alpha}_i$ for each $l\in[p]$, whose optimal uniform convergence rate is $O(n^{-\frac{m}{2m+r}})$, up to $\log n$ terms stone1982optimal. Corollary (ref) shows that this optimal rate can be attained by the local PCA estimator when $K\lesssim p$, $\bar{\varsigma}=\underline{\varsigma}$, and an optimal number of nearest neighbors is used to balance the two terms in the above error bound. This result appears to be new to the literature, to the best of our knowledge.
Also note that our setup in this paper allows the latent functions $\eta_l$'s to be heterogeneous across $l$ in general. When $\eta_l(\boldsymbol{\alpha}_i)$ admits a special structure such as $\eta_l(\boldsymbol{\alpha}_i)=\eta(\boldsymbol{\alpha}_i,\bm\varpi_l)$ for some bivariate function $\eta$ and some random variables $\bm\varpi_l$, the optimal convergence rate for estimating the global function $\eta$ may be different than the optimal rate for the cross-sectional nonparametric regression described above, depending on the smoothness of $\eta$ and the dimensions of the arguments $\boldsymbol{\alpha}_i$ and $\bm\varpi_l$. See Gao-Lu-Zhou_2015_AoS for a discussion of the $L^2$-type minimax optimal rate in the special case of graphon estimation.
So far we have focused on the case where the data matrix $\bm{X}$ is complete, but in many applications $\bm{X}$ may have missing entries. The proposed local PCA method is robust with respect to mild missing value issues and thus can still be applied to some policy evaluation problems such as synthetic controls. Detailed discussion is deferred to Section (ref).
The analysis so far has focused on Equation (ref), which assumes $\bm{x}_i$ takes a purely nonlinear factor structure. However, high-rank components may exist in $\bm{x}_i$, and it is the residuals that take a possibly nonlinear factor structure, as described by Equation (ref) below:
where $\bm{W}_i=(\bm{w}_{i,1}, \cdots, \bm{w}_{i,q})\in\mathbb{R}^{p\times q}$ is a matrix of observed covariates. It can be viewed as a linear regression of $\bm{x}_i$ on $q$ regressors $\bm{w}_{i,1}, \cdots, \bm{w}_{i,q}$ with possibly nonlinear fixed effects $\boldsymbol{\eta}(\boldsymbol{\alpha}_i)$. In this context $\bm{W}_i$ needs to have sufficiently high-rank variation, otherwise they will be too collinear with the unspecified low-rank component $\boldsymbol{\eta}(\boldsymbol{\alpha}_i)$ and $\boldsymbol{\vartheta}$ cannot be identified. This is similar to the identification condition for panel data regression with interactive fixed effects.
The main analysis of this paper can be applied once a consistent estimator of $\boldsymbol{\vartheta}$ is available. It can be obtained using the idea of partially linear regression. Specifically, we propose the following procedure:
Under additional regularity conditions on $\bm{W}_i$, it can be shown that $\widehat{\boldsymbol{\vartheta}}$ converges to $\boldsymbol{\vartheta}$ sufficiently fast and the main results established previously still hold for $\mathcal{N}_i$, $\widehat{\bm{F}}_{\langle i \rangle}$ and $\widehat{\boldsymbol{\Lambda}}_{\langle i \rangle}$. Formal analysis is available in Section SA-2 of the SA and is omitted here to conserve space.
We conduct a Monte Carlo investigation of the finite sample performance of the proposed method. We consider two nonlinear factor models for continuous data and one for discrete data:
The latent variables $\alpha_i$ and $\varpi_l$ follow the uniform distributions on $[0, 1]$.
We consider $2\,000$ simulations with $n=p=1\,000$. To implement the proposed method, we use the first half of rows of $\bm{X}$ for $K$-NN matching and the second half for principal component analysis. We take the pseudo-max distance and set the number of nearest neighbors $K=\mathsf{c}\times n^{2/3}$ for $\mathsf{c}=0.5, 1$ and $1.5$. To avoid extracting too weak local factors, we let the number of local principal components $\mathsf{d}_i=2$ if $\upsilon_{2, \langle i \rangle}/\upsilon_{3,\langle i \rangle}\geq \log\log K$ and $\mathsf{d}_i=1$ otherwise, for each $i\in[n]$. We also compare the proposed method (LPCA) with a natural alternative---global principal component analysis (GPCA). Specifically, we use the entire matrix $\bm{X}$ to extract principal components and form the mean prediction for each entry. The number of factors is determined by using the eigenvalue ratio test Ahn_2013_ECMA based on doubly demeaned data. We report the following results for each method: (1) maximum absolute error (MAE), i.e., $\max_{i\in[n], l\in\mathcal{C}^\ddagger}|\widehat{\eta}_{il}-\eta_{il}|$; and (2) the prediction error for three missing entries in the last row of $\bm{X}$. For (2), in each simulated dataset we take the three units with the values of $\alpha_i$ equal to the $0.1$-, $0.5$-, and $0.9$-quantiles of the sample respectively, and then replace their corresponding values in the last row of $\bm{X}$ with zeros (“missing values"). This operation mimics the data pattern in some policy evaluation settings: a unit whose individual feature is at the low, medium or high level relative to the whole distribution gets treated in the last period and thus has a missing value in the last row of $\bm{X}$. Finally, all these measurements of performance are averaged across $2\,000$ repetitions.
\FloatBarrier
\FloatBarrier
The proposed method performs well across the three settings. Compared with GPCA, LPCA usually has smaller or at least comparable maximum absolute errors and prediction errors for the three missing entries. It can be seen that such improvement is more pronounced when the nonlinearity of $\eta_l$ is more severe (e.g., Model 1) and a sufficiently small $K$ is used. The results are relatively stable across different choices of $K$ except that in Model 1---a highly nonlinear one---using a very large $K$ leads to much increase in MAE. In the SA we report additional simulation evidence for datasets with smaller $p$ and for LPCA with a larger fraction of rows for $K$-NN matching (and thus a smaller fraction for obtaining the principal components). The results are largely close to the main findings above.
The proposed method is readily applicable to certain matrix completion problems. For instance, in the classical synthetic control design Abadie_2020_JEL, the researcher is interested in the effect of a policy that only affects one single unit, and all other units in the data remain untreated throughout the observation period. Formally, we can define the potential outcome under the treatment $y_{il}(1)$ and that in the absence of the treatment $y_{il}(0)$ for each unit $i\in[n]$ and time $l\in[p]$. The observed outcome $y_{il}=d_{il}y_{il}(1)+(1-d_{il})y_{il}(0)$ where $d_{il}=\mathds{1}(i=1,l>p_0)$ for some $p_0<p$. That is, a policy intervention was implemented at time $p_0+1$ and affected the first unit only. The key task in this problem is to predict the missing counterfactual outcome $y_{1l}(0)$ for the treated unit in each post-treatment period $l> p_0$. Usually, the matrix of the potential outcomes $\bm{Y}^N$, whose typical entry is given by $y_{il}(0)$, is assumed to admit a linear factor structure, which justifies the key idea of synthetic controls that uses a linear combination of untreated units to predict the counterfactual of the treated.
The proposed local PCA method relaxes the linearity requirement by allowing for a possibly nonlinear factor structure. If the previous procedure could be applied to the matrix $\bm{Y}^N$, we would immediately have an estimate of the mean of the counterfactual outcome of the treated in the post-treatment period, and Corollary (ref) would apply. However, several entries at the “bottom-left” corners of $\bm{Y}^N$ are unobserved. Let $$x_{il}=y_{il}\mathds{1}(i\neq 1 \text{ or } l\leq p_0).$$ Thus, the matrix $\bm{X}$ is identical to $\bm{Y}^N$ except that the missing outcomes of the treated in $\bm{Y}^N$ are replaced with zeros. Then we apply the proposed procedure to the observed matrix $\bm{X}$. Assume the number of post-treatment periods $p-p_0$ is fixed, which is common in the classical synthetic control analysis. Thus $\bm{X}$ can be viewed as a slightly perturbed version of $\bm{Y}^N$. The following theorem shows that the proposed mean estimator $\widehat{\bm{F}}_{\langle i \rangle}\widehat{\boldsymbol{\Lambda}}_{\langle i \rangle}'$ is robust to such “small” perturbation.
This result relies on the robustness of leading singular vectors with respect to small perturbation of $\bm{Y}^N$. Similar results still hold for missing data patterns other than synthetic controls. The key requirement is that there are not many missing entries in each row and in each column. See the proof of Theorem (ref) in Appendix (ref) for more details. When missingness leads to relatively large perturbation of $\bm{Y}^N$, local PCA is not robust and thus may not be directly applied. Feng_2023_wp discusses the scenario where there are many missing entries in one or a few columns (or rows) and show that the extracted local factors and loadings obtained from the complete part of the matrix can be used for prediction of missing entries.
To illustrate the proposed method, we reanalyze the effect of 2012 Kansas tax cuts on economic growth (see Rickman-Wang_2018_RSUE for more details). The second quarter of 2012 (2012Q2)---when the governor Brownback signed the tax cut bill into law---is used as the starting time of the policy. As described above, this can be viewed as a synthetic control problem where the goal is to predict the counterfactual economic performance of Kansas after 2012Q2 had the tax cut policy not been implemented.
We have a data matrix on quarterly GDP per capita of 50 states ($n=50$) from 1990Q1 to 2016Q1 ($p=105$). We take the first difference of the original data and analyze the effect of tax cuts on GDP per capita growth rates. The growth rates of Kansas after 2012Q2 are replaced with missing values. The same strategy as described in Section (ref) is used to implement local PCA: we take the pseudo-max distance, set $K=n^{2/3}\approx 14$ and choose the number of local factors based on the strengths of leading singular values. The first $40$ periods are used for $K$-NN matching, and the remaining periods are used for PCA. For comparison, we also implement the classical synthetic control (SC) method that relies on a simplex constraint on weights.
The result is shown in Figure (ref) below. It turns out that in 9 out of 16 post-treatment periods the observed GDP growth rate sequence is below the prediction from LPCA. Note that LPCA estimates the mean value of GDP growth, and thus the diffrerence between the two growth paths may arise from both estimation error and an idiosyncratic noise in each period. To reduce the impact of the idiosyncratic noise, we also compute the average growth rate over the entire post-treatment periods for the counterfactual Kansas, which is $0.53$ percentage points higher than that of the observed Kansas. By constrast, SC yields a predicted average post-treatment growth rate that is $0.19$ percentage points lower than the observed Kansas, which is not plausible given the poor fit of SC in the pre-treatment period.
To see the impact of the tax cut policy on the level of GDP per capita, we take the GDP per capita of Kansas in 2012Q1 as the initial value and translate the predicted growth rates into a counterfactual GDP per capita trajectory for Kansas. In this case, we compare LPCA with SC based on the data of GDP per capita levels rather than that using the growth rates data, which is also common in the literature. The results are given in Figure (ref). In this case both methods show that overall the tax cut policy has a pronounced negative shock on the GDP growth trajectory of Kansas, but the effect given by SC is relatively small in magnitude, especially in later post-treatment periods.
\FloatBarrier
\FloatBarrier
This paper studied optimal estimation of large-dimensional nonlinear factor models. The observed variables were modeled as possibly nonlinear functions of some latent variables with the functional forms unspecified. A local principal component analysis method that combines $K$ nearest neighbors matching and principal component analysis was proposed to estimate the nonlinear factor structure and recover the information on latent variables and latent functions. The large-sample properties of the proposed estimators were established, including a sharp bound on the matching discrepancy of nearest neighbors under generic distance functions, the sup-norm error bounds for estimated local factors and loadings, and the uniform convergence rate for the estimator of the factor structure. The method can also be applied to some matrix completion problems such as synthetic controls, which does not require the mean matrix to be exactly low-rank.