Causal Inference in Possibly Nonlinear Factor Models

Causal Inference in Possibly Nonlinear Factor Models

Keywords: causal inference, latent confounders, nonlinear factor model, low-rank method, heterogeneous treatment effects, doubly-robust estimator, high-dimensional data

\thispagestyle{empty}

\onehalfspacing \setcounter{page}{1}

\pagestyle{plain}

Introduction

Understanding effects of policy interventions is central in many disciplines Heckman-Vytlacil_2007_HandbookChapter,Angrist-Pischke_2008_book,Imbens-Rubin_2015_book,Abadie-Cattaneo_2018_ARE,Hernan-Robins-2020_book. When observational data are used, researchers usually confront the challenge that the treatment is nonrandomly assigned based on some characteristics that are not directly observed. The confounding effects of these variables (confounders) make it difficult to uncover the true causal relation between the outcome and the treatment. Commonly used econometric methods that assume selection on observables are inappropriate in this situation. This paper proposes a treatment effects model in which a large set of observed covariates, as the noisy measurements of the underlying confounders, are available. The key assumption is that the observed measurements and unobserved confounders are linked via an unknown, possibly nonlinear factor model. The former, though not affecting the potential outcome and the treatment assignment directly, provide information on the latter, thus making it possible to resolve the confounding issue. Exploiting this underlying factor structure, I develop a novel inference method for counterfactual analysis, which can be used in many applications such as synthetic control designs, recommender systems, diffusion index forecasts, and network analysis.

As an example, consider the effect of a scholarship on the academic performance of newly admitted college students. One may be concerned about the confounding effect of the unobserved precollege ability, since it may correlate with both a student's likelihood of getting a scholarship and her future academic performance. If the researcher is able to observe the same student taking multiple tests in different subjects or time periods at the precollege stage, these past test scores may play the role of the noisy measurements of the unobserved ability. The nonlinear factor structure allows for a flexible latent relationship between ability and test outcomes, which may vary across subjects or time in a complex way.

The key building block (one of the main contributions of this paper) is a carefully designed local principal subspace approximation procedure that allows for flexible functional forms in the factor model. The procedure begins with $K$-nearest neighbors ($K$-NN) matching for each unit on the observed noisy measurements. The number of nearest neighbors, $K$, diverges as the sample size increases, which differs from other matching techniques that use only a fixed number of matches Abadie-Imbens_2006_ECMA. Within each local neighborhood formed by the $K$ matches, the underlying possibly nonlinear factor structure is approximated by a linear factor structure and can be estimated using principal component analysis (PCA). Theoretical properties of this local PCA method are derived in this context. Under mild conditions on the unknown factor structure, the nearest neighbors and estimated factor loadings characterize the unobserved confounders and can be used to match comparable units in the subsequent treatment effects analysis.

Building upon this observation, I develop a novel inference procedure for a large class of causal parameters. It has three appealing features. First, as a dimension reduction technique, the proposed method allows users to obtain low-dimensional information on latent confounders from large-dimensional noisy measurements. It only requires some but not all measurements to be informative about latent confounders, and it is unnecessary to know their identities a priori (see Remark (ref) below). Second, the proposed method does not impose a functional form assumption on the relationship between latent confounders and noisy measurements, thus making the final inference more robust. In particular, the nonlinearity of this relationship is allowed but not assumed, and the classical linear factor model can be covered as a special case. Third, the output of local PCA can be readily used as input to many classical econometric estimation such as local polynomial kernel regression Fan-Gijbels_1996_book. Thus, the local PCA, as a useful pre-processing step, can be combined with many other econometric applications and is of independent interest.

To fix ideas, suppose that the treatment occurs at some point in time (staggered adoption can also be allowed as described in Section SA-4.1 of the Supplemental Appendix). The assignment is correlated with unit-specific latent features $\boldsymbol{\alpha}_i\in\mathbb{R}^{\mathsf{d}_\alpha}$ for $1\leq i\leq n$. The untreated outcome, observed in $T_0$ periods prior to the treatment, is a time-heterogeneous, possibly nonlinear function of $\boldsymbol{\alpha}_i$, say, $\eta_t(\boldsymbol{\alpha}_i)$, plus some noise, which is usually termed a (possibly) nonlinear factor structure Yalcin-Amemiya_2001_SS. The latent features $\{\boldsymbol{\alpha}_i\}_{i=1}^n$ play the role of confounders in this context and are akin to fixed effects in the panel data literature Arellano_2003_book. Geometrically, the set of latent functions $\{\eta_t(\cdot)\}_{t=1}^{T_0}$ generates a low-dimensional subspace embedded in a high-dimensional space when the number of pre-treatment periods $T_0$ is large but the number of latent confounders $\mathsf{d}_{\alpha}$ is small. Suppose that different values of the latent confounders induce non-negligible differences in outcomes in many pre-treatment periods. In this case, the $K$ nearest neighbors of each unit as appropriately measured by the observed outcome should also be close in terms of the latent confounders. Such nearest neighbors form a local neighborhood for a unit and are approximately lying on a subspace that can be characterized by a linear combination of basis functions of $\boldsymbol{\alpha}_i$ and estimated by local PCA. Consequently, the underlying nonlinear factor structure is locally approximated by principal subspaces, up to errors governed by the number of nearest neighbors and the number of local principal components extracted. The availability of many repeated measurements of the latent confounders (pre-treatment outcomes in this example) is crucial for the validity of this approximation. It affects the matching discrepancy of nearest neighbors and the estimation precision of local principal components.

As in linear factor models Bai_2003_ECMA, the values of the latent variables cannot be exactly recovered without additional normalizations. Nevertheless, the $K$ nearest neighbors and local principal components from the above approximation procedure suffice to control for the latent confounders in the subsequent analysis. In fact, they can be readily used as inputs in commonly used nonparametric kernel regression. The local region used in the estimation is defined by nearest neighbors, and the extracted local principal components play the role of generated regressors that provide further approximation to unknown conditional expectation functions of interest. The number of nearest neighbors implicitly governs the bandwidth of the regression, which determines the consistency of final estimators and is the main tuning parameter in the proposed estimation procedure. By contrast, the number of local principal components extracted is analogous to the order of the basis in local polynomial regression and is often fixed in practice.

In the causal inference context, I propose using various local regression methods to estimate conditional means of potential outcomes and conditional treatment probabilities (generalized propensity scores), which form the basis of regression imputation and propensity score weighting estimators. In contrast with standard nonparametric regression analysis, the conditioning variables in this scenario are indirectly obtained from the observed measurements, and the noise in their factor structure restricts one's ability to select a bandwidth. Using a small or fixed number of nearest neighbors does not necessarily lead to a small bandwidth and thus is not helpful for further bias reduction. Consequently, the possibly large smoothing bias of the nonparametric ingredients may render the final inference on causal parameters invalid. To deal with this issue, I follow the Neyman-orthogonalization strategy that has been extensively applied in the recent double/debiased machine learning literature Belloni-Chernozhukov-Hansen_2014_RES,Farrell_2015_JoE,Chernozhukov-et-al_2018_EJ,Chernozhukov-et-al_2020_wp. In treatment effects models, the widely used doubly-robust scores Robins-Rotnitzky_1995_JASA,Cattaneo_2010_JoE are estimating equations constructed based on the efficient influence function and are automatically Neyman orthogonal. Owing to this property, valid inference can be conducted that only requires mild restrictions on the local principal subspace approximation. In the proposed estimation procedure, the set of observed measurements is split for the purpose of local PCA, whereas the final inference stage does not require sample splitting, which differs from other debiased learning methods based on cross-fitting.

Based on the ideas above, I develop a novel estimation and inference procedure for treatment effects analysis. The results cover a large class of estimands, including counterfactual distributions and functionals thereof, and provide the basis for analyzing many causal quantities of interest such as average, quantile, and distributional treatment effects. Moreover, the underlying local PCA method has broad applicability, providing a new tool for the analysis of panel and network data and other data with similar structures. Some useful results are established. First, under mild geometric conditions on the underlying subspace, a sharp bound is derived for the implicit discrepancy of latent variables induced by nearest neighbors matching. Second, uniform convergence of the estimated local factors and loadings is established, taking into account the possibly heterogeneous strength of factors due to the nonlinearity of the model. These results can be applied to study, for example, linear regression models with nonlinear fixed effects. Detailed technical results are available in Section SA-2 of the Supplemental Appendix (SA), which is of independent interest. Typical applications, including staggered adoption, recommender systems, inference with measurement error, diffusion index forecasts, and network analysis, are discussed in Section SA-4 of the SA.

The paper is organized as follows. The rest of this section discusses the related literature. In Section (ref), I set up a multi-valued treatment effects model and describe the nonlinear factor structure of the large-dimensional measurements of latent confounders. Section (ref) gives a detailed description of the entire estimation procedure, accompanied by a step-by-step empirical illustration using the data of Acemoglu-et-al_2016_JFE. Section (ref) presents the main theoretical results and some Monte Carlo evidence. Section (ref) discusses uniform inference on counterfactual distributions as well as other useful extensions. Section (ref) concludes. The Supplemental Appendix contains all theoretical proofs, additional technical results, methodological discussions, and typical applications. Replications of the simulation study and empirical illustration are available at \url{https://github.com/yingjieum/replication-Feng_2021}.

Related Literature

This paper contributes to several strands of literature. First, since the observed covariates may be viewed as an array of noisy measurements of the latent confounders, my theoretical framework is closely related to nonlinear models with measurement errors. Much effort has been devoted to the identification of such models (see Schennach_2016_ARE for a review). For example, factor models can be utilized to construct repeated measurements of unobserved variables, which allows for the identification of their distribution under suitable normalizations. A general treatment following this strategy is available in Cunha-Heckman-Schennach_2010_ECMA, using and extending results in Hu-Schennach_2008_ECMA. My paper takes a different route. A large-dimensional nonlinear factor model is exploited to directly extract the geometric relation among different units in terms of the latent variables, which is then used to control for their confounding effects in the treatment effects analysis. Some measurements are allowed to be uninformative about the latent confouders, and to identify the causal effect of interest, it is unnecessary to recover the exact values (or distributions) of latent confounders. Conceptually, the extracted information from the observables plays a similar role as a control function, conditional on which the treatment assignment is no longer confounded. See Wooldridge_2015_JHR for a review of control function methods in econometrics and Altonji-Mansfield_2018_AER for an application of the idea of using the transformation of observables to control for unobservables in the context of estimating group effects.

Second, my study builds on and extends some results on large-dimensional factor analysis and panel regression with fixed effects Bai_2009_ECMA,Bai-Wang_2016_ARE,Wang-Fan_2017_AoS. In particular, my proposed method generalizes the idea of linear factor-augmented prediction---sometimes referred to as diffusion index forecasts in macroeconometrics Stock-Watson_2002_JASA,Bai-Ng_2006_ECMA---to nonlinear factor models. The differences are that the proposed method does not rely on a linear factor structure and that my primary goal is inference on treatment effects or other causal quantities rather than the prediction of outcomes. Recent work by Chernozhukov-et-al_2019_wp develops an inference method for linear panel regression models, where both slopes and intercepts have linear factor structures. By contrast, my paper focuses on a heterogeneous treatment effects model, and the panel-like structure of the observed measurements is exploited to control for latent confounders rather than being of direct interest. Section (ref) below extends the main analysis by specifying a more general structure for observed measurements, which can be viewed as a linear panel regression model with high-rank regressors. The slope coefficients are homogeneous across both dimensions, whereas the intercept admits a possibly nonlinear factor structure. Another recent study by Bonhomme-Lamadon-Manresa_2021_wp develops two-step grouped fixed-effects estimators that discretize latent heterogeneity by $K$-means clustering. By contrast, my paper relies on a general identification condition (see Remark (ref) for details) and uses local principal subspace approximation strategy. The proposed method can achieve more flexible approximation of smooth functions of latent features, and an intermediate result (Theorem (ref)) also characterizes the uniform convergence rates of nonparametric estimators of individual-specific features such as conditional expectation of an outcome given an individual's latent confounders.

Third, the idea of local approximation of nonlinear subspaces 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. For example, the local tangent space alignment (LTSA) algorithm of Zhang-Zha_2004_SIAM exploits the idea of local PCA, and the local linear embedding (LLE) algorithm of Roweis-Saul_2000_Science aims at learning local self-reproducing weights. Other recent advances include Peng-Lu-Wang_2015_NN,Zhang-et-al_2015_PR,Arias-Lerman-Zhang_2017_JMLR, among others. These methods are used to construct global nonlinear subspaces that preserve the local geometry of the data for the purpose of classification, clustering or data visualization. Unlike these studies, this paper focuses on estimation and inference of causal parameters in the treatment effects model rather than recovering the latent nonlinear subspaces. Also, statistical properties of the proposed local PCA procedure are formally characterized and is of independent interest for other applications.

Finally, this study contributes to the existing literature on causal inference and program evaluation (see Abadie-Cattaneo_2018_ARE for a review). For example, it is connected with the fast-growing literature on synthetic control (see Abadie_2020_JEL and references therein) and staggered adoption designs Athey-Imbens_2018_wp. The classical synthetic control method and many variants thereof are often motivated by assuming a linear factor structure for the pre-treatment data. By contrast, my paper allows for a possibly nonlinear factor structure and does not rely on the strong assumption of linear factor models. Using the geometric relation among different units characterized by nearest neighbors and local factor loadings, I derive formal large-sample properties of the proposed estimators under mild side conditions.

Treatment Effects Model with Latent Variables

Suppose that a random sample $\{(y_{i}, s_{i}, \bm{x}_{i}, \bm{w}_i, \bm{z}_i)\}_{i=1}^n$ is available, where $y_{i}\in\mathbb{R}$ is the outcome of interest, $s_i\in\mathcal{J}=\{0,\cdots, J\}$ denotes the treatment status, and $\bm{x}_{i}\in\mathbb{R}^{T}$, $\bm{w}_i\in\mathbb{R}^{T\mathsf{d}_w}$ and $\bm{z}_i\in\mathbb{R}^{\mathsf{d}_z}$ are vectors of covariates. $\bm{x}_i$, $\bm{w}_i$ and $\bm{z}_i$ play different roles in later analysis: $\bm{x}_i$ and $\bm{w}_i$ are used to obtain information on a vector of unobserved confounders $\boldsymbol{\alpha}_i\in\mathbb{R}^{\mathsf{d}_{\alpha}}$, whereas $\bm{z}_i$ itself is a set of observed confounders that can be controlled for directly. Some covariates may be used for the two purposes simultaneously, and thus $\bm{z}_i$ may share some variables in common with $\bm{x}_i$ and $\bm{w}_i$. The asymptotic theory in this paper is developed assuming $n$ and $T$ simultaneously increase to infinity whereas $\mathsf{d}_w$, $\mathsf{d}_z$ and $\mathsf{d}_{\alpha}$ are fixed.

I follow the standard potential outcomes framework. Let $y_i(\jmath)$ denote the potential outcome of unit $i$ at treatment level $\jmath\in\mathcal{J}$. Construct an indicator variable $d_{i}(\jmath)=\mathds{1}(s_{i}=\jmath)$ for each $\jmath\in\mathcal{J}$. The observed outcome can be written as $y_{i}=\sum_{\jmath=0}^{J}d_{i}(\jmath)y_{i}(\jmath)$. Many interesting parameters can be defined in this framework, and the key challenge is to overcome the missing data issue. For example, when $s_i$ is binary, i.e., $s_i\in\{0,1\}$, the identification of average treatment effects on the treated (ATT) relies on $\mathbb{E}[y_i(0)|s_i=1]$, but $y_i(0)$ is unobservable for the treated group. This hurdle is often overcome by imposing an unconfoundedness condition so that the treatment assignment becomes independent of potential outcomes after conditioning on a set of observed covariates. By contrast, this paper assumes that \[ y_i(\jmath)\protect\mathpalette{\protect\independenT}{\perp} d_i(\jmath')\,|\,\bm{z}_i, \boldsymbol{\alpha}_i,\quad \forall \jmath,\jmath'\in\mathcal{J}. \] Recall that $\bm{z}_i$ is observed, but $\boldsymbol{\alpha}_i$ is unobserved and thus cannot be directly controlled for. As described later in Section (ref), the noisy measurements $\bm{x}_i$ contain information on $\boldsymbol{\alpha}_i$ and help restore unconfoundedness in the treatment effects analysis.

For each treatment level $\jmath\in\mathcal{J}$, the outcome of interest is characterized by a possibly nonlinear, reduced-form model:

where $\varsigma_{i,\jmath}$ is the conditional expectation of the potential outcome at treatment level $\jmath$ given the observed $\bm{z}_i$ and unobserved $\boldsymbol{\alpha}_i$, and $\psi_{\mathsf{y}}^{-1}(\cdot):\mathbb{R}\mapsto\mathbb{R}$ is a (known) link function associated with the outcome equation.

On the other hand, introduce a (known) link function $\boldsymbol{\psi}^{-1}_{\mathsf{s}}(\cdot):(0,1)^{J+1}\mapsto\mathbb{R}^{J}$ associated with the treatment equation and set $\jmath=0$ as the base level. The assignment mechanism is described by

where $\bm{d}_i=(d_i(0),\cdots, d_i(J))'$, $\bm{p}_i=(p_{i,0},\cdots, p_{i,J})'$, $\bm{v}_i=(v_{i,0},\cdots, v_{i,J})'$, $\bm{\rho}(\cdot)=(\rho_1(\cdot),\cdots,\rho_J(\cdot))'$, and $\boldsymbol{\Gamma}=(\boldsymbol{\gamma}_1, \cdots, \boldsymbol{\gamma}_J)'$ for $\boldsymbol{\gamma}_\jmath\in\mathbb{R}^{\mathsf{d}_z}$, $\jmath=1, \ldots, J$. Notice that each $p_{i,\jmath}$ for $\jmath\in\mathcal{J}$ is the conditional probability of treatment level $\jmath$, which would be the usual propensity score if $\boldsymbol{\alpha}_i$ were observable.

An important feature of this model is that $\bm{z}_i$ and $\boldsymbol{\alpha}_i$ enter the two equations simultaneously, implying that they play the role of confounders in the potential outcomes framework. For simplicity, $\varsigma_{i,\jmath}$ and $p_{i,\jmath}$ are assumed to take generalized partially linear forms: the unobserved $\boldsymbol{\alpha}_i$ enters the model nonparametrically through the unknown functions $\mu_\jmath(\cdot)$ and $\rho_\jmath(\cdot)$, whereas the observed $\bm{z}_i$ is controlled for in an additive-separable way.

Introducing the link functions $\psi_{\mathsf{y}}^{-1}$ and $\boldsymbol{\psi}_{\mathsf{s}}^{-1}$ is convenient in practice, but it is less relevant to the core idea of this paper and notationally cumbersome. Thus, the discussion of this general case is deferred to Section (ref). For the moment, I make the first simplification of the general model by specifying identity links:

Structure of Large-Dimensional Measurements

The observed covariates $\bm{x}_i=(x_{i1},\cdots, x_{iT})'$ play the role of noisy measurements of latent confounders $\boldsymbol{\alpha}_i$. This paper considers a general covariates-adjusted nonlinear factor model for $\bm{x}_i$. Specifically, partition $\bm{w}_i\in\mathbb{R}^{T\mathsf{d}_w}$ into $T$-vectors of covariates: $\bm{w}_i=(\bm{w}_{i,1}',\cdots, \bm{w}_{i,\mathsf{d}_w}')'$ where $\bm{w}_{i,\ell}=(w_{i1,\ell}, \cdots, w_{iT,\ell})'$ for $\ell=1,\cdots, \mathsf{d}_w$. The measurements $\bm{x}_i$ are characterized by the following model:

where $\mathcal{F}$ is a $\sigma$-field generated by unobserved random elements $\{\boldsymbol{\alpha}_i\}_{i=1}^n$ and $\{\eta_t(\cdot)\}_{t=1}^T$.

Equation (ref) is indeed a linear regression model with an unknown possibly nonlinear factor component. The regressors $\{\bm{w}_{i,\ell}\}_{\ell=1}^{\mathsf{d}_w}$ need to be sufficiently high-rank (enough variation across both $i$ and $t$) for the identification of $\{\vartheta_\ell\}_{\ell=1}^{\mathsf{d}_w}$. Since incorporating $\{\bm{w}_{i,\ell}\}_{\ell=1}^{\mathsf{d}_w}$ is notationally cumbersome and less relevant to the core idea of this paper, the discussion is deferred to Section (ref). For the moment, I make the second simplification by setting $\vartheta_\ell=0$ for all $\ell=1,\cdots, \mathsf{d}_w$:

Let $\boldsymbol{\eta}(\cdot)=(\eta_{1}(\cdot), \cdots, \eta_{T}(\cdot))'$ and $\bm{u}_{i}=(u_{i1},\cdots,u_{iT})'$. Define $T\times n$ matrices $\bm{X}=(\bm{x}_1,\cdots, \bm{x}_n)$, $\boldsymbol{\eta}=(\boldsymbol{\eta}(\boldsymbol{\alpha}_1),\cdots, \boldsymbol{\eta}(\boldsymbol{\alpha}_n))$ and $\bm{u}=(\bm{u}_1,\cdots, \bm{u}_n)$. Equation (ref) can be written in matrix form: $\bm{X}=\boldsymbol{\eta}+\bm{u}$.

Throughout the paper, the latent variables $\{\boldsymbol{\alpha}_i\}_{i=1}^n$ and the latent functions $\{\eta_t(\cdot)\}_{t=1}^T$ are understood as random elements, but the main analysis is conducted conditional on them. In this sense, they are analogous to fixed effects in the panel data literature. The number of latent variables $\mathsf{d}_\alpha$ is assumed to be known in the theoretical analysis. In practice, however, it is often unknown and may need to be determined by the researcher using, for example, selection techniques developed in the factor analysis literature Bai-Ng_2002_ECMA,Ahn_2013_ECMA. See Remark (ref) for more discussion. A formal procedure for determining $\mathsf{d}_\alpha$ is left for future research.

This setup indeed encompasses many examples in the literature. Suppose that $\eta_t(\boldsymbol{\alpha}_i)=\alpha_i+\varpi_t$ for some $\varpi_t\in\mathbb{R}$. Then, Equation (ref) reduces to the classical two-way fixed effects model in panel data analysis. If, instead, we assume $\eta_t(\boldsymbol{\alpha}_i)=\bm{\varpi}_t'\boldsymbol{\alpha}_i$ for some $\bm{\varpi}_t\in\mathbb{R}^{\mathsf{d}_\alpha}$, Equation (ref) reduces to an interactive fixed effects model Bai_2009_ECMA. In fact, the two-way fixed effects, interactive fixed effects and many other popular methods in empirical studies implicitly restrict the latent mean structure $\boldsymbol{\eta}$ to be exactly low-rank. In contrast, this paper allows $\boldsymbol{\eta}$ to be full rank due to the potential nonlinearity of the latent functions $\{\eta_t(\cdot)\}_{t=1}^T$, while the variation of the large-dimensional $\bm{x}_{i}$ may still be explained by a few low-dimensional components in a possibly nonlinear way.

Notation

Latent functions. For a generic sequence of functions $\{h_t(\cdot)\}_{t=1}^M$ defined on a compact support, let $\nabla^{\ell}\bm{h}_t(\cdot)$ be a vector of $\ell$th-order partial derivatives of $h_t(\cdot)$, and define $\mathscr{D}^{[\kappa]}\bm{h}_t(\cdot)=(\nabla^{0}\bm{h}_t(\cdot)', \cdots, \nabla^{\kappa}\bm{h}_t(\cdot)')'$, i.e., a column vector that stores all partial derivatives of $h_t(\cdot)$ up to order $\kappa$. The derivatives on the boundary are understood as limits with the arguments ranging within the support. When $\ell=1$, $\nabla\bm{h}_t(\cdot):=\nabla^1\bm{h}_t(\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}\|_2=\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.

Asymptotics. For sequences of numbers or random variables, $a_n\lesssim b_n$ denotes $\limsup_n|a_n/b_n|$ is finite, and $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$. $\rightsquigarrow$ denotes convergence in distribution.

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{q}]=\sum_{j=1}^{\mathsf{d}}q_j$ and $\bm{v}^{\bm{q}}=v_1^{q_1}v_2^{q_2}\cdots v_{\mathsf{d}}^{q_{\mathsf{d}}}$.

Outline of Estimation Procedure

This section describes the main procedure for counterfactual analysis, which consists of three steps. First, relevant information on $\boldsymbol{\alpha}_i$ is extracted based on Equation (ref). Second, the conditional means $\{\varsigma_{i,\jmath}\}_{i=1}^n$ of potential outcomes and conditional treatment probabilities $\{p_{i,\jmath}\}_{i=1}^n$ are estimated by local least squares where the extracted information from the first step plays the role of kernel functions and generated regressors. Third, estimators of causal parameters of interest are constructed based on doubly-robust score functions. See Algorithm \hyperlink{t2}{1} for a short summary. The main tuning parameter in this procedure is the number of nearest neighbors $K$, which governs the bandwidth of nonparametric regression in the second step. The number of principal components to be extracted $\mathsf{d}_{\lambda}$ can be either fixed or selected the way described in Remark (ref).

In addition to methodological discussions, each step will be accompanied by an empirical illustration using the data of Acemoglu-et-al_2016_JFE, which analyzes the effect of the announcement of the appointment of Tim Geithner as Treasury Secretary on November 21, 2008 on stock returns of financial firms that were connected to him. This study can be viewed as an example of the synthetic control design in the program evaluation literature (see Abadie_2020_JEL for a review). Specifically, the treatment of interest is the appointment of Geithner, which starts at a particular date (referred to as “event day $0$" hereafter). All firms remain untreated prior to the appointment. Starting at event day $0$, a subgroup of firms that are connected to Geithner are treated ($s_i=1$), while the other group remains untreated ($s_i=0$). Variables used in this analysis and the parameter of interest are listed in the following.

The sample consists of $583$ firms in total ($n=583$) and $22$ of them are treated (“connected to Geithner"). To be comparable with the results in Acemoglu-et-al_2016_JFE, the observed measurements $\bm{x}_i$ only include stock returns for $250$ days that ends $30$ days prior to the Geithner announcement ($T=250$). Note that the proposed method is not restricted to synthetic control studies illustrated by this example. See Section SA-4 of the SA for other applications.

Step 1: Latent Variables Extraction

The goal is to extract information on latent confounders by employing Equation (ref). The main ideas are sketched below. Section SA-2 of the SA provides discussion of a more general local principal subspace approximation procedure.

Row-wise Splitting. Split the row index set $\mathcal{T}=\{1,\cdots, T\}$ of $\bm{X}$ into two non-overlapping subsets randomly: $\mathcal{T}=\mathcal{T}^\dagger\cup\mathcal{T}^\ddagger$ with $T^\dagger=|\mathcal{T}^\dagger|$, $T^\ddagger=|\mathcal{T}^\ddagger|$ and $T^\dagger\asymp T^\ddagger\asymp T$. Accordingly, the data matrix $\bm{X}$ is divided into two submatrices $\bm{X}^{\dagger}$ and $\bm{X}^{\ddagger}$ with row indices in $\mathcal{T}^\dagger$ and $\mathcal{T}^\ddagger$ respectively. $\bm{u}^{\dagger}$ and $\bm{u}^{\ddagger}$ are defined similarly. This step is needed only when local PCA is implemented.

$K$-Nearest Neighbors Matching. This step makes use of the subsample labeled by $\dagger$, i.e., the submatrix of $\bm{X}$ with row indices in $\mathcal{T}^{\dagger}$. For a generic unit $i\in\{1,\cdots, n\}$, search for a set of indices $\mathcal{N}_i$ for its $K$ nearest neighbors (including $i$ itself) in terms of a distance metric $\mathfrak{d}(\cdot, \cdot)$:

Usual choices include Euclidean distance $\mathfrak{d}_2(\bm{X}^{\dagger}_{\cdot i}, \bm{X}^{\dagger}_{\cdot j})=\frac{1}{\sqrt{T^\dagger}}\|\bm{X}^{\dagger}_{\cdot i}-\bm{X}^{\dagger}_{\cdot j}\|_2$ and pseudo-max distance $\mathfrak{d}_\infty(\bm{X}^{\dagger}_{\cdot i}, \bm{X}^{\dagger}_{\cdot j})=\max_{l\neq i, j} |\frac{1}{T^\dagger}(\bm{X}^{\dagger}_{\cdot i}-\bm{X}^{\dagger}_{\cdot j})'\bm{X}^{\dagger}_{\cdot l}|$. The latter, proposed by Zhang-Levina-Zhu_2017_BIMA, has appealing features. In particular, it may accommodate (conditional) heteroskedasticity of errors in the nonlinear factor model, and under Assumption (ref) below, matching on the noisy measurements using $\mathfrak{d}_\infty(\cdot,\cdot)$ translates into a sharp bound on the matching discrepancy of the underlying latent variables. From now on, attention is restricted to results based on $\mathfrak{d}(\cdot,\cdot)=\mathfrak{d}_\infty(\cdot,\cdot)$. Properties of Euclidean distance are discussed in Section SA-2 of the SA. Moreover, when the noisy measurements differ in scale or importance for revealing information on the latent variables, it may be desirable to rescale or reweight different measurements when searching for nearest neighbors. Such transformations can be viewed as particular choices of the distance metric. See Remark (ref) below for more discussion.

The number of nearest neighbors $K$ is the main tuning parameter of the entire estimation procedure. In practice, following the discussion below Theorem (ref), one may use, for example, cross validation or some plug-in rules, to select an optimal $K$ that minimizes the mean squared error of the estimators of $\{\varsigma_{i,\jmath}\}_{i=1}^n$ or $\{p_{i,\jmath}\}_{i=1}^n$ (see Step 2 in Section (ref)). Under mild conditions, this choice can be used to construct a valid inference procedure in the last step.

As a conceptual illustration, Figure (ref) shows an artificial two-dimensional surface embedded in a three-dimensional space. $K$-NN matching for a particular unit $i$ (colored red) generates a local neighborhood (the circled region). Note that the distance $\mathfrak{d}_\infty(\cdot,\cdot)$ is defined based on averaging information across the $t$ dimension. If the errors in $\bm{u}_i$ are independent or weakly dependent across $t$, their impact on the distance becomes negligible as the dimensionality $T$ grows large. On the other hand, if the (noise-free) latent factor structure is not too singular (see Assumption (ref) below), any two points found close on the surface should also be close in terms of the underlying latent variables. Therefore, the nearest neighbors obtained by matching on the observed measurements are similar in terms of the unobserved confounders, which is the key building block of subsequent analysis.

\FloatBarrier

\FloatBarrier

Using the data of Acemoglu-et-al_2016_JFE, I implement $K$-NN matching for each unit based on stock returns in the first 125 days with $K=127$. This relatively large $K$ is selected based on a data-driven procedure described in Step 2 in Section (ref). Due to the noise in the measurements, choosing a small $K$ may not help reduce the resultant matching discrepancy (see discussions below Theorem (ref)). To have a sense of the performance of $K$-NN matching, I calculate for each unit the maximum distance of matched pairs divided by the standard deviation of the distance across all pairs, which can be viewed as a normalized matching discrepancy in terms of the observed returns. Table (ref) reports some summary statistics for treated and control groups respectively. Matching performs well for treated units, whereas some control units are matched with someone relatively far away. In later analysis, I will check the robustness of the results by varying the number of nearest neighbors or dropping a few control units with large discrepancy.

\FloatBarrier

\FloatBarrier

Local Principal Component Analysis. This step makes use of the subsample labeled by $\ddagger$, i.e., the submatrix of $\bm{X}$ with row indices in $\mathcal{T}^{\ddagger}$. Given a set of nearest neighbors $\mathcal{N}_i$ from the previous step, define a $T^\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$. For these nearest neighbors, the unknown function $\eta_t(\cdot)$ can be locally approximated by a linear combination of some basis functions of latent variables. Then, $\bm{X}_{\langle i \rangle}$ admits a linear factor structure up to approximation errors:

where $\bm{u}_{\langle i \rangle}=(\bm{u}^{\ddagger}_{\cdot j_1(i)}, \cdots, \bm{u}^{\ddagger}_{\cdot j_K(i)})$. $\bm{F}_{\langle i \rangle}\boldsymbol{\Lambda}_{\langle i \rangle}'+\bm{r}_{\langle i \rangle}$ is the possibly nonlinear factor component. The $K\times \mathsf{d}_\lambda$ matrix $\boldsymbol{\Lambda}_{\langle i \rangle}$ can be viewed as approximation basis functions of latent confounders (evaluated at the data points), the $T^\ddagger\times \mathsf{d}_\lambda$ matrix $\bm{F}_{\langle i \rangle}$ collects the corresponding coefficients, and $\bm{r}_{\langle i \rangle}$ is the resultant approximation error. The user-specified parameter $\mathsf{d}_\lambda$ governs the number of approximation terms. $\bm{F}_{\langle i \rangle}$ and $\boldsymbol{\Lambda}_{\langle i \rangle}$, referred to as factor and loading matrices respectively, can be identified up to a rotation and estimated by PCA Bishop_2006_bookpattern:

such that $\frac{1}{T^\ddagger}\tilde\bm{F}_{\langle i \rangle}'\tilde\bm{F}_{\langle i \rangle}=\bm{I}_{\mathsf{d}_{\lambda}}$ and $\frac{1}{K}\tilde\boldsymbol{\Lambda}_{\langle i \rangle}'\tilde\boldsymbol{\Lambda}_{\langle i \rangle}$ is diagonal. Let $\widehat{\boldsymbol{\lambda}}_{\ell,\langle i \rangle}$ be the column in $\widehat{\boldsymbol{\Lambda}}_{\langle i \rangle}$ that corresponds to a generic unit $\ell$.

The idea underlying (ref), i.e., applying PCA locally to neighbors of each unit, 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 across the $t$ dimension, for units within the same local neighborhood, the nonlinear factor components $\bm{F}_{\langle i \rangle}\boldsymbol{\Lambda}_{\langle i \rangle}+\bm{r}_{\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 is independent (or weakly dependent) across $t$. Note that splitting is not necessary if a researcher believes that $K$-NN matching suffices for later analysis.

The decomposition (ref) is primarily of theoretical interest. In practice, there is no need to specify a particular approximation basis $\boldsymbol{\Lambda}_{\langle i \rangle}$ for implementing PCA as in (ref). Also, the number of local principal components to be extracted ($\mathsf{d}_{\lambda}$) plays a similar role as the degree of the basis in local polynomial regression and can be set as a fixed number (independent of $n$ and $T$). Two strategies may be employed:

Note that whichever strategy is used, we can at most achieve the approximation power such that the order of approximation bias does not exceed that of the noise. See more discussion about determining $\mathsf{d}_{\lambda}$ in Remark (ref).

The idea of local PCA is illustrated in Figure (ref). Units around the red dot are approximately lying on a (local) linear tangent plane. 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 factors can be differentiated from the noise, a local nonlinear principal subspace can be constructed for a higher-order approximation of the underlying surface.

\FloatBarrier

\FloatBarrier

Using the data of Acemoglu-et-al_2016_JFE, I implement local PCA for each unit. Recall that for each firm a set of nearest neighbors has been obtained using the stock returns in the first $125$ days. PCA can be conducted for this subgroup of firms using their stock returns in the next 125 days. Figure (ref) shows several leading eigenvalues corresponding to the neighborhood for a particular unit (“AMERICAN EXPRESS CO."), suggesting that extracting one or two local principal components is a reasonable choice. In the subsequent analysis, I set $\mathsf{d}_{\lambda}=1$. Results based on $\mathsf{d}_{\lambda}=2$ are similar and omitted to conserve space.

\FloatBarrier

\FloatBarrier

Step 2: Factor-Augmented Regression

For the outcome equation (ref), the predicted value $\widehat{\varsigma}_{i,\jmath}$ for unit $i$ is given by

It can be viewed as a local least squares regression with generated regressors $\widehat{\boldsymbol{\lambda}}_{\ell,\langle i \rangle}$. The treatment equation (ref) can be treated exactly the same way. By regressing $d_\ell(\jmath)$ on $\bm{z}_{\ell}$ and $\widehat{\boldsymbol{\lambda}}_{\ell,\langle i \rangle}$ for $\ell\in\mathcal{N}_i$, one can obtain the predicted value $\widehat{p}_{i,\jmath}$. Note that as discussed in Section (ref) below, other approaches such as nonparametric logit or probit regression can also be employed to estimate these propensity scores.

In practice, I suggest taking $\mathsf{d}_\lambda$ as a fixed number and choosing the tuning parameter $K$ accordingly. For instance, we can focus on the subgroup at the treatment level $\jmath$, and let $\mathsf{d}_\lambda=\mathsf{d}_{\alpha}+1$ to achieve the same approximation power of local linear estimation. Then, $K$ can be chosen possibly through two strategies:

For the purpose of illustration, I implement a local regression of stock returns at date $t$ on the leading factor loading extracted previously, for each $t=-20,\cdots,0, 1$, where $t=0$ denotes the day when the treatment starts. Figure (ref) shows the fitted values in black and the observed daily returns in grey for the $22$ treated firms, and the result for American Express Co. is displayed in Figure (ref). Recall that the fitted values are the estimates of conditional expectations of stock returns without treatment given the latent variables. Clearly, after day $0$, many sequences of stock returns increase sharply compared with the corresponding fitted values.

\FloatBarrier

\FloatBarrier

Step 3: Counterfactual Analysis

The final step is to estimate the counterfactual means of potential outcomes, which forms the basis of estimators for other causal parameters. Specifically, consider $\theta_{\jmath,\jmath'}:=\mathbb{E}[y_{i}(\jmath)|s_{i}=\jmath']$. Let $p_{\ell}=\mathbb{P}(s_{i}=\ell)$ for any $\ell\in\mathcal{J}$. Under unconfoundedness conditional on $\bm{z}_i$ and $\boldsymbol{\alpha}_i$ (see Assumption (ref) below), \[ \theta_{\jmath,\jmath'}=\mathbb{E}\left[\frac{d_{i}(\jmath')\varsigma_{i,\jmath}}{p_{\jmath'}}+\frac{p_{i,\jmath'}}{p_{\jmath'}}\frac{d_{i}(\jmath)(y_{i}-\varsigma_{i,\jmath})}{p_{i,\jmath}}\right]. \] An estimator of $\theta_{\jmath,\jmath'}$ is given by

where $\widehat{\varsigma}_{i,\jmath}$, $\widehat{p}_{i,\jmath}$ and $\widehat{p}_{i,\jmath'}$ are obtained in the second step, and $\widehat{p}_\ell=\frac{1}{n}\sum_{i=1}^{n}d_i(\ell)$ for $\ell=\jmath,\jmath'$. For the purpose of inference, a simple plug-in variance estimator for $\widehat{\theta}_{\jmath,\jmath'}$ is

One may expect $\sqrt{n}\widehat{\sigma}_{\jmath,\jmath'}^{-1}(\widehat{\theta}_{\jmath,\jmath'}-\theta_{\jmath,\jmath'})\rightsquigarrow \mathsf{N}(0,1)$, which will be shown in Theorem (ref) below. Confidence intervals and hypothesis testing procedures can be constructed accordingly.

Estimators of other parameters may be constructed in a similar way or based on $\{\widehat{\theta}_{\jmath,\jmath'}\}_{\jmath,\jmath'\in\mathcal{J}}$. For example, the average treatment effect on the treatment group $\jmath'=\ell$ compared to the baseline treatment status $\jmath=0$ can be estimated by $\widehat{\theta}_{\ell,\ell}-\widehat{\theta}_{0,\ell}$ where $\widehat{\theta}_{\ell,\ell}=\sum_{i=1}^nd_i(\ell)y_i/\sum_{i=1}^nd_i(\ell)$.

As an illustration, I estimate the average treatment effect of Geithner connections on cumulative returns from day 0 to day 1 (CAR[0,1]) for firms with connections. Since the number of treated units is relatively small, the propensity score is obtained by taking a simple local average within each local neighborhood. For the outcome equation, I run a local regression of stock returns of firms with no connections on the factor loading extracted previously ($\mathsf{d}_\lambda=1$). Different choices of $K$ are considered, which correspond to $K=Cn^{2/3}$ where $C=0.5, 1, 2$, and I also use the strategy outlined in Section SA-5.3 of the SA to obtain a DPI choice $\widehat{K}_{\mathtt{DPI}}$ based on an initial choice $K=2n^{2/3}$. Assuming there exists one latent variable, such choices coincide with the mean squared error (MSE) optimal rate of the underlying nonparametric estimators. The above procedure is applied to the full sample and a base sample. The latter, as defined in Acemoglu-et-al_2016_JFE, excludes firms whose returns are highly correlated with Citigroup.

Results are reported in the first two columns of Table (ref). I also include two results based on synthetic matching in Acemoglu-et-al_2016_JFE and one result based on a penalized synthetic control method in Abadie-Lhour_2020_wp. To make these results comparable, I report the 95% confidence intervals for hypothesis testing of the average treatment effect (on the treated) being equal to zero (numbers in brackets in Table (ref)), but note that the underlying assumptions and inference methodology of the other two papers are different from mine. The estimated average cumulative abnormal return for the connected firms using the proposed method ranges from $0.065$ to $0.105$ and significantly differs from zero at the $0.05$ level. I also check the robustness of the results by excluding firms in the control group with large normalized matching discrepancy (top $10\%$ in Table (ref)). Results are similar and omitted to conserve space.

The analysis so far has controlled for latent variables only. Three additional covariates are available in the dataset of Acemoglu-et-al_2016_JFE: firm size (log of total assets), profitability (return on equity), and leverage (total debt to total capital) as of 2008. They can be incorporated into the local least squares regression in Step 2 as additional regressors $\bm{z}_i$. Results are reported in the third and fourth columns of Table (ref). The estimated effect is slightly smaller than that without additional covariates, but still significant at the 0.05 level.

\FloatBarrier

\FloatBarrier

Main Results

Assumptions

I begin with unconfoundedness and overlap conditions commonly used in the causal inference literature. Note that the conditioning variables $\boldsymbol{\alpha}_i$ in this scenario are unobservable.

The next assumption imposes mild regularity conditions on the treatment effects model and the latent structure of $\bm{x}_i$.

Parts (a), (b), and (c) concern the regularities of the treatment effects model characterized by Equations (ref) and (ref). The conditional expectations of potential outcomes and propensity scores are sufficiently smooth functions, and other standard conditions are imposed on the conditioning variables and errors. Regarding the latent structure of $\bm{x}_i$ described in Equation (ref), part (d) ensures that all latent functions belong to a H\"{o}lder class of order $\bar{m}$, and part (e) are standard conditions on errors commonly used in factor analysis and graphon estimation. The constant $\bar{m}$ governs the smoothness of unknown functions, and $\nu$ controls the tails of error terms. They are assumed to be the same across Equations (ref)-(ref) only for simplicity.

Recall that the main task of the first step is to learn $\boldsymbol{\alpha}_i$ from $\bm{x}_i$. $\boldsymbol{\alpha}_i$ is not identifiable, but it is unnecessary to identify it in an exact sense since only predictions based on $\boldsymbol{\alpha}_i$ matter for the main analysis. Intuitively, the goal is to extract local geometric relations among latent $\boldsymbol{\alpha}_i$'s, which are reflected by index sets for nearest neighbors and local factor loadings described in Section (ref). The next three assumptions detail the restrictions on the latent nonlinear structure generated by $\{\eta_t(\cdot)\}_{t=1}^T$, ensuring that the relations learned from $\bm{x}_i$'s can be translated into that for $\boldsymbol{\alpha}_i$'s.

Assumption (ref) below can be intuitively understood as an “identification” condition for $\boldsymbol{\alpha}_i$. It implies that the difference in $\boldsymbol{\alpha}_i$ can be revealed by the observed measurements $\bm{x}_i$ as the dimensionality of $\bm{x}_i$ grows large, though exact identification of $\boldsymbol{\alpha}_i$ is impossible without further restrictions. Due to the row-wise sample splitting, I will write $\boldsymbol{\eta}^\dagger(\cdot):=(\eta_t(\cdot))_{t\in\mathcal{T}^\dagger}$, and $\boldsymbol{\eta}^{\ddagger}(\cdot):=(\eta_t(\cdot))_{t\in\mathcal{T}^\ddagger}$, which are $T^\dagger\times 1$ and $T^\ddagger\times 1$ column vectors of latent functions respectively.

When this condition holds, units within the same neighborhood obtained by matching on the large-dimensional $\bm{x}_i$ (see Step 1 in Section (ref)) are similar in terms of the latent $\boldsymbol{\alpha}_i$. Otherwise, units with quite different values of $\boldsymbol{\alpha}_i$ could be matched.

The idea underlying this assumption is related to the completeness condition widely used in econometric identification problems Newey-Powell_2003_ECMA,Chernozhukov-Hansen_2005_ECMA,Hu-Schennach_2008_ECMA. 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, condition (ref) amounts to saying that there is enough variation observed on the latent surface for different values of the latent variables.

The next assumption concerns the non-collinearity of derivatives of latent functions, which facilitates the quantification of indirect matching and local PCA.

Assumption (ref) ensures that the derivatives of latent functions up to order $m-1$ are not too collinear. It is analogous to the non-degenerate factors condition commonly imposed in linear factor analysis Bai_2003_ECMA, Bai_2009_ECMA.

The last assumption on the latent structure, which I refer to as non-collapsing, permits accurate translation from the matching discrepancy of observables to that of unobservables. Specifically, let $\mathscr{P}_{\boldsymbol{\alpha}_0}[\cdot]$ be the projection operator onto the $\mathsf{d}_\alpha$-dimensional space (embedded in $\mathbb{R}^{T^\dagger}$) spanned by the local tangent basis at $\boldsymbol{\eta}^\dagger(\boldsymbol{\alpha}_0)$, i.e., $\nabla\boldsymbol{\eta}^\dagger(\boldsymbol{\alpha}_0)$. Take an orthogonalized basis of this tangent space. Denote by $\mathscr{P}_{\boldsymbol{\alpha}_0,\ell}[\cdot]$ the projection operator onto the $\ell$th direction of the tangent space. Then, $\mathscr{P}_{\boldsymbol{\alpha}_0,\ell}[\boldsymbol{\eta}^\dagger(\boldsymbol{\alpha})]$ is the projection of $\boldsymbol{\eta}^\dagger(\boldsymbol{\alpha})$ onto the $\ell$th direction of the tangent space at $\boldsymbol{\eta}^\dagger(\boldsymbol{\alpha}_0)$.

Though seemingly involved at first glance, condition (ref) has an intuitive geometric interpretation. Note that the latent functions $\{\eta_t(\cdot)\}_{t\in\mathcal{T}^\dagger}$ generate a $\mathsf{d}_\alpha$-dimensional surface embedded in $\mathbb{R}^{T^\dagger}$, and thus (ref) simply says that if the whole surface is projected onto the tangent space at any data point, the dimensionality of the projection does not drop, as implied by the name “non-collapsing”.

Assumptions (ref)-(ref) are a group of lower-level conditions. The cornerstone of the proposed method is a local principal subspace approximation procedure, which has broad applicability and can be justified under higher-level conditions. In particular, the analysis below can be easily adapted to cover semi-strong factor models Wang-Fan_2017_AoS. See Section SA-2.2 of the SA for details.

Theoretical Results

Throughout the analysis below, I write $\delta_{KT}=(K^{1/2}\wedge T^{1/2})/\sqrt{\log (n\vee T)}$ and $h_{K,\alpha}=(K/n)^{1/\mathsf{d}_\alpha}$. Recall that $K$ is the number of nearest neighbors and $T$ is the dimensionality of $\bm{x}_i$. The asymptotic analysis is conducted assuming both $K$ and $T$ diverge as $n\rightarrow\infty$. As explained before, I use $\mathfrak{d}(\cdot,\cdot)=\mathfrak{d}_\infty(\cdot,\cdot)$. Moreover, $\mathsf{d}_\lambda$ is the number of extracted leading local principal components. I assume that $\mathsf{d}_\lambda=\binom{m-1+\mathsf{d}_\alpha}{\mathsf{d}_\alpha}$ so that the local approximation terms up to order $(m-1)$ for $\{\eta_t(\cdot)\}$ are extracted.

Latent Structure

I begin with the covariate equation (ref) and provide some important intermediate results that may be of independent interest. More detailed results are given in Section SA-2 of the SA. The first theorem concerns the discrepancy of latent variables induced by matching on the observed measurements.

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)^{1/\mathsf{d}_\alpha}$ 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 Gyorfi-et-al_2002_bookchapter. The second term $\sqrt{\log n/T}$ arises from the existence of $\bm{u}_i$. By construction of the distance metric, averaging across the $t$ dimension can shrink the impact of $\bm{u}_i$ to the order of $T^{-1/2}$ up to a log penalty.

Note that if $\boldsymbol{\alpha}_i$ were observed, matching could be directly implemented on it with the number of matches $K$ fixed, and large-sample properties of the resultant matching estimators have been established in Abadie-Imbens_2006_ECMA. 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 in Theorem (ref). Using a fixed (or small) number of nearest neighbors is unable to further reduce bias and thus is not recommended in this scenario.

The rate restriction in Theorem (ref) is exploited in application of maximal inequalities, which ensures uniform convergence of sample averages across $t$. It becomes more relaxed when more stringent moment conditions (larger $\nu$) hold. In addition, as explained earlier, Assumption (ref) is used to derive a sharp bound on the indirect matching discrepancy when $\mathfrak{d}_\infty(\cdot, \cdot)$ is used. When it does not hold, a loose bound may still be established. See Theorem SA-2.2 and Remark SA-2.1 in the SA for details.

Next, I consider the properties of the estimated factors and loadings from local PCA. Note that the decomposition of $\bm{X}_{\langle i \rangle}$ given by Equation (ref) is still arbitrary since the approximation basis $\boldsymbol{\Lambda}_{\langle i \rangle}$ is undefined. From a practical perspective, users do not need specify $\boldsymbol{\Lambda}_{\langle i \rangle}$ explicitly, and the PCA procedure has automatically imposed normalization so that the estimated factors and loadings are uniquely defined. For the purpose of theoretical analysis, however, $\boldsymbol{\Lambda}_{\langle i \rangle}$ needs to be appropriately defined so that it aligns with the estimand of $\widehat{\boldsymbol{\Lambda}}_{\langle i \rangle}$ and possesses approximation power for general smooth functions. Formally, let \[ \boldsymbol{\alpha}\in\mathcal{A}\mapsto \boldsymbol{\lambda}(\boldsymbol{\alpha}):=(\lambda_1(\boldsymbol{\alpha}), \cdots, \lambda_{\mathsf{d}_{\lambda}}(\boldsymbol{\alpha}))' \] be a $\mathsf{d}_\alpha$-variate monomial basis of degree no greater than $m-1$ centered at $\boldsymbol{\alpha}_i$ (including the constant). A typical element of $\boldsymbol{\lambda}(\boldsymbol{\alpha})$ is then given by $(\boldsymbol{\alpha}-\boldsymbol{\alpha}_i)^{\bm{q}}$ for $[\bm{q}]\leq m-1$. Define $\boldsymbol{\Lambda}_{\langle i \rangle}=(\boldsymbol{\lambda}(\boldsymbol{\alpha}_{j_1(i)}), \cdots, \boldsymbol{\lambda}(\boldsymbol{\alpha}_{j_K(i)}))'$. Heuristically, $K$-NN matching has detected a group of units with similar latent features, and $\boldsymbol{\Lambda}_{\langle i \rangle}$ further characterizes their local relations that may be used for higher-order approximation,

Note that by Theorem (ref), the distance between $\boldsymbol{\alpha}_i$ and its nearest neighbors $\{\boldsymbol{\alpha}_{j_k(i)}\}_{k=1}^K$ is diminishing as $n$ diverges. Consequently, the loadings of different factors in Equation (ref) are possibly shrinking in magnitude at heterogeneous rates. The next theorem, as the key building block of the main results, shows that $\boldsymbol{\Lambda}_{\langle i \rangle}$ can be estimated up to a rotation provided that the leading approximation terms included in $\boldsymbol{\Lambda}_{\langle i \rangle}$ are sufficiently strong.

The estimation errors of $\widehat{\boldsymbol{\Lambda}}_{\langle i \rangle}$ consists of two parts. The first term reflects the estimation variance. Since latent variables are not observed, this rate of convergence relies on both $K$ and $T$, as in linear factor analysis Bai_2003_ECMA. The second term is simply the resulting approximation error. The rate condition $\delta_{KT}h_{K,\alpha}^{m-1}\rightarrow\infty$ ensures that the leading terms $\bm{F}_{\langle i \rangle}\boldsymbol{\Lambda}_{\langle i \rangle}'$ can be differentiated from the remainder in Equation (ref), and the other two are used in application of maximal inequalities. If $T\asymp n$, the first rate condition reduces to $(n/K)^{\frac{2m-2}{\mathsf{d}_{\alpha}}}=o(K/\log n)$, and the second and third ones can be combined and simplified to $n^{\frac{4}{\nu}}\lesssim K/\log n$. In this simple case, for $K=n^{A}$, $A>\max\{\frac{4}{\nu}, \frac{2m-2}{2m-2+\mathsf{d}_{\alpha}}\}$ suffices. In particular, if $\nu$ is sufficiently large, this restriction can be satisfied by setting, for example, $K\asymp n^{\frac{2m}{2m+\mathsf{d}_{\alpha}}}$ (or equivalently, $h_{K,\alpha}\asymp n^{-\frac{1}{2m+\mathsf{d}_{\alpha}}}$), which coincides with the MSE-optimal choices of tuning parameters in the nonparametrics literature.

The result of Theorem (ref) indeed concerns the convergence of $\widehat{\boldsymbol{\Lambda}}_{\langle i \rangle}$ in terms of sup-norm, which is also uniform over the local neighborhoods indexed by $\langle i \rangle$. The key proof strategy is a leave-one-out trick used in studies of principal components, e.g., Abbe-et-al_2020_AoS. It helps construct sup-norm bounds on estimated singular values. Similar results can also be established for the estimated factors. The rate of convergence may be heterogeneous across columns of $\widehat{\bm{F}}_{\langle i \rangle}$, reflecting the differing magnitude of underlying approximation terms. See Theorem SA-2.4 of the SA for details. Note that the uniform convergence in Theorem (ref) is convenient for later analysis, but the rate conditions required may be stronger than needed for pointwise or $L_2$ convergence.

Before moving to the next step, I show the uniform convergence of the estimated common components, which may be of independent interest for panel data analysis. Let $\boldsymbol{\eta}^{\ddagger}$ be the submatrix of $\boldsymbol{\eta}$ with row indices in $\mathcal{T}^{\ddagger}$. Write $\boldsymbol{\eta}_{\langle i \rangle}=(\boldsymbol{\eta}^{\ddagger}_{\cdot j_1(i)}, \cdots, \boldsymbol{\eta}^{\ddagger}_{\cdot j_K(i)})$ and $\widehat{\boldsymbol{\eta}}_{\langle i \rangle}=\widehat{\bm{F}}_{\langle i \rangle}\widehat{\boldsymbol{\Lambda}}_{\langle i \rangle}'$.

This theorem shows that the nonlinear factor components can be consistently estimated, and the convergence is uniform over both dimensions. It plays an important role in latent variables extraction when additional high-rank regressors are used as in Equation (ref).

Counterfactual Analysis

I first show the uniform convergence of the estimated conditional means of potential outcomes and propensity scores obtained through local factor-augmented regressions.

The convergence above should be read as uniform over all the data points indexed by $i$, which respects the fact that $\boldsymbol{\alpha}_i$ is not directly observed and we obtain information on it for the $n$ units in the dataset. In this sense, it slightly differs from some semiparametric analysis where uniformity over the whole support is established (or assumed directly). Again, the estimation errors reflect both variance and bias, including the impact of the generated regressors $\widehat{\boldsymbol{\Lambda}}_{\langle i \rangle}$.

Now, I am ready to apply previous results to inference on the counterfactual means of potential outcomes. The following theorem establishes the asymptotic normality of the proposed estimator.

As discussed before, the doubly-robust score function helps relax the condition on the convergence rates of $\widehat{\varsigma}_{i,\jmath}$ and $\widehat{p}_{i,\jmath}$. In line with the results in the double/debiased machine learning literature, the fourth rate condition essentially requires that the product of two estimation errors be of smaller order than $n^{-1/2}$, which can be satisfied, for example, when the convergence in Theorem (ref) is faster than $n^{-1/4}$. As discussed below Theorem (ref), when $T\asymp n$ and $\nu$ is sufficiently large, the first three restrictions can be satisfied by $K\asymp n^{\frac{2m}{2m+\mathsf{d}_{\alpha}}}$. The resultant convergence rates of $\widehat{\varsigma}_{i,\jmath}$ and $\widehat{p}_{i,\jmath}$ coincide with the usual MSE-optimal rates (up to a log term) in the nonparametrics literature, which suffices to satisfy the faster-than-$n^{1/4}$ requirement if $m>\mathsf{d}_{\alpha}/2$.

Note that this paper focuses on large-$K$ asymptotics, which is analogous to a large bandwidth in kernel estimation or a small number of approximation terms in series estimation. If $K$ is small relative to the sample size, a non-negligible undersmoothing bias may arise in the distributional approximation, and bias-robust inference may be needed. See Cattaneo-Jansson_2018_ECMA,Cattaneo-Jansson-Ma_2019_REStud,Matsushita-Otsu_2019_wp for more discussions of undersmoothing bias and possible solutions.

Numerical Results

I conducted a Monte Carlo investigation of the finite sample performance of the proposed method. I consider a binary treatment design $\mathcal{J}=\{0,1\}$. The potential outcomes are $y_{i}(0)=\alpha+\alpha^2+\epsilon_{i,0}$ and $y_{i}(1)=2\alpha+\alpha^2+1+\epsilon_{i,1}$. The treatment is $s_i=\mathds{1}(v_i\leq p_i)$ where $p_i=\exp((\alpha-0.5)+(\alpha-0.5)^2)/(1+\exp((\alpha-0.5)+(\alpha-0.5)^2))$. The observed covariates are generated based on $x_{it}=\eta_t(\alpha_i)+u_{it}$ with $\eta_t(\alpha_i)=(\alpha_i-\varpi_t)^2$ in Model 1 and $\eta_t(\alpha_i)=\sin(\pi(\alpha_i+\varpi_t))$ in Model 2. $\epsilon_{i,0}, \epsilon_{i,1}\sim \mathsf{N}(0,1)$ and $\alpha_i,\varpi_t, v_i\sim\mathsf{U}(0,1)$. $u_{it}\sim\mathsf{N}(0,1)$ and is i.i.d over $i$ and $t$. $\{\epsilon_{i,0}\}$, $\{\epsilon_{i,1}\}$, $\{\alpha_i\}$, $\{\varpi_t\}$, $\{v_i\}$ and $\{u_{it}\}$ are independent.

I consider 5,000 simulated datasets with $n=T=1,000$ each. For each simulated dataset, a point estimate of the counterfactual mean $\theta_{0,1}=\mathbb{E}[y_i(0)|s_i=1]$ is obtained. I report bias (BIAS), standard deviation (SD), root mean squared error (RMSE), coverage rate (CR) of nominal 95% confidence interval and its average length (AL) in Table (ref). The results in the first three rows (“local linear") are based on local PCA with two extracted principal components ($\mathsf{d}_\lambda=2$) combined with a two-fold row-wise sample splitting. For simplicity, the first half ($T^\dagger=500$) is used for nearest neighbors matching, and the second half ($T^\ddagger=500$) is used for local PCA. The number of nearest neighbors is taken to be $K=Cn^{4/5}$ for $C=0.5$, $1$, $1.5$ respectively. This rate coincides with the MSE-optimal choice in the (cross-sectional) nonparametric regression. Results reported in Row 4-6 (“local constant”) are based on the simple local average estimator described in Section (ref) without row-wise sample splitting. The number of nearest neighbors is taken to be $K=Cn^{2/3}$ for $C=0.5$, $1$, $1.5$ respectively. Using the strategy described in Section (ref), I also obtain the DPI choice of $K$, i.e., $\widehat{K}_{\mathtt{DPI}}$, based on an initial choice $K=1.5 n^{4/5}$ for local linear estimation and $K=1.5 n^{2/3}$ for local constant estimation. It turns out that the results are robust to the choice of $K$, though local constant approximation may have larger bias in some cases.

\FloatBarrier

\FloatBarrier

Extensions

Some useful extensions are discussed in this section. The first subsection extends the previous results to uniform inference on counterfactual distributions. The second concerns including additional regressors into the nonlinear factor model. The third extends the linear factor-augmented regression to generalized partially linear models.

Uniform Inference

In many applications, the outcome of interest is a certain transformation of the original potential outcome via a function $g(\cdot)\in\mathcal{G}$, and uniform inference over the function class $\mathcal{G}$ is desired. In general, the goal can be achieved in two steps: (i) strengthen the asymptotic expansion in Theorem (ref)(a) to be uniform, that is, the remainder needs to be negligible uniformly over $g\in\mathcal{G}$; (ii) show that the influence function as a process indexed by $\mathcal{G}$ weakly converges to a limiting process. The general treatment of such issues can be found in, e.g., Barrett-Donald_2003_ECMA, Chernozhukov_2013_ECMA,Donald-Hsu_2014_JoE.

I will focus on counterfactual distributions, the analysis of which relies on a particular function class $\mathcal{G}=\{y\mapsto \mathds{1}(y\leq \tau):\tau\in\mathcal{Y}\}$. Each $g(\cdot)\in\mathcal{G}$ corresponds to a particular value $\tau\in\mathcal{Y}$. Therefore, I will write $y_{i,\tau}(\jmath)=\mathds{1}(y_i(\jmath)\leq \tau)$ and $y_{i,\tau}=\mathds{1}(y_i\leq \tau)$. Accordingly, Equation (ref) becomes $$y_{i,\tau}(\jmath)=\varsigma_{i,\jmath,\tau}+\epsilon_{i,\jmath,\tau},\quad\varsigma_{i,\jmath,\tau}=\bm{z}_i'\boldsymbol{\beta}_{\jmath,\tau}+\mu_{\jmath,\tau}(\boldsymbol{\alpha}_i),$$ where $\varsigma_{i,\jmath,\tau}=\mathbb{P}(y_i(\jmath)\leq \tau|\bm{z}_i,\boldsymbol{\alpha}_i)$. For each $\tau$, the second step of the estimation procedure in Section (ref) is implemented to obtain an estimator $\widehat{\varsigma}_{i,\jmath,\tau}$ of $\varsigma_{i,\jmath,\tau}$. The parameter of interest is $\theta_{\jmath,\jmath'}(\tau)=\mathbb{E}[\mathds{1}(y_i(\jmath)\leq \tau)|s_i=\jmath']$, the counterfactual distribution function of $y_i(\jmath)$ for the group with\ treatment status $\jmath'$. From the perspective of uniform inference, $\theta_{\jmath,\jmath'}(\cdot)$ is a parameter in $\ell^\infty(\mathcal{Y})$, a function space of bounded functions on $\mathcal{Y}$ equipped with sup-norm. To establish the limiting distribution of the proposed estimator, I slightly strengthen the smoothness condition used in Assumption (ref).

The following theorem shows that the (rescaled) counterfactual distribution process weakly converges to a limiting Gaussian process indexed by $\tau\in\mathcal{Y}$, which forms the basis of uniform inference. See vandeVaart_book_1996 for underlying technical details.

Under proper regularity conditions, the weak convergence above can be applied to construct inference procedures for other quantities such as quantile treatment effects by the functional delta method. See Section SA-6.2 of the SA for details.

The limiting Gaussian process can be approximated based on a practically feasible multiplier bootstrap procedure widely used in the literature. To be specific, take an i.i.d sequence of random variables $\{\omega_i\}_{i=1}^n$ independent of the data with mean zero and variance one. Define a uniformly consistent estimator of $\varphi_{i,\jmath,\jmath'}(\cdot)$: $$\widehat{\varphi}_{i,\jmath,\jmath'}(\cdot)= \frac{d_i(\jmath')(\widehat{\varsigma}_{i,\jmath,\cdot}-\widehat{\theta}_{\jmath,\jmath'}(\cdot))}{\widehat{p}_{\jmath'}}+ \frac{\widehat{p}_{i,\jmath'}}{\widehat{p}_{\jmath'}}\frac{d_i(\jmath)(y_{i,\cdot}-\widehat{\varsigma}_{i,\jmath,\cdot})}{\widehat{p}_{i,\jmath}}.$$ The following corollary shows that conditional on the data, $\frac{1}{\sqrt{n}}\sum_{i=1}^n\omega_i\widehat{\varphi}_{i,\jmath,\jmath'}(\cdot)$ weakly converges to the same limiting process $\mathsf{Z}_{\jmath,\jmath'}(\cdot)$ as in Theorem (ref). In practice, one only needs to simulate this feasible approximation process by taking random draws of $\{\omega_i\}_{i=1}^n$.

To showcase the uniform inference procedure, I use the data of Acemoglu-et-al_2016_JFE to check the (first-order) stochastic dominance (SD) of $\theta_{1,1}(\cdot)$ over $\theta_{0,1}(\cdot)$, where $\theta_{1,1}(\cdot)$ and $\theta_{0,1}(\cdot)$ respectively denote the cumulative distribution functions (CDFs) of potential stock returns of firms connected to Geithner if they were connected and not connected with him. The main ideas are outlined here. By definition of SD, the null hypothesis is $\theta_{1,1}(\tau)\leq \theta_{0,1}(\tau)$ for all $\tau\in\mathcal{Y}$. An intuitive test statistic is $\sqrt{n}\sup_{\tau\in\mathcal{Y}}(\theta_{1,1}(\tau)-\theta_{0,1}(\tau))$. The null hypothesis is rejected if the test statistic is greater than a certain critical value. Given the asymptotic expansions of $\widehat{\theta}_{1,1}(\cdot)$ and $\widehat{\theta}_{0,1}(\cdot)$, the critical value can be obtained by simulating the supremum of the approximation process, i.e., $\sup_{\tau\in\mathcal{Y}}(\frac{1}{\sqrt{n}}\sum_{i=1}^n(\widehat{\varphi}_{i,1,1}(\tau)-\widehat{\varphi}_{i,0,1}(\tau)))$. In practice, the supremum over the whole support is simply replaced by maximum over a set of user-specified evaluation points.

For each $\tau\in\mathcal{Y}$, implement the estimation procedure described in Section (ref). Varying the values of $\tau$, I obtain two estimated distribution functions for firms with connections, as shown in Figure (ref). The treated outcome $Y(1)$ is the potential cumulative return with connections to Geithner and the untreated outcome $Y(0)$ refers to that without connections. To better understand the estimation uncertainty, $95\%$ confidence bands for the two estimated CDFs are plotted, which are based on simulating the maximum absolute value of the corresponding (studentized) approximation processes. In each case, the value of $\tau$ is restricted to range from $0.1$-quantile to $0.9$-quantile of the estimated distribution.

\FloatBarrier

\FloatBarrier

It turns out that the estimated CDF for the treated outcome is well below that for the untreated outcome. Formally, I take all observed values of cumulative returns as the evaluation points, and then simulate the maximum of the approximation process by taking $500$ draws of random weights $\{\omega_i\}_{i=1}^n$. In this simple example, the test statistic equals $0$, which is well below the critical value $7.5$ for a confidence level of $0.95$ obtained through simulation. Thus, SD of $\theta_{1,1}(\cdot)$ over $\theta_{0,1}(\cdot)$ cannot be rejected. It implies that the positive effects of political connections in this example are felt over the entire distribution of the stock returns of financial firms connected with Geithner, which is a stronger conclusion than that based simply on the mean in Section (ref).

High-Rank Covariates

The analysis so far is based on Equation (ref), assuming $\bm{x}_i$ takes a purely nonlinear factor structure. However, high-rank components may exist in $\bm{x}_i$, and it is the latent structure of the residuals that contains relevant information on $\boldsymbol{\alpha}_i$, as described by Equation (ref). It can be viewed as a generalization of linear regression with interactive fixed effects. Intuitively, due to the existence of the unknown $\eta_t(\boldsymbol{\alpha}_i)$, the regressors $\bm{w}_{i,\ell}$ have to be sufficiently high-rank, otherwise they will be too collinear with the latent component and $\{\vartheta_\ell\}_{\ell=1}^{\mathsf{d}_w}$ cannot be identified. This is similar to the identification condition for semiparametric partially linear regression.

The main analysis of this paper can still be applied once we have some consistent estimators of $\vartheta_\ell$'s. They can be obtained using the idea of partially linear regression. Specifically, I refer to Step 1 in Section (ref) as a general local principal subspace approximation procedure, which will be applied to other sequences in addition to $\{\bm{x}_i\}$. A slightly revised estimation procedure can be used to extract the latent variables:

Under additional regularity conditions on $\{\bm{w}_{i,\ell}\}_{\ell=1}^{\mathsf{d}_w}$, it can be shown that $\widehat{\boldsymbol{\vartheta}}$ converges to $\boldsymbol{\vartheta}$ sufficiently fast and the main results established previously still hold. Formal analysis is available in Section SA-6 of the SA and is omitted here to conserve space.

Generalized Partially Linear Forms

The analysis so far focuses on the simplified models (ref) and (ref). It can be extended to generalized partially linear models described by (ref) and (ref). We can employ the local quasi-maximum likelihood method. For Equation (ref), consider a quasi-log-likelihood function $\mathcal{L}_{\mathsf{y}}(\varsigma, y)$ such that $\frac{\partial}{\partial \varsigma}\mathcal{L}_{\mathsf{y}}(\varsigma,y)=\frac{y-\varsigma}{V_{\mathsf{y}}(\varsigma)}$ for some positive function $V_{\mathsf{y}}$. An estimator of $\varsigma_{i,\jmath}$ is given by \[

\] For each $i$, the fitting is restricted to its local neighborhood $\mathcal{N}_i$. Equation (ref) can be treated similarly by appropriately choosing a quasi-likelihood $\mathcal{L}_{\mathsf{s}}(\cdot,\cdot)$ associated with $\{V_{\mathsf{s},\jmath}(\cdot)\}_{\jmath\in\mathcal{J}}$ satisfying $\frac{\partial}{\partial\zeta_{\jmath}}\mathcal{L}_{\mathsf{s}}(\boldsymbol{\psi}_{\mathsf{s}}(\bm{\zeta}),\bm{d})=\frac{d(\jmath)-\psi_{\mathsf{s},\jmath}(\bm{\zeta})}{V_{\mathsf{s},\jmath}(\bm{\zeta})}$. The predicted conditional treatment probability is given by $\widehat{\bm{p}}_{i}=(\widehat{p}_{i,0},\cdots, \widehat{p}_{i,J})'=\boldsymbol{\psi}_{\mathsf{s}}(\widehat{\boldsymbol{\gamma}}\bm{z}_i+\bm{\widehat{\rho}}(\boldsymbol{\alpha}_i))$. The asymptotic properties of these estimators can be derived under additional regularity conditions on the quasi-likelihood and link functions. See Section SA-6 of the SA for details.

Alternatively, one may exploit other standard methods in the semiparametrics literature, e.g., profiled quasi-maximum likelihood, to estimate the parametric components $\boldsymbol{\beta}_\jmath$'s and $\boldsymbol{\gamma}_\jmath$'s, though it is computationally more burdensome. See Hardle-Muller-Sperlich-Werwatz_2012_book for implementation details.

Conclusion

This paper has developed a causal inference method for treatment effects models with some confounders not directly observed. Relevant information on these latent confounders is extracted from a large set of noisy measurements that admits an unknown, possibly nonlinear factor structure. Such information is then used to match comparable units in the subsequent counterfactual analysis. Large-sample properties of the proposed estimators are established. The results cover a large class of causal parameters, including average treatment effects and counterfactual distributions. The method is illustrated with an empirical application studying the effect of political connections on stock returns of financial firms.