# Estimating Counterfactual Matrix Means
with Short Panel Data††thanks: First version: December 11, 2023. Authors
listed alphabetically. We thank Dmitry Arkhangelsky, Susan Athey, Stéphane
Bonhomme, Kirill Borusyak, Jiafeng Chen, Rebecca Diamond, Matthew Gentzkow,
Bryan Graham, Christian Hansen, Guido Imbens, David Ritzwoller, Vasilis
Syrgkanis, and the participants of the Stanford Econometrics Lunch and the
Stanford Causal Panel Data Conference for their thoughtful comments and
valuable feedback.

Lihua Lei Stanford Graduate School of Business and Department of Statistics
(by courtesy); Email<EMAIL_ADDRESS>Brad Ross Stanford Graduate School
of Business; Email<EMAIL_ADDRESS>
###### Abstract

We develop a more flexible approach for identifying and estimating average
counterfactual outcomes when several but not all possible outcomes are
observed for each unit in a large cross section. Such settings include event
studies and studies of outcomes of “matches” between agents of two types, e.g.
workers and firms or people and places. When outcomes are generated by a
factor model that allows for low-dimensional unobserved confounders, our
method yields consistent, asymptotically normal estimates of counterfactual
outcome means under asymptotics that fix the number of outcomes as the cross
section grows and general outcome missingness patterns, including those not
accommodated by existing methods. Our method is also computationally
efficient, requiring only a single eigendecomposition of a particular
aggregation of any factor estimates constructed using subsets of units with
the same observed outcomes. In a semi-synthetic simulation study based on
matched employer-employee data, our method performs favorably compared to a
Two-Way-Fixed-Effects-model-based estimator.

Keywords: _panel data, missing not-at-random, low-rank, factor model,
interactive fixed effects, event study, bipartite network data_

##  1 Introduction

Researchers are often interested in estimating average counterfactual outcomes
in a population when “short” panel data with missing outcomes are available,
namely, observations of several but not all outcomes for each unit in a large
cross section. For example, in event study settings, units receive a treatment
at different times, units’ outcomes are observed over several time periods pre
and post-treatment, and researchers are interested in estimating average
counterfactual outcomes of treated units had they not been treated
(Ashenfelter Card, 1985; Bertrand ., 2004; Angrist Pischke, 2009). As another
example, several empirical literatures estimate average outcomes of
counterfactual “matches” between pairs of agents of two different “types,”
e.g. wages of workers at various firms (Abowd ., 1999), test scores of
students taught by various teachers (Jackson ., 2014), and earnings and health
outcomes of people living in various places (Finkelstein ., 2016; Chetty
Hendren, 2018; Card ., 2023).

In such settings, exogenous variation in which outcomes are observed for
different units is often hard to come by, and outcome missingness is plausibly
determined by unobserved, confounding variables that also affect the outcomes
themselves. To make progress, since only a few outcomes are observed per unit
(Andrews ., 2008; Bonhomme ., 2023), researchers often assume unobserved
confounders are low-dimensional and use multiple observations per unit to
estimate an explicit model of how those confounders affect
outcomes.111Usually, these models also require “strict exogeneity,” namely
that, conditional on low-dimensional confounders, outcomes are independent of
missingness. In keeping with much of applied practice, this paper does the
same. However, Ashenfelter Card (1985) and Bonhomme . (2019) discuss strict
exogeneity’s plausibility in event study and match outcome contexts,
respectively.

A canonical model applied in these settings is the Two-Way Fixed Effects
(TWFE) model, which assumes outcomes only differ systematically within and
across units in levels.222In event study settings, a TWFE model of control
potential outcomes across time periods implies that units’ control potential
outcomes follow “parallel trends” over time (Borusyak ., 2023; Ghanem .,
2022). While the TWFE model enables outcome means to be identified and
estimated under a myriad of outcome missingness patterns and a small number of
observations per unit (Jochmans Weidner, 2019; Borusyak ., 2023), it
meaningfully restricts how unobserved confounders can affect outcomes (Ghanem
., 2022).

In event study settings, a large literature has sought to allow for richer
confounding relationships than those implied by parallel-trends-type
assumptions by using a low-rank factor model of outcomes (see Section 1.1 for
references).333Factor models are also frequently called “interactive fixed
effect” models. In contrast with the TWFE model, the factor model allows
potentially multidimensional unobserved confounders to have differential
effects on different outcomes. However, existing factor model-based methods
cannot be applied generally, both because they only work under certain outcome
missingness patterns, and because many require a large number of observed
outcomes per unit.

In this paper, we seek to bridge the gap between the versatility of TWFE-based
methods and the expressivity of factor model-based methods. In particular, we
develop an approach for identifying, estimating, and conducting valid
inference on counterfactual outcome means under factor models in short panels
with general outcome missingness patterns.

(a) Event Study With Staggered Entry

(b) Bipartite Match Outcomes

Figure 1.1: Each of these figures illustrates which outcomes, corresponding to
columns, are observed for which units, corresponding to units, in two real-
world panel datasets. A given coordinate is colored dark red if the outcome
corresponding to the coordinate’s column is observed for the unit
corresponding to the coordinate’s row and is colored light red otherwise.
Figure 1(a) is generated using data from the longitudinal congestion pricing
study studied in Ater . (), while Figure 1(b) is constructed from the VWH
Italian matched employer-employee dataset, as detailed in Section 5.1.

Figure 1.1 provides two empirical examples of complex outcome missingness
patterns under which the methods developed in this paper can be applied. As an
example of an event study setting, Figure 1(a) displays the control potential
outcome observation pattern from Ater, Shany, Ross, Turkel Vasserman (), which
evaluates the effect of a government-run usage-based congestion pricing
incentive on the driving behavior of several thousand Israeli drivers from
January of 2020 through July of 2021. As can be seen in Figure 1(a), drivers
were recruited for observation in a staggered fashion over time and received
treatment after a non-random 20 weeks of monitoring. As an example of a
bipartite matching setting, Figure 1(b) shows which of a subset of Italian
workers in the Veneto Worker Histories (VWH) dataset received wages in each of
seven provinces and three firm types within each province between 1991 and
2001;444 The VWH dataset was developed by the Economics Department in
Università Ca’ Foscari Venezia under the supervision of Giuseppe Tattara, and
can be accessed at https://www.frdb.org/en/dati/dati-inps-carriere-lavorative-
in-veneto/. see Section 5.1 for details.

In both cases, outcome missingness patterns are credibly determined by
unobserved characteristics of units that also affect outcomes. In Figure
1(a)’s setting, it is possible that drivers who would benefit more from the
treatment like remote workers joined the study earlier, while in Figure 1(b)’s
setting, workers likely live in the regions with job opportunities that would
pay them the highest wages. In addition, the TWFE model’s restriction that
unit-level heterogeneity does not affect how outcomes differ within units is
plausibly unrealistic in both examples. In Figure 1(a)’s setting, drivers in
remote-work-compatible occupations plausibly changed their commuting patterns
differently in response to COVID-19 lockdowns than drivers whose occupations
required continued in-person work. In Figure 1(b)’s setting, if industries are
unequally distributed across space, workers whose skills are
disproportionately valued by the industries in some region might receive
higher wages there than if they lived in a region without those industries.
Finally, in both examples, only a small number of outcomes are observed per
unit relative to the number of units in each sample.

Our identification strategy for counterfactual outcome means in short panels
with missing outcomes is based on a novel type of matrix we call an
_Aggregated Projection Matrix (APM)_. Surprisingly, any basis for an APM’s
null space serves as a collection of the factor vectors corresponding to all
outcomes with respect to the same basis. To construct an APM, we first group
units into cohorts that share the same sets of observed outcomes. Next, we use
the observations for each cohort to recover the collections of factors
corresponding to each cohort’s observed outcomes up to cohort-specific bases
using any of the myriad of existing approaches for identifying factors in
short panels without missing data, e.g. those discussed in Section 4.3.
Finally, we aggregate the projection matrices onto the spans of the matrices
whose non-zero rows are the recovered cohort-specific factors in a particular
manner into an APM. Having identified the factors corresponding to all
outcomes via the APM’s null space, we then leverage an idea proposed in
various forms in Imbens . (2021), Brown Butts (2022), and Agarwal, Shah Shen
(2023) to recover cohort-specific outcome means by learning a linear
relationship between the factors corresponding to the observed and missing
outcomes in the target cohort.

We prove that our approach identifies the collection of factors corresponding
to every outcome so long as a particular graph is connected, where the graph
is constructed based on the overlap in observed outcomes between pairs of
cohorts. Importantly, our method accommodates more general missingness
patterns than other approaches for short panels, and it automatically stitches
together different pieces of identifying information each used in isolation by
existing methods. It also does not make assumptions on how units select into
cohorts based on their low-dimensional unobserved confounders. Moreover,
unlike approaches designed for long panels with many observed outcomes per
unit, our approach does not require exact recovery of those unobserved
confounders, and is thus not vulnerable to the incidental parameter problem
(Neyman Scott, 1948; T. Lancaster, 2000).

We also translate this identification strategy into a plug-in estimator. In
particular, we compute an estimated APM using estimates of cohort-specific
factors constructed using the data corresponding to each cohort. We then use
the rows of the matrix of eigenvectors corresponding to the smallest
eigenvalues of the estimated APM as estimates of the factor vectors
corresponding to all outcomes. As such, our estimator is quite simple to
compute. In an asymptotic regime in which the number of outcomes remains fixed
as the cross-sectional dimension of the panel grows, we show that this
estimator is consistent and asymptotically normal, and that a weighted
bootstrap procedure provides valid asymptotic inference. These results rely on
an exact, first-order expansion of the operator mapping a symmetric matrix
into the projection matrix onto the space spanned by some subset of its
eigenvectors. To our knowledge, this expansion is new.

Finally, we demonstrate the empirical performance of our estimator via a semi-
synthetic simulation study based on the VWH dataset of wages earned by workers
at different types of firms in the Veneto region of Italy. To define outcomes,
we cluster firms within each province into three types by their weekly wage
distributions as in Bonhomme . (2019) and define a worker’s outcome
corresponding to a given type of firm in a given province as the average
weekly wage they would earn were they to work for that type of firm in that
province. Importantly, the outcome missingness pattern for this setting, which
we illustrate in Figure 1(b), is complex enough to preclude counterfactual
outcome mean estimation using existing factor-model-based methods.

To assess the accuracy of our factor-model-based procedure relative to an
estimator based on a TWFE model of counterfactual outcome means, we mask an
observed outcome for some cohort of units in our data and resample units from
these masked data to construct synthetic datasets. We then compute various
error metrics of both estimators of the artificially hidden cohort outcome’s
mean across resampled synthetic datasets. Across many masked cohort outcomes,
our procedure consistently delivers outcome mean estimates with lower bias and
root mean squared error than the TWFE-based estimator. In line with Bonhomme .
(2019), our simulation results suggest that complementarities between workers
and firms do affect wages.

###  1.1 Related Work

Recently, a variety of methods have been proposed that can be used to estimate
and conduct inference on a target cohort outcome mean by aggregating accurate
imputations of individual-level factor structure estimates for units in the
target cohort constructed using data on units’ observed outcomes, e.g.
Gobillon Magnac (2016); Xu (2017); Athey . (2021); J. Bai Ng (2021);
Arkhangelsky . (2021); Ben-Michael . (2022); Xiong Pelger (2023); Agarwal,
Dahleh . (2023); Chernozhukov . (2023); Farias . (2021); Choi Yuan (2023);
Choi . (2023); Arkhangelsky Hirshberg (2023); Fernández-Val . (2021). These
methods require the number of outcomes and maximum number of observed outcomes
per cohort to grow as the number of units grows, since outcome-specific factor
vectors can be estimated consistently using variation across many units and
unit-specific loading vectors can be estimated consistently using variation
across many outcomes. However, under our asymptotic regime of interest in
which only a finite number of outcomes are observed per unit, we cannot
recover each unit’s loadings without bias that persists as the number of units
grows and does not average out to zero across units in the population,
precluding averages across imputations. Further, this bias can contaminate the
outcome-specific factor estimates themselves even if the cross-sectional panel
dimension is large (Neyman Scott, 1948; T. Lancaster, 2000). On the other
hand, when the missingness pattern is complicated, a restriction on the
missingness mechanism, such as missing at random, or knowledge of it, such as
the probability of jointly observing each pair of outcomes, is often needed
even in long panels (Xiong Pelger, 2023). However, this information is
typically unavailable and hard to come by in settings like those illustrated
in Figure 1.1.

Instead, we build on an approach suggested in various forms in Imbens .
(2021), Brown Butts (2022), and Agarwal, Shah Shen (2023): given just factor
vectors corresponding to all outcomes, we can construct a linear combination
of the target cohort’s observed outcome means that equals the target cohort
outcome mean. Imbens . (2021) call any such linear combination a bridge
function to highlight the connections they make between this setting and the
proximal causal inference literature (Miao ., 2018; Deaner, 2018).

The bridge function-based identification strategy was originally developed in
the context of settings with a block missingness pattern, as illustrated in
Figure 1(a). In these settings, there exists a “reference” cohort of units for
whom all outcomes are observed, so any method for identifying and estimating a
factor model without missing data can be applied to the data from the
reference cohort to recover the factor vectors for all outcomes with respect
to the same basis; we discuss several such methods in Section 4.3.

The bridge function approach can still be applied to identify some cohort
outcome means under even more general missingness patterns so long as a
reference cohort exists. In event study settings with staggered treatment
adoption, so long as there is a “never-treated” group of units for whom _all_
control potential outcomes are observed,555Athey . (2021) also require a
nontrivial share of never-treated units for their estimator to be consistent
in the staggered treatment adoption setting. these units can be used as a
reference cohort, and the approaches developed in Callaway Karami (2023),
Brown Butts (2022), Brown . (2023), and Callaway Tsyawo (2023) identify and
yield consistent estimates of outcome means when the number of outcomes
remains fixed as the cross section’s size grows. Agarwal, Dahleh . (2023) show
that the bridge function approach can be applied to identify some cohort
outcome means in even more general settings like the one illustrated in Figure
1(b) by finding block missingness patterns embedded within the broader outcome
missingness pattern.

However, several important challenges remain in estimating and conducting
valid inference on cohort outcome means under factor models that are common in
empirical settings such as those illustrated in Figure 1.1. First, many
embedded block missingness patterns may exist within a broader missingness
pattern to identify some cohort outcome means, and it is unclear how one
should combine information gleaned from each of them in a computationally
efficient manner to improve estimate precision. Second and more importantly,
embedded block missingness patterns cannot always be found to identify all
cohort outcome means. Our identification and estimation approach automatically
combines information about factor vectors from all cohorts to identify and
estimate cohort outcome means without requiring embedded block missingness
patterns to exist.

The graph connectivity assumption underlying our identification argument bears
resemblance to those made in several papers in the rich literature on fixed-
effect-like models of bipartite match outcomes under strict exogeneity, e.g.
Abowd . (2002); Jochmans Weidner (2019); Bonhomme . (2019); Hull (2018); see
Bonhomme (2020) for a review. As it happens, we show in Appendix D.2 that,
while our connectivity requirement enables identification for general numbers
of factors and loadings, in the special case where the factors and loadings
are unidimensional, our assumption is essentially equivalent to theirs.

Our consistency and asymptotic linearity results are also related to a large
literature applying perturbation-theoretic results to characterize the
concentration of eigenspaces of random matrices.666See Yu . (2015) for clear
and concise statements of several results from the literature useful for
statistical applications, and Z. Bai Silverstein (2010) for a textbook
treatment To characterize the asymptotic distributions of eigenspaces of
estimated matrices as we do in this paper, the concentration guarantees
implied by canonical zeroth-order approximation bounds like the Davis–Kahan
theorem (see Yu . (2015) for a convenient version) are too coarse. In our
paper, we instead characterize the asymptotic distributions of estimated
eigenspaces and the estimators on which they depend by deriving a non-
asymptotic, first-order eigenspace projection operator expansion based on
Kato’s integral (Kato, 1949). Several other papers have also applied Kato’s
integral to derive concentration results for eigenspaces of random matrices
with (approximately) independent entries, e.g. Oliveira (2010); Mao . (2021);
Lei (2020). Relatedly, Simons (2023) applies an asymptotic linearization of
the eigenspace operator from Jg. Sun (1991) to derive asymptotically valid
hypothesis tests concerning eigenspaces of an estimated, potentially non-
symmetric matrix.

##  2 Setup and Intuition

###  2.1 Setup

We describe our setting of interest formally as follows. Researchers observe a
large, i.i.d. sample of $N$ units, and each unit has $T$ outcomes associated
with them, where $Y_{it}^{*}$ denotes unit $i$’s outcome $t$. Importantly, not
all outcomes are observed for each unit. To describe which outcomes are
observed, we group units into $C$ _cohorts_ , where $C_{i}$ denotes unit $i$’s
cohort. For the units in cohort $c$, we only observe the outcomes with indices
$t$ in a subset of all outcome indices
$\mathcal{T}_{c}\subseteq\\{1,\dotsc,T\\}$. To distinguish between observed
and missing outcomes, we define $Y_{it}=Y_{it}^{*}$ if outcome $t$ is observed
for unit $i$ (i.e. if $t\in\mathcal{T}_{C_{i}}$) and $Y_{it}=\emptyset$
otherwise. Given these unbalanced panel data, researchers are often interested
in estimating aggregations of cohort-level outcome means:

$\mu_{ct}\coloneqq\operatorname{\mathbb{E}}[Y_{it}^{*}~{}|~{}C_{i}=c].$ (2.1)

We refer to $\mu_{ct}$ as a counterfactual outcome mean since outcome $t$ may
not be observed for the units in cohort $c$. To model the small number of
observed outcomes per unit in our asymptotic theory in Section 4.2, we will
focus on the case where the number of outcomes $T$ remains bounded as the
sample size $N$ grows.

The setup described above encapsulates several kinds of causal panel data
analyses. In event study settings like Figure 1(a)’s, $Y_{it}^{*}$ is unit
$i$’s control potential outcome had they not yet been treated by period $t$,
units belong to cohort $c$ if they first received treatment at time $c$, and
researchers typically estimate aggregations of average treatment effects on
the units in each cohort like dynamic effects across post-treatment
periods.777See L. Sun Abraham (2021) and Callaway Sant’Anna (2021) for more
in-depth discussions of this model of event study settings with staggered
treatment adoption, and De Chaisemartin d’Haultfoeuille (2020) for treatment
effect definitions relevant to settings in which treatment is not an absorbing
state. In bipartite matching settings like Figure 1(b)’s, $Y_{it}^{*}$ is
“row-type” unit $i$’s outcome if matched with “column-type” unit $t$, row-type
units belong to the same cohort if they were matched with the same set of
column-type units, and researchers typically estimate aggregations of
$\mu_{ct}$s that characterize the degree to which row-type and column-type
unit heterogeneity affect differences in average observed outcomes across
column-type units.888In Appendix D.1, we discuss a variant of the
decomposition proposed in Finkelstein . (2016) for this purpose that only
requires estimates of the outcome means identified in this paper. We leave
identifying the higher-order outcome moments required by outcome variance
decompositions in the literature (e.g. Abowd . (1999)) to future work.

We assume outcomes $Y_{it}^{*}$ are generated according to a rank-$r$ factor
model, i.e. that outcomes are determined by the inner product of a fixed
vector of outcome-specific factors
$\gamma_{t}\in\operatorname{\mathbb{R}}^{r}$ and a vector of unit-specific
loadings $\lambda_{i}\in\operatorname{\mathbb{R}}^{r}$ plus an error term
$\operatorname{\varepsilon}_{it}$ that has zero mean given unit $i$’s loadings
and cohort membership $C_{i}$:

$\displaystyle
Y_{it}^{*}=\gamma_{t}^{\prime}\lambda_{i}+\operatorname{\varepsilon}_{it},\quad\operatorname{\mathbb{E}}[\operatorname{\varepsilon}_{it}~{}|~{}\lambda_{i},C_{i}]=0.$
(2.2)

As in other fixed-effect-like models like the TWFE model or the static model
studied in Bonhomme . (2019), the assumption that
$\operatorname{\varepsilon}_{it}$ is mean-independent of cohort membership
$C_{i}$ conditional on loadings $\lambda_{i}$ implies that only the
coordinates of $\lambda_{i}$ can serve as confounders that affect both units’
cohort memberships and their outcomes. By considering the factor vectors
$\gamma_{t}$ to be fixed, we essentially condition our inference on the
factors $\\{\gamma_{t}:t=1,\dotsc,T\\}$. If we consider the factors
$\gamma_{t}$ to be random, the residual mean condition in (2.2) can be read as
$\operatorname{\mathbb{E}}[\operatorname{\varepsilon}_{it}~{}|~{}\lambda_{i},C_{i},\gamma_{1},\dotsc,\gamma_{T}]=0$.

Throughout this paper, we assume the rank of the factor model $r<T$ is
known.999Papers that use empirical strategies based on TWFE models or parallel
trends assumptions are also making a choice about the dimension of unobserved
confounders. In keeping with other methods based on fixed-effect-like models
of short panels with unobserved heterogeneity (e.g. Bonhomme . (2019)), we
suggest that researchers using our method assess the robustness of their
results to different values of $r$, insofar as the cohort outcome means they
would like to estimate are identified given that choice of $r$, as discussed
in Section 4.1.101010The problem of determining the number of factors $r$ from
short panel data without assuming a the distributions of residuals
$\operatorname{\varepsilon}_{it}$ belong to a known parametric family has only
recently been studied to our knowledge (Ahn ., 2013; Fortin ., 2022, 2023),
unlike the well-established literature on determining $r$ in panels for which
both $N$ and $T$ are large (J. Bai Ng, 2002; Onatski, 2009; Ahn Horenstein,
2013; Gagliardini ., 2019). The simulation study in Section 5.2 that assesses
the performance of our estimation procedure on held-out cohort outcomes also
provides a blueprint for assessing the relative performance of our estimator
given different choices of $r$.

###  2.2 Intuition

Broadly, our approach for identifying and estimating a target cohort outcome
mean $\mu_{c^{*}t^{*}}$ proceeds in two stages. First, we recover all of the
factor vectors across outcomes with respect to a common basis. Then, we use
the factor vectors corresponding to the target cohort’s observed outcomes and
the target outcome to extrapolate from the target cohort’s observed outcome
means to the target outcome mean. To highlight the additional challenge posed
by the first stage, we proceed by first providing intuition for the second
stage in a setting where the first stage happens to be straightforward.

(a) Block Missingness Pattern

(b) Three-Cohort Staircase Pattern

Figure 2.1: This figure illustrates two stylized outcome missingness patterns,
where, similarly to Figure 1.1, a coordinate is colored dark red if the
outcome corresponding to the coordinate’s column is observed for the unit
corresponding to the coordinate’s and is colored light red otherwise. We use
numbers to indicate the cohorts of units for whom different blocks of outcomes
are observed.

In particular, we focus on a simple, three-cohort outcome missingness pattern
illustrated in Figure 1(b) and consider estimating a target cohort outcome
mean $\mu_{c^{*}t^{*}}$, where $c^{*}=1$ is the target cohort and
$t^{*}=t_{1}$ is the target outcome. Suppose also that we have access to the
factors $\gamma_{t}$ for all outcomes. Then, under regularity conditions,
there exist many vectors of coefficients $\beta$ we can construct such that a
linear combination of the factor vectors for the target cohort $c^{*}$’s
observed outcomes with those coefficients equals the factor vector
$\gamma_{t^{*}}$ corresponding to the target outcome $t^{*}$:

$\sum_{t\in\mathcal{T}_{c^{*}}}\beta_{t}\gamma_{t}=\gamma_{t^{*}}.$ (2.3)

As it happens, the same linear combination of the target cohort’s observed
outcome means $\operatorname{\mathbb{E}}[Y_{it}~{}|~{}C_{i}=c^{*}]$ for
$t\in\mathcal{T}_{c^{*}}$ will equal the target cohort’s target outcome mean
(Agarwal, Shah Shen, 2023; Imbens ., 2021; Brown Butts, 2022):

$\sum_{t\in\mathcal{T}_{c^{*}}}\beta_{t}\underbrace{\operatorname{\mathbb{E}}[Y_{it}~{}|~{}C_{i}=c^{*}]}_{=\operatorname{\mathbb{E}}[Y_{it}^{*}~{}|~{}C_{i}=c^{*}]}\overbrace{=}^{\text{by
\eqref{eq:factor_model}}}\sum_{t\in\mathcal{T}_{c^{*}}}\beta_{t}\gamma_{t}^{\prime}\operatorname{\mathbb{E}}[\lambda_{i}~{}|~{}C_{i}=c^{*}]\overbrace{=}^{\text{by
\eqref{eq:bridge_function_extrapolation}}}\gamma_{t^{*}}^{\prime}\operatorname{\mathbb{E}}[\lambda_{i}~{}|~{}C_{i}=c^{*}]\overbrace{=}^{\text{by
\eqref{eq:factor_model}}}\mu_{c^{*}t^{*}}.$ (2.4)

In keeping with Imbens . (2021), we refer to any linear transformation
$\beta^{\prime}$ where $\beta$ satisfies (2.3) as a _bridge function_.

Turning to the first stage of our approach, we now discuss how to recover the
factors $\gamma_{t}$ used to construct bridge functions $\beta^{\prime}$ via
solutions to (2.3). A well known fact in the rich literature on factor models
is that, even in settings where all outcomes are observed, only a common
linear transformation $Q\gamma_{t}$ of each factor vector $\gamma_{t}$ can be
identified, where $Q$ is an _unknown_ , $r\times r$ invertible matrix we refer
to as a basis.111111This basis indeterminacy is inevitable because, for any
invertible $Q$, the distributions of $Y_{it}^{*}$ and
$\gamma_{t}^{\prime}QQ^{-1}\lambda_{i}+\operatorname{\varepsilon}_{it}$ are
observationally equivalent; see e.g. Anderson Rubin (1956) for an early
reference and Anderson (2009) for a textbook treatment.

Luckily, the bridge function approach renders this basis indeterminacy
immaterial when our target cohort is $c^{*}=1$ and our target outcome is
$t^{*}=t_{1}$. To see why, note that in this case, there exists a “reference”
cohort $2$ of units for whom we observe _both_ the target outcome $t_{1}$ and
some of the outcomes also observed for the units in the target cohort
(Agarwal, Dahleh ., 2023; Brown Butts, 2022). In other words, the set of units
and outcomes inside the black dashed rectangle in Figure 1(b) constitutes a
block missingness pattern like the one shown in Figure 1(a) embedded within
the broader outcome missingness pattern in Figure 1(b). Then, under a variety
of additional assumptions discussed in Section 4.3, a myriad of methods can
still be applied using just the data from cohort $2$ to recover transformed
factor vectors $Q\gamma_{t_{1}}$ corresponding to the target outcome and
$Q\gamma_{t}$ for $t\in\mathcal{T}_{1}\cap\mathcal{T}_{2}$ corresponding to
the overlapping observed outcomes between cohorts $1$ and $2$, where $Q$ is a
common basis. Since the set of valid bridge functions $\beta^{\prime}$
satisfying (2.3) is invariant to multiplying all of the factor vectors
$\gamma_{t}$ by the _same_ basis matrix $Q$, the second stage of our
identification approach remains unaffected.

Unfortunately, when our target outcome is instead $t^{*}=t_{2}$, no reference
cohort exists for whom we observe both the target outcome $t_{2}$ and any
observed outcomes for our target cohort $c^{*}=1$. As such, there is no subset
of the data we can directly use to recover the factors corresponding to the
target outcome and the target cohort’s observed outcomes with respect to the
same basis.121212The lack of a reference cohort also precludes the imposition
of a common normalization that ensures the factor vectors are uniquely
defined; see J. Bai Ng (2013) for a detailed discussion of such
normalizations. Instead, again under different sets of additional assumptions
(see Section 4.3), we can only use the observed outcomes of units in each
cohort $c$ to recover transformed factor vectors $\tilde{\gamma}_{ct}\coloneqq
Q_{c}\gamma_{t}$ corresponding to cohort $c$’s observed outcomes
$t\in\mathcal{T}_{c}$ with respect to a _cohort-specific_ basis
$Q_{c}$.131313As will be made clear in Section 4, our approach will not
actually require the identification of $\tilde{\Gamma}_{c}=\Gamma
Q_{c}^{\prime}$ for some fixed basis matrices $Q_{c}$; we introduce them here
for simplicity of exposition. In other words, at best, we can recover cohort-
specific factor vectors that are “misaligned.” To be able to find bridge
functions that satisfy (2.3), we need to construct an aligned set of factor
vectors expressed with respect to a common basis that correspond to both the
target outcome $t_{2}$ and the target cohort’s observed outcomes
$t\in\mathcal{T}_{1}$.

To describe our solution to this factor vector misalignment problem, we define
some additional notation. First, stack the factor vectors $\gamma_{t}$ row-
wise into a $T\times r$ matrix $\Gamma$, and stack the cohort-specific
transformed factor vectors $\tilde{\gamma}_{ct}$ recovered from the data for
cohort $c$ into the $T\times r$ matrix $\tilde{\Gamma}_{c}$, where the $t$-th
row of $\tilde{\Gamma}_{c}$ equals $\tilde{\gamma}_{ct}=Q_{c}\gamma_{t}$ if
outcome $t$ is observed for the units in cohort $c$ and a vector of zeros
otherwise. In addition, let $E_{c}$ be the $T\times T$ diagonal matrix whose
$t$-th diagonal entry is one if outcome $t$ is observed for cohort $c$, i.e.
$t\in\mathcal{T}_{c}$, and zero otherwise. Finally, for any matrix $M$, we
denote the projection matrix onto $M$’s column space by $\Pi(M)\coloneqq
M(M^{\prime}M)^{+}M^{\prime}$, where $(M^{\prime}M)^{+}$ is the Moore-Penrose
pseudo-inverse of $M^{\prime}M$.

Armed with this notation, we can now describe the two insights that underpin
our solution to the factor vector alignment problem introduced above. First,
since the column space of $\tilde{\Gamma}_{c}=E_{c}\Gamma Q_{c}^{\prime}$ does
not depend on $Q_{c}$, the projection matrix $\Pi(\tilde{\Gamma}_{c})$ is the
same as the projection matrix $\Pi(E_{c}\Gamma)$ onto the column space of the
matrix $E_{c}\Gamma$ whose non-zero rows are exactly the true factor vectors
corresponding to cohort $c$’s observed outcomes.141414By definition, we have
that $\Pi(\tilde{\Gamma}_{c})=E_{c}\Gamma
Q_{c}^{\prime}(Q_{c}\Gamma^{\prime}E_{c}\Gamma
Q_{c})^{-1}Q_{c}\Gamma^{\prime}E_{c}=E_{c}\Gamma(\Gamma^{\prime}E_{c}^{\prime}\cdot
E_{c}\Gamma)^{-1}\Gamma^{\prime}E_{c}=\Pi(E_{c}\Gamma)$. Thus, the projection
matrix $\Pi(\tilde{\Gamma}_{c})$ derived from cohort $c$’s data provides a
unique representation of the available information about the factor vectors
corresponding to cohort $c$’s observed outcomes free from contamination by the
cohort-specific $Q_{c}$.

Given the cohort-specific projection matrices $\Pi(\tilde{\Gamma}_{c})$, we
then must aggregate them in such a way so as to recover the factor vectors
$\gamma_{t}$ defined with respect to some common basis. Along these lines, our
second key observation is that the column space of $\Gamma$ is a subset of the
null space of the matrix
$E_{c}-\Pi(\tilde{\Gamma}_{c})=E_{c}-\Pi(E_{c}\Gamma)$ for each cohort
$c$.151515Since $\Pi(E_{c}\Gamma)=E_{c}\Pi(E_{c}\Gamma)=\Pi(E_{c}\Gamma)E_{c}$
and $E_{c}^{2}=E_{c}$,
$(E_{c}-\Pi(E_{c}\Gamma))\Gamma=E_{c}(I-\Pi(E_{c}\Gamma))E_{c}\Gamma=\operatorname{\mathbf{0}}_{T\times
r}$. To leverage this insight, we define we define the _Aggregated Projection
Matrix_ (APM) operator as follows for any $C$ matrices
$\Gamma_{1},\dotsc,\Gamma_{C}\in\operatorname{\mathbb{R}}^{T\times r}$:

$A(\Gamma_{1},\dotsc,\Gamma_{C})\coloneqq\sum_{c=1}^{C}\left(E_{c}-\Pi(\Gamma_{c})\right).$
(2.5)

Despite its name, $A(\Gamma_{1},\dotsc,\Gamma_{C})$ is not itself a projection
matrix in general. However, since
$A(\tilde{\Gamma}_{1},\dotsc,\tilde{\Gamma}_{c})$ is by definition the sum of
the matrices $E_{c}-\Pi(\tilde{\Gamma}_{c})$ across cohorts $c$, the column
space of $\Gamma$ must lie in the null space of
$A(\tilde{\Gamma}_{1},\dotsc,\tilde{\Gamma}_{c})$ as well. In Section 4.1, we
show that, perhaps surprisingly, so long as there is sufficient overlap
between the observed outcomes across pairs of cohorts, the null space of
$A(\tilde{\Gamma}_{1},\dotsc,\tilde{\Gamma}_{c})$ is in fact _exactly_ the
column space of $\Gamma$. Thus, the rows of any basis matrix for the null
space of $A$ can serve as valid factor vectors when applying the bridge
function identification strategy.161616One might imagine instead explicitly
recovering the matrix $Q_{c_{1}}Q_{c_{2}}^{-1}$ that aligns cohort $c_{2}$’s
factors with cohort $c_{1}$’s by regressing cohort $c_{1}$’s factor vectors
$\tilde{\gamma}_{c_{1}t}$ on cohort $c_{2}$’s factor vectors
$\tilde{\gamma}_{c_{2}t}$ corresponding to the overlapping observed outcomes
$t\in\mathcal{T}_{c_{1}}\cap\mathcal{T}_{c_{2}}$ between the two cohorts as in
J. Bai Ng (2021). One could then chain multiplications of these “pairwise
aligning” matrices together to align cohort-specific factor vectors with each
other, evocative of the identification and estimation strategy proposed in
Hull (2018) under a TWFE-like model of outcomes. However, constructing these
aligning chains is nontrivial in realistic applications with more than a few
cohorts like those illustrated in Figure 1.1, and it is not clear how one
would aggregate information about factors across potentially large numbers of
possible aligning chains. In Section 4.1, we discuss how the procedure we
propose next makes use of information from all available aligning chains to
identify the column space of $\Gamma$ without needing to enumerate them
explicitly.

##  3 Estimation and Inference Procedure

Having provided intuition for the conceptual underpinnings of our approach, we
now introduce our general estimation and inference procedure, which has four
high-level steps. First, we consider the data corresponding to units in each
cohort separately and use them to estimate the factors corresponding to the
observed outcomes for each cohort. Second, we aggregate these cohort-specific
factor estimates into a particular matrix whose eigendecomposition yields
valid estimates of the factors corresponding to all outcomes. Third, given
those factor estimates, for each cohort, we construct the minimum norm bridge
function that consistently extrapolates from that cohort’s observed outcome
means to all target outcome means. Finally, we conduct valid inference on
cohort outcome means by taking advantage of the asymptotic normality of our
estimator.

We now describe our procedure in more detail, beginning with step two outlined
above. We assume for now that for each cohort $c$, we have access to estimates
$\hat{\gamma}_{ct}$ of the factor vectors corresponding to each outcome
$t\in\mathcal{T}_{c}$ observed for cohort $c$. In Section 4.3, we discuss a
variety of methods that exist for constructing valid cohort-specific factor
vector estimates under different sets of additional assumptions beyond the
factor model (2.2). Importantly, most estimators cohort-specific factor
vectors $\tilde{\gamma}_{ct}$ require minimal computational overhead for small
$T$, and they can be computed in parallel across cohorts using subsets of the
data associated with each cohort. We then stack these cohort-specific
estimates into the $T\times r$ matrix $\hat{\Gamma}_{c}$ whose $t$-th row is
$\hat{\gamma}_{ct}$ if $t\in\mathcal{T}_{c}$, i.e. if outcome $t$ is observed
for the units in cohort $c$ and a vector of zeros otherwise.

Given cohort-specific estimated factor matrices
$\hat{\Gamma}_{1},\dotsc,\hat{\Gamma}_{C}$, the second step of our procedure
consists of constructing the estimated APM $\hat{A}\coloneqq
A(\hat{\Gamma}_{1},\dotsc,\hat{\Gamma}_{C})$ and computing its
eigendecomposition. Because $T$ is typically small in settings in which we
envision our method will be applied, this step can be done extremely quickly
using the myriad of optimized eigendecomposition routines available in most
programming languages.171717 Examples include the eigvecs function in Julia’s
LinearAlgebra module and the eigen function available in base R. Let
$\hat{\Gamma}$ be the $T\times r$ matrix whose columns are eigenvectors
corresponding to the $r$ smallest eigenvalues of $\hat{A}$;181818Recall that
eigenvectors are unique up to magnitudes, signs, and permutations between
indices corresponding to repeated eigenvalues. as we show in Section 4.2, row
$t$ of $\hat{\Gamma}$ is a consistent estimate of the factor vector
$\gamma_{t}$ with respect to a particular error metric.

Describing the third step of our procedure requires several more definitions.
Collect unit $i$’s observed outcomes into the vector
$Y_{i}\coloneqq(Y_{i1},\dotsc,Y_{iT})^{\prime}$, and note that although
$Y_{i}$ has undefined entries corresponding to unobserved outcomes, all
entries of the vector $E_{C_{i}}Y_{i}$ are well-defined since the entries
corresponding to the unobserved outcomes are zero. We also let
$N_{c}\coloneqq\sum_{i=1}^{N}\mathbbm{1}\\{C_{i}=c\\}$ denote the number of
units in cohort $c$. We can then define our estimator of the vector
$\mu_{c^{*}}\in\operatorname{\mathbb{R}}^{T}$ whose $t$-th entry is the target
cohort $c^{*}$’s mean value of outcome $t$, $\mu_{c^{*}t}$:

$\hat{\mu}_{c^{*}}\coloneqq\hat{\Gamma}\left(E_{c^{*}}\hat{\Gamma}\right)^{+}\cdot\frac{1}{N_{c^{*}}}\sum_{i=1}^{N}\mathbbm{1}\\{C_{i}=c^{*}\\}E_{C_{i}}Y_{i}.$
(3.1)

Given that $\hat{\mu}_{c^{*}}$ is simply a vector of least-norm solutions to
underdetermined linear equations, it can also be computed efficiently in most
programming languages.191919Given that $T$ is small, the simplest method is to
compute the singular value decomposition $E_{c^{*}}\hat{\Gamma}=USV^{\prime}$,
e.g. using the svd function available in Julia’s LinearAlgebra module or base
R, and then compute
$\hat{\mu}_{c^{*}}=\hat{\Gamma}VS^{+}U^{\prime}\frac{1}{N_{c^{*}}}\sum_{i=1}^{N}\mathbbm{1}\\{C_{i}=c^{*}\\}E_{C_{i}}Y_{i}$,
where $S^{+}$ is simply the diagonal matrix $S$ with its positive diagonal
entries replaced by their inverses.

Often, researchers are actually interested in estimating a known, vector-
valued function $\theta\coloneqq h(\mu,\eta)$ of cohort outcome means
$\mu\coloneqq(\mu_{1}^{\prime},\dotsc,\mu_{C}^{\prime})^{\prime}$, along with
a vector of nuisance parameters $\eta$ that are estimable from the data. For
example, in event study settings, it is common to report average effects of a
treatment across different numbers of time periods relative to units’
treatment times. As introduced in Callaway Sant’Anna (2021) and L. Sun Abraham
(2021) and reviewed in Example D.1 of Appendix D.1, we can express these
estimands as weighted averages of the differences between each cohort’s
average control potential outcome means $\mu_{c}$ and their observed, treated
potential outcome means (part of $\eta$), where the weights are determined by
the relative sizes of each cohort (also part of $\eta$).

In the context of bipartite match data on patients’ health outcome when living
in different geographic areas, Finkelstein . (2016) suggest an approach to
attributing shares of the differences in average health outcomes across
regions to people and places that can be expressed in the form $h(\mu,\eta)$,
albeit based on a TWFE model of match outcomes. In Example D.2 of Appendix
D.1, we discuss nonparametric analogs of their estimands of the form
$h(\mu,\eta)$ that do not depend on a particular model of outcomes, but that
can be estimated under the factor model (2.2) using our procedure.

Given $\theta$’s relevance in applied contexts, as an extension of the third
step of our procedure, we compute a plug-in estimate $\hat{\theta}\coloneqq
h(\hat{\mu},\hat{\eta})$ of $\theta$, where the vector
$\hat{\mu}\coloneqq(\hat{\mu}_{1}^{\prime},\dotsc,\hat{\mu}_{C}^{\prime})^{\prime}$
collects the cohort outcome mean vector estimates $\hat{\mu}_{c}$ defined in
(3.1) across cohorts, and $\hat{\eta}$ is an estimate of the nuisance
parameter vector $\eta$. For convenience, we summarize the three steps of our
estimation procedure in Algorithm 1.

Data: $\\{(C_{i},Y_{i})\\}_{i=1}^{N}$

1 for  _$c\in\\{1,\dotsc,C\\}$_ do

Compute cohort-specific factor estimates $\hat{\gamma}_{ct}$ for
$t\in\mathcal{T}_{c}$ (see Section 4.3)

Construct a $T\times r$ matrix $\hat{\Gamma}_{c}$ with row $t$ equalling
$\hat{\gamma}_{ct}$ if $t\in\mathcal{T}_{c}$, and
$\operatorname{\mathbf{0}}_{r}$ otherwise

end for

2Construct estimated APM $\hat{A}=P(\hat{\Gamma}_{1},\dotsc,\hat{\Gamma}_{C})$
as in (2.5)

3 Compute $T\times r$ matrix $\hat{\Gamma}$ of eigenvectors corresponding to
$\hat{A}$’s $r$ smallest eigenvalues

4 for  _$c\in\\{1,\dotsc,C\\}$_ do

Compute outcome mean estimate vector $\hat{\mu}_{c}$ for cohort $c$ as in
(3.1)

end for

5Compute estimate $\hat{\eta}$ of nuisance parameters $\eta$ necessary for
estimating $\theta$

6 Compute plug-in target parameter estimate
$\hat{\theta}=h(\hat{\mu},\hat{\eta})$, where
$\hat{\mu}=(\hat{\mu}_{1}^{\prime},\dotsc,\hat{\mu}_{C}^{\prime})^{\prime}$

Algorithm 1 Estimation

For the fourth and final step of our procedure, we construct simultaneous
$1-\alpha$ confidence intervals for the $p$ coordinates of $\theta$.202020A
set of confidence intervals for the coordinates of $\theta$ has simultaneous
coverage $1-\alpha$ if the probability that all coordinates of $\theta$ lie
inside their respective intervals simultaneously is at least $1-\alpha$. For
simplicity and numerical stability, we construct these intervals via a
Bayesian bootstrap (Rubin, 1981).212121In principle, any bootstrap procedure
that weights or resamples units could also be used, but to avoid pathological
cases where no units from some cohort are sampled, we prefer weighted
bootstrap procedures. To describe how we construct these intervals, we
introduce more notation. Given a vector of $N$ non-negative weights
$W\coloneqq(W_{1},\dotsc,W_{N})^{\prime}$ that sum to one, we will assume that
both our cohort-specific factor matrix estimators $\hat{\Gamma}_{ct}$ and our
nuisance parameter estimator $\hat{\eta}$ can be adapted to accommodate non-
uniform sampling weights $W$, which we denote $\hat{\Gamma}_{c}(W)$ and
$\hat{\eta}(W)$, respectively; see Section 4.3 for details on weighted
versions of $\hat{\Gamma}_{c}$. We also let $\hat{A}(W)\coloneqq
A(\hat{\Gamma}_{1}(W),\dotsc,\hat{\Gamma}_{C}(W))$ denote the APM constructed
from the weighted, estimated cohort-specific factor matrices
$\hat{\Gamma}_{1}(W),\dotsc,\hat{\Gamma}_{C}(W)$, we let $\hat{\Gamma}(W)$
denote the equivalent of $\hat{\Gamma}$ constructed from the weighted
estimated APM $\hat{A}(W)$, and we let $\hat{\mu}_{c}(W)$ denote the weighted
equivalent of $\hat{\mu}_{c}$ defined as follows:

$\displaystyle\hat{\mu}_{c}(W)\coloneqq\hat{\Gamma}(W)\left(E_{c}\hat{\Gamma}(W)\right)^{+}\sum_{i=1}^{N}\frac{W_{i}\mathbbm{1}\\{C_{i}=c\\}}{\sum_{j=1}^{N}W_{j}\mathbbm{1}\\{C_{j}=c\\}}E_{C_{i}}Y_{i}.$
(3.2)

A weighted version of our plug-in estimator of $\hat{\theta}$ is
$\hat{\theta}(W)\coloneqq h(\hat{\mu}(W),\hat{\eta}(W))$, where $\hat{\mu}(W)$
is the weighted analog of $\hat{\mu}$. Before continuing, we note that under
uniform weights $W_{1}=\dotsc=W_{N}=\frac{1}{N}$, all of the quantities
defined previously in this paragraph equal their unweighted counterparts
defined above.

Data: $\\{(C_{i},Y_{i})\\}_{i=1}^{N}$, number of bootstrap samples $M$.

1 for  _$m\in\\{1,\dotsc,M\\}$_ do

Sample $\xi_{m1},\dotsc,x_{mN}\operatorname{\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}}\mathrm{Exp}(1)$

Construct weight vector $W_{m}\coloneqq(W_{m1},\dotsc,W_{mN})^{\prime}$ with
$W_{mi}\coloneqq\xi_{mi}/\sum_{j=1}^{N}\xi_{mj}$

for  _$c\in\\{1,\dotsc,C\\}$_ do

Compute weighted cohort outcome mean estimate vector $\hat{\mu}_{c}(W_{m})$ as
in (3.2) (using Algorithm 1 with weighted analogs)

end for

Compute weighted nuisance parameter estimates $\hat{\eta}(W_{m})$

Compute weighted target parameter estimate
$\hat{\theta}^{*}_{m}=h(\hat{\mu}(W_{m}),\hat{\eta}(W_{m}))$

end for

2for  _$j\in\\{1,\dotsc,p\\}$_ do

Compute estimate $\hat{\sigma}_{j}$ of $\hat{\theta}$’s standard error, e.g.
as in (3.3)

end for

3Compute estimated critical value $\hat{q}_{1-\alpha}$ as in (3.4)

4 Compute simultaneous $1-\alpha$ confidence intervals $\hat{\mathcal{C}}_{j}$
as in (3.5)

Algorithm 2 Weighted Bootstrap Inference

Given this notation, our inference procedure proceeds as follows. First, for
each iteration $m$ of a large number $M$ of repetitions,222222We recommend
$M\geq 500$. we take $N$ i.i.d. draws $\xi_{m1},\dotsc,\xi_{mN}$ from the
$\mathrm{Exponential}(1)$ distribution,232323Other non-negative distributions
are also possible as long as they satisfy regularity conditions; see Section
3.6.2 in van der Vaart Wellner (1996) for more examples. construct a vector
$W_{m}\coloneqq(W_{m1},\dotsc,W_{mN})^{\prime}$ of $N$ normalized weights
$W_{mi}\coloneqq\xi_{mi}/\sum_{j=1}^{N}\xi_{mj}$, and compute a weighted
target parameter estimate
$\hat{\theta}^{*}_{m}\coloneqq\hat{\theta}(W_{m})$.242424Since $T$ is small,
each computation of $\hat{\theta}(W_{m})$ should be quite fast, as discussed
when describing Algorithm 1 above. Further, $\hat{\theta}(W_{m})$ can be
computed in parallel across iterations $m$, boosting computational efficiency
further. However, in settings where $T\ll N$ but $N$ and $T$ are both very
large, the simplicity of this weighted bootstrap procedure may be outweighed
its computational burden. In such cases, one could instead construct a
multiplier bootstrap inference procedure like the one proposed in Belloni .
(2017) using the influence function expression given in Theorem 4.4, Theorem
4.5, and Corollary D.1. Next, for each coordinate $j\in\\{1,\dotsc,p\\}$ of
$\theta$, we compute an estimate $\hat{\sigma}_{j}$ of the standard error of
$\hat{\theta}_{j}$, e.g.

$\hat{\sigma}_{j}\coloneqq\frac{q_{0.75}(\hat{\theta}^{*}_{1j},\dotsc,\hat{\theta}^{*}_{Mj})-q_{0.25}(\hat{\theta}^{*}_{1j},\dotsc,\hat{\theta}^{*}_{Mj})}{q_{0.75}(Z)-q_{0.25}(Z)},$
(3.3)

where we let $q_{\zeta}(x_{1},\dotsc,x_{M})$ denote the $\zeta$th quantile
across scalars $x_{1},\dotsc,x_{M}$, and we let $q_{\zeta}(Z)$ denote the
$\zeta$th quantile of the standard Gaussian distribution.252525We suggest this
interquartile-range-based estimate of estimator standard errors because it is
more robust to outliers than other standard error estimators like the standard
deviation over weighted bootstrap draws; see the discussion in Remark 3.2 in
Chernozhukov . (2013) for details. Finally, we let $\hat{q}_{1-\alpha}$ denote
the following estimated critical value:

$\hat{q}_{1-\alpha}\coloneqq q_{1-\alpha}(z^{*}_{1},\dotsc,z^{*}_{M}),\quad
z^{*}_{m}\coloneqq\max_{j\in\\{1,\dotsc,p\\}}\left\lvert\hat{\theta}^{*}_{mj}\big{/}\hat{\sigma}_{j}\right\rvert,$
(3.4)

and we define our simultaneous $1-\alpha$ confidence intervals
$\hat{\mathcal{C}}_{j}$ for $\theta_{j}$ across $j=1,\dotsc,p$ as follows:

$\hat{\mathcal{C}}_{j}\coloneqq\left[\hat{\theta}_{j}-\hat{q}_{1-\alpha}\hat{\sigma}_{j},\hat{\theta}_{j}+\hat{q}_{1-\alpha}\hat{\sigma}_{j}\right].$
(3.5)

For convenience, we summarize the steps of our inference procedure in
Algorithm 2.

##  4 Theoretical Properties

###  4.1 Identification Given Cohort-Specific Factors

To highlight the contributions of this paper, we assume for now that we have
identified the column spaces of the cohort-specific factor matrices
$E_{c}\Gamma$:

######  Assumption 4.1.

The projection matrices $\Pi(E_{c}\Gamma)$ are identified for all cohorts
$c=1,\dotsc,C$.

In Section 4.3, we will discuss several sets of additional assumptions
standard in the rich literature on factor models under which Assumption 4.1
holds.

To characterize which outcomes our approach can and cannot identify given
identified cohort specific factor projection matrices $\Pi(E_{c}\Gamma)$, we
define an object we call the _Observed Outcome Overlap Graph_
$\mathcal{G}_{r}$. The graph $\mathcal{G}_{r}$ consists of $C$ vertices, one
corresponding to each cohort, and an edge between two distinct cohorts $c_{1}$
and $c_{2}$ if the space spanned by the factor vectors corresponding to the
two cohorts’ overlapping observed outcomes
$\mathcal{T}_{c_{1}}\cap\mathcal{T}_{c_{2}}$ has at least $r$ dimensions, or
formally, if

$\mathrm{rank}\left(E_{c_{1}}E_{c_{2}}\Gamma\right)=r.$ (4.1)

While (4.1) is a requirement on the unobserved factor matrix $\Gamma$, we note
that a necessary condition that only depends on the known sets of observed
outcomes $\mathcal{T}_{c}$ across cohorts is that cohorts $c_{1}$ and $c_{2}$
have at least $r$ overlapping observed outcomes:
$\lvert\mathcal{T}_{c_{1}}\cap\mathcal{T}_{c_{2}}\rvert\geq r$. As we show in
Appendix D.2, this condition is also sufficient when $r=1$. In Figure E.1, we
illustrate $\mathcal{G}_{1}$ for the two empirical examples introduced in
Section 1.

Our key requirement to identify all of the outcome means for a given cohort
can be stated succinctly in terms of the connectedness of $\mathcal{G}_{r}$:

######  Assumption 4.2.

The observed outcome overlap graph $\mathcal{G}_{r}$ is connected.

Before continuing, three remarks concerning Assumption 4.2 are in order.
First, Assumption 4.2 implies that every cohort must have at least $r$
observed outcomes, and that every factor must affect at least one outcome in
every cohort (a proof is provided in Appendix A.1):

######  Lemma 4.1.

Under Assumption 4.2, $\mathrm{rank}(E_{c}\Gamma)=r$ for all cohorts
$c=1,\dotsc,C$.

A consequence of Lemma 4.1 is that our identification, estimation, and
inference results hold for the subset of units for whom at least $r$ outcomes
are observed, i.e. conditional on $\lvert\mathcal{T}_{C_{i}}\rvert\geq
r$.262626Such an assumption is analogous to restrictions of the samples in
empirical work on bipartite matching settings to the row-type units who are
matched with at least two column-type units.

Second, if $\mathcal{G}_{r}$ is not connected, our identification, estimation,
and inference results apply instead to the subset of cohorts belonging to the
connected component of $\mathcal{G}_{r}$ that contains the target cohort
$c^{*}$ and the subset of outcomes observed for at least one of those cohorts.
Third, we show in Appendix D.2 that when $r=1$, Assumption 4.2 is equivalent
to typical connectedness assumptions made in the literature on fixed-effect-
type models of bipartite network data under strict exogeneity (see Abowd .
(2002); Jochmans Weidner (2019); Bonhomme . (2019); Hull (2018) for examples,
and Bonhomme (2020) for a review).

Given identification of the cohort-specific factor projection matrices
$\Pi(E_{c}\Gamma)$, we denote the population APM by substituting the
population cohort-specific factor matrices $E_{c}\Gamma$ into the APM
definition (2.5) as follows:

$A\coloneqq A(E_{1}\Gamma,\dotsc,E_{c}\Gamma).$ (4.2)

We are now equipped to state our main identification result:

######  Theorem 4.2.

Suppose Assumptions 4.1 and 4.2 hold; then the column space of $\Gamma$ is
identified by the null space of the APM $A$.

The containment of the column space of $\Gamma$ in the null space of the APM
$A$ can be shown succinctly even without Assumption 4.2 (see Footnote 15).
However, the containment of the null space of the APM $A$ in the column space
of $\Gamma$ under Assumption 4.2 requires a more involved argument, which we
provide in Appendix A.2.

Given identification of the column space of the factor matrix $\Gamma$, we can
then identify the target cohort’s outcome means via the bridge function
approach articulated in Section 3:

######  Corollary 4.3.

Suppose Assumptions 4.1 and 4.2 hold, and let $\tilde{\Gamma}$ be any basis
for the null space of the APM $A$. Then the vector $\mu_{c^{*}}$ of outcome
means for cohort $c^{*}$ is identified as follows:

$\mu_{c^{*}}=\tilde{\Gamma}\left(E_{c^{*}}\tilde{\Gamma}\right)^{+}\operatorname{\mathbb{E}}[E_{c^{*}}Y_{i}~{}|~{}C_{i}=c^{*}].$
(4.3)

We provide a proof of Corollary 4.3 in Appendix A.3.

Before continuing, we note that, because we impose the factor model functional
form (2.2), our identification argument for $\mu_{c}$ does not rely on any
support assumptions about the cohort-specific loading distributions
$\lambda_{i}~{}|~{}C_{i}=c$, unlike some other approaches to estimating
counterfactual outcome means in the presence of unobserved heterogeneity. For
example, in the context of event study settings, methods in the Synthetic
Control family like Abadie . (2010) and Arkhangelsky . (2021) typically
require the average target cohort’s loadings to lie in the convex hull of the
loadings of the units in the donor pool. In the context of bipartite match
outcomes, Bonhomme . (2019) require every latent type of row-type unit to be
matched with every latent type of column-type unit. In Section 4.3, none of
the additional sets of assumptions we discuss require overlapping support in
the cohort-specific loading distributions for Assumption 4.1 to hold either.
As such, our approach is robust to large discrepancies in the unobserved
confounding variables across cohorts, so long as the low-rank factor model
(2.2) is correctly specified.

###  4.2 Asymptotic Linearity of $\hat{\mu}_{c}$ Given Cohort-Specific Factors

Having shown how to identify the column space of $\Gamma$ and, by extension,
$\mu_{c}$, we now turn to showing that the plug-in estimator described in
Section 3 yields asymptotically linear estimates of these parameters. First,
we introduce some convenient notation. Let $\mathrm{vec}(M)$ denote the
vectorization of the matrix $M$, i.e. the vector containing the columns of the
matrix $M$ stacked in order, and let $\hat{\operatorname{\mathbb{E}}}_{N}$
denote the expectation operator with respect to the empirical measure
$\hat{P}_{N}$ with sample size $N$.

To highlight the contributions of this paper, as in Section 4.1, we will
assume in this section that we are equipped with asymptotically linear
estimators $\Pi(\hat{\Gamma})$ of the column spaces of $E_{c}\Gamma$,
represented uniquely by $\Pi(E_{c}\Gamma)$:

######  Assumption 4.3.

For each cohort $c=1,\dotsc,C$, there exists a function $\phi_{c}$ of $C_{i}$
and $Y_{i}$ such that the vectorization of the estimated cohort-specific
factor projection matrix $\Pi(\hat{\Gamma}_{c})$ satisfies the following
asymptotic expansion as $N\rightarrow\infty$:

$\sqrt{N}\mathrm{vec}\left(\Pi(\hat{\Gamma}_{c})-\Pi(E_{c}\Gamma)\right)=\sqrt{N}\hat{\operatorname{\mathbb{E}}}_{N}[\phi_{c}(C_{i},Y_{i})]+o_{p}\left(1\right),$
(4.4)

$\operatorname{\mathbb{E}}[\phi_{c}(C_{i},Y_{i})]=\operatorname{\mathbf{0}}_{T}$,
and
$\operatorname{\mathbb{E}}\left[\left\lVert\phi_{c}(C_{i},Y_{i})\right\rVert_{2}^{2}\right]<\infty$.

In Section 4.3, we will discuss several sets of additional assumptions
standard in the rich literature on factor models under which estimators
satisfying Assumption 4.3 exist.

Before continuing, we note that, at first glance, Assumption 4.3 might appear
to be chosen simply to imply asymptotic linearity of the empirical APM
$\hat{A}$ without much effort, since the population APM $A$ and its plug-in
counterpart $\hat{A}$ can be directly expressed as sums of $\Pi(E_{c}\Gamma)$
and $\Pi(\hat{\Gamma}_{c})$ across $c$, respectively. While Assumption 4.3 is
certainly convenient for this reason, perhaps more importantly, it is also a
weaker requirement on the cohort-specific factor matrix estimators
$\hat{\Gamma}_{c}$ than assuming the kind of asymptotic linearity guarantees
typically shown under stronger assumptions in the literature on factor
estimation, as discussed in more detail in Section 4.3 and Appendix D.3.

Next, we can show that the plug-in estimator $\hat{\Gamma}$ of $\Gamma$ is
asymptotically linear in the following sense, where $\otimes$ denotes the
Kronecker product:

######  Theorem 4.4.

Suppose that Assumptions 4.1, 4.2, and 4.3 hold. Then as $N\rightarrow\infty$,

$\displaystyle\sqrt{N}\mathrm{vec}\left(\Pi(\hat{\Gamma})-\Pi(\Gamma)\right)$
$\displaystyle=\sqrt{N}\hat{\operatorname{\mathbb{E}}}_{N}\left[H\sum_{c=1}^{C}\phi_{c}(C_{i},Y_{i})\right]+o_{p}\left(1\right)$
(4.5)

where $H$ is the following $T^{2}\times T^{2}$ matrix:

$H\coloneqq A^{+}\otimes\Pi(\Gamma)+\Pi(\Gamma)\otimes A^{+},$ (4.6)

$\operatorname{\mathbb{E}}[H\sum_{c=1}^{C}\phi_{c}(C_{i},Y_{i})]=\operatorname{\mathbf{0}}_{T^{2}}$,
and $\operatorname{\mathbb{E}}[\lVert
H\sum_{c=1}^{C}\phi_{c}(C_{i},Y_{i})\rVert_{2}^{2}]<\infty$.

We provide a proof of Theorem 4.4 in Appendix A.4. Our proof relies on an
exact, first-order expansion of the operator mapping a symmetric matrix into
the projection matrix onto the space spanned by some subset of its
eigenvectors derived using Kato’s integral (Kato, 1949). Since the projection
matrix onto the space spanned by a matrix’s columns is a unique representation
of that matrix’s column space, this expansion allows us to directly bound the
error incurred by $\Pi(\hat{\Gamma})$ as an estimator of the column space of
$\Gamma$ under minimal assumptions on the eigenvalues of $A$. Because it is
exact, our expansion appears to be new as far as we are aware. Since this
expansion may be of independent interest, we provide a self-contained
description and proof in Appendix B.

Having shown that our plug-in estimator of the column space of $\Gamma$ is
asymptotically linear, we now turn to showing that $\hat{\mu}_{c}$ is
therefore also an asymptotically linear estimator of $\mu_{c}$. To do so, we
require the following additional regularity conditions on cohort sizes and
outcome moments, where
$Y_{i}^{*}\coloneqq(Y_{i1}^{*},\dotsc,Y_{iT}^{*})^{\prime}$ denotes the vector
that collects unit $i$’s outcomes:

######  Assumption 4.4.

For all cohorts $c=1,\dotsc,C$, $\operatorname{\mathbb{P}}(C_{i}=c)>0$ and
$\operatorname{\mathbb{E}}\left[\lVert
Y_{i}^{*}\rVert_{2}^{2}~{}\middle|~{}C_{i}=c\right]<\infty$.

As it happens, we can express $(E_{c}\Gamma)^{+}$ as a known function of
$E_{c}$ and $\Pi(\Gamma)$, which enables us to derive a first-order asymptotic
expansion of our plug-in estimator $\hat{\mu}_{c}$ with respect to
$\Pi(\hat{\Gamma})$.272727Such a result may be no surprise given that the
Moore-Penrose pseudo-inverse of a tall matrix is invariant to an invertible
transformation of the matrix’s rows. This result implies $\hat{\mu}_{c}$ has
exactly our desired properties:

######  Theorem 4.5.

Suppose that Assumptions 4.1, 4.2, 4.3, and 4.4 hold. Then as
$N\rightarrow\infty$,

$\displaystyle\sqrt{N}(\hat{\mu}_{c^{*}}-\mu_{c^{*}})$
$\displaystyle=\sqrt{N}\hat{\operatorname{\mathbb{E}}}_{N}\left[\psi_{c^{*}}(C_{i},Y_{i})\right]+o_{p}\left(1\right),$
(4.7) $\displaystyle\psi_{c^{*}}(C_{i},Y_{i})$ $\displaystyle\coloneqq
M_{c^{*}}H\sum_{c=1}^{C}\phi_{c}(C_{i},Y_{i})+\frac{\mathbbm{1}\\{C_{i}=c^{*}\\}}{\operatorname{\mathbb{P}}(C_{i}=c^{*})}\left(\Gamma\left(E_{c^{*}}\Gamma\right)^{+}E_{c^{*}}Y_{i}-\mu_{c^{*}}\right),$

where

$\displaystyle M_{c^{*}}$
$\displaystyle\coloneqq(2\mu_{c^{*}}^{\prime}E_{c^{*}}R_{c^{*}})\otimes
E_{c^{*}}+(\mu_{c^{*}}^{\prime}E_{c^{*}}(I-R_{c^{*}}))\otimes(I+R_{c^{*}}),$
(4.8) $\displaystyle R_{c^{*}}$
$\displaystyle\coloneqq\Gamma(\Gamma^{\prime}E_{c^{*}}\Gamma)^{-1}\Gamma^{\prime}(I-E_{c^{*}}),$

$\operatorname{\mathbb{E}}[\psi_{c^{*}}(C_{i},Y_{i})]=\operatorname{\mathbf{0}}_{T}$,
and
$\operatorname{\mathbb{E}}\left[\lVert\psi_{c^{*}}(C_{i},Y_{i})\rVert_{2}^{2}\right]<\infty$.

We provide a proof of Theorem 4.5 in Appendix A.5. We note that both
$M_{c^{*}}$ and $R_{c^{*}}$ are invariant to changes of $\Gamma$’s basis.

Given $\hat{\mu}_{c}$’s asymptotic linearity, we can derive the asymptotic
properties of the plug-in estimator $\hat{\theta}$ of some target estimand
$\theta=h(\theta,\eta)$ and the inference procedure for $\theta$ introduced in
Section 3 as corollaries, where $h$ is a known, smooth function of $\theta$
and a vector of nuisance parameters $\eta$ that are consistently estimable at
a parametric rate. In particular, in Appendix D.1, we prove that
$\hat{\theta}$ is asymptotically normal and that our Bayesian-bootstrap-based
simultaneous confidence intervals $\hat{\mathcal{C}}_{j}$ for $\theta$
described in Section 3 have valid simultaneous coverage of the coordinates of
$\theta$ under minimal additional regularity conditions.

In Appendix D.1, we also formally describe two examples of target parameters
that satisfy those regularity conditions. For event study settings, we show
how the dynamic treatment effect parameters discussed in Callaway Sant’Anna
(2021) and L. Sun Abraham (2021) fit into our framework. For studying match
outcomes, we introduce a decomposition of the difference in average observed
match outcomes between two column-type units that attributes shares of that
contrast to differences in the column-type units’ effects on match outcomes
and selection of row-type units into certain matches, analogous to the target
parameter in Finkelstein . (2016).

###  4.3 Identifying and Estimating Cohort-Specific Factors

Having established that, given sufficiently accurate estimates of cohort-
specific factor projection matrices $\Pi(E_{c}\Gamma)$, our approach yields
consistent estimates of and allows us to conduct valid inference on functions
of outcome means, we now turn to constructing such estimates of
$\Pi(E_{c}\Gamma)$. Importantly, without additional assumptions,
$\Pi(E_{c}\Gamma)$ cannot be identified (Anderson Rubin, 1956). As such, in
this section, we provide a non-exhaustive discussion of several common sets of
assumptions from the rich literature on estimating factor models without
missing outcomes that enable identification and estimation of
$\Pi(E_{c}\Gamma)$ when the number of observed outcomes
$\lvert\mathcal{T}_{c}\rvert$ for cohort $c$ remains fixed as $N_{c}$ grows.

#### Uncorrelated, Homoskedastic Outcomes.

First, we discuss a minimal set of assumptions that allows the canonical
Principal Components (PC) estimator to yield consistent estimates of
$\Pi(E_{c}\Gamma)$ when we require $T$ to stay finite as $N$ grows. In
particular, let $V_{c}$ denote the $T\times T$ matrix of cohort $c$’s observed
outcomes’ second moments:

$V_{c}\coloneqq\operatorname{\mathbb{E}}[E_{c}Y_{i}Y_{i}^{\prime}E_{c}~{}|~{}C_{i}=c],$
(4.9)

and let $\tilde{\Gamma}_{c,\text{PC}}$ be any matrix whose columns are
eigenvectors of $V_{c}$ corresponding to $V_{c}$’s $r$ largest eigenvalues.
Further, let $\hat{V}_{c}(W)$ denote a weighted empirical counterpart of
$V_{c}$ for some vector of weights $W=(W_{1},\dotsc,W_{N})^{\prime}$ that are
non-negative and sum to one:

$\hat{V}_{c}(W)\coloneqq
N\hat{\operatorname{\mathbb{E}}}_{N}[W_{i}E_{c}Y_{i}Y_{i}^{\prime}E_{c}~{}|~{}C_{i}=c],$
(4.10)

and define the (weighted) PC estimator of $E_{c}\Gamma$ to be any matrix
$\hat{\Gamma}_{c,\text{PC}}(W)$ whose columns are eigenvectors of
$\hat{V}_{c}(W)$ corresponding to $\hat{V}_{c}(W)$’s $r$ largest
eigenvalues.282828As discussed in Footnote 17 in Section 3,
eigendecompositions can be computed efficiently using optimized routines
available in most programming languages. When
$W_{1}=\dotsc=W_{N}=\frac{1}{N}$, we abuse notation slightly and let
$\hat{V}_{c}\coloneqq\hat{V}_{c}(W)$ and
$\hat{\Gamma}_{c,\text{PC}}\coloneqq\hat{\Gamma}_{c,\text{PC}}(W)$.

Importantly, the PC estimator does not require any auxiliary data like
covariates or instruments to be consistent, and, since it is equivalent to an
eigendecomposition, efficient algorithms exist for its computation. The costs
of its simplicity are the assumptions it requires. First, all loadings must
have nontrivial variances, which rules out the case where all units have the
same loadings. Second, the residuals $\operatorname{\varepsilon}_{it}$ must be
uncorrelated across units $i$ and outcomes $t$, and, while they can have
arbitrarily varying variances across units $i$, they must be constant across
outcomes $t$.

While these assumptions are weaker than the requirement from the classical
factor model literature that $\operatorname{\varepsilon}_{it}$ be drawn
independently from the same mean-zero Gaussian distribution (see e.g. Anderson
(1963)), Theorem 4 in J. Bai (2003) proves that they are both sufficient (as
shown in Connor Korajczyk (1986)) and necessary for the PC estimator to be
consistent when $T$ remains fixed as $N$ grows. In Appendix D.3, we state
these assumptions formally and prove directly that
$\Pi(\hat{\Gamma}_{c,\text{PC}})$ is an asymptotically linear estimator of
$\Pi(E_{c}\Gamma)$ as required by Assumption 4.3 under slightly weaker
conditions than are typically imposed in the literature. To do so, we again
rely on our first-order expansion of the eigenspace projection operator
described in Appendix B.

#### Uncorrelated, Heteroskedastic Outcomes.

If we are instead only willing to believe that the residuals
$\operatorname{\varepsilon}_{it}$ are uncorrelated across units $i$ and
outcomes $t$ but can have arbitrarily heterogeneous variances across both
units and outcomes, then the PC estimator $\Pi(\hat{\Gamma}_{c,\text{PC}})$
will be inconsistent for $\Pi(E_{c}\Gamma)$ if $T$ remains fixed as $N$ grows
(Chamberlain Rothschild, 1983; J. Bai, 2003).292929In fact, the PC estimator
is numerically equivalent to minimizing a least squares objective
$N^{-1}\sum_{i=1}^{N}\left\lVert
E_{c}(Y_{i}-\Gamma\lambda_{i})\right\rVert_{2}^{2}$ with respect to $\Gamma$
and every $\lambda_{i}$ so, as shown in Chamberlain Rothschild (1983),
$\hat{\Gamma}_{c,\text{PC}}$ suffers from the incidental parameter problem
under residual heteroskedasticity (Neyman Scott, 1948; T. Lancaster, 2000).
Instead, the literature suggests estimating $\Pi(E_{c}\Gamma)$ via an
optimization-based approach called Factor Analysis (FA) (see e.g. Chapter 14
of Anderson (2009)). While the global maximizers of the FA objective yield
asymptotically linear estimators of $\Pi(E_{c}\Gamma)$ assuming the population
parameters maximize the population objective (Anderson, 2009), the objective
is non-concave, and even state-of-the-art algorithms for verifiably computing
its global maximum are difficult to scale for even moderate $r$ and $T$
(Bertsimas ., 2017; Khamaru Mazumder, 2019).

Instead, in Appendix D.4, we sketch a new, computationally efficient procedure
for estimating $\Pi(E_{c}\Gamma)$ that can be applied when residuals
$\operatorname{\varepsilon}_{it}$ are uncorrelated but heteroskedastic so long
as at least $2r+1$ outcomes are observed per cohort, like internal-
instruments-based approaches to identifying factors under uncorrelated but
heteroskedastic residuals (Madansky, 1964; Heckman Scheinkman, 1987;
Freyberger, 2018; Harding ., 2022). The procedure is based on the fact that,
if we split the observed outcomes for a given cohort into two disjoint
subsets, the matrix of covariances between pairs of observed outcomes from the
two subsets is determined solely by the factor structure. As such, the top $r$
left singular vectors of this matrix yield estimates of the factor vectors
corresponding to each subset of outcomes with respect to the same basis.
Repeating this procedure for a particular sequence of partitions of a cohort’s
observed outcomes yields factor vector estimates that can then be “stitched
together” via another application of our APM-based factor estimation
procedure. This procedure is guaranteed to recover $\Pi(E_{c}\Gamma)$ with
only $r+1$ eigendecompositions.

#### Other Identifying Assumptions Via Auxiliary Data.

With access to auxiliary data, we can also relax the assumption of
uncorrelatedness of residuals across outcomes using several different
approaches. One approach is to assume that some vector of at least $r$ unit-
and-outcome-specific covariates are also outcome-specific linear functions of
$\lambda_{i}$ plus covariate residuals that are uncorrelated with
$\lambda_{i}$ and $\operatorname{\varepsilon}_{it}$. Then one can apply the
Common Correlated Effects method (Pesaran, 2006; Westerlund ., 2019),
Essential Regression (Bing ., 2022), transfer learning (Duan ., 2023), or the
Diversified Projections method (Fan Liao, 2022) to estimate
$\Pi(E_{c}\Gamma)$. Alternatively, given access to at least $r$ unit-and-
outcome-specific external instruments that are correlated with the coordinates
of $\lambda_{i}$ but not $\operatorname{\varepsilon}_{it}$, one can apply the
GMM estimators proposed in Ahn . (2013) and Robertson Sarafidis (2015) to
estimate $\Pi(E_{c}\Gamma)$.

##  5 Empirical Illustration

###  5.1 Setting

Our empirical illustration of our methods is based on data from the Veneto
Worker Histories (VWH) dataset, which is a dataset derived from the Italian
social security administration containing the full history of weekly wages for
every person who ever lived or worked in any of the seven provinces in the
Veneto region of Italy from 1975 to 2001.303030As noted in Footnote 4, the VWH
dataset was developed by the Economics Department in Università Ca’ Foscari
Venezia under the supervision of Giuseppe Tattara, and can be accessed at
https://www.frdb.org/en/dati/dati-inps-carriere-lavorative-in-veneto/. Each
observation corresponds to a worker working for a firm in a given year, and
contains information about the number of weeks that person worked at that firm
in that year as well as the total wages they were paid for doing so.

We will focus on characterizing the performance of our procedure as an
alternative to TWFE-based methods for assessing to what degree worker’s
locations causally affect their wages as opposed to differences in workers’
observed wages across locations being driven by purely worker sorting, in the
spirit of Card . (2023). As discussed in Example D.2 inspired by Finkelstein .
(2016) in Appendix D.1, an important input to any such analysis are predicted
counterfactual wages for workers had they instead worked in locations we do
not observe them working at in the data. The workhorse methods for
constructing these predicted counterfactual wages are based on the TWFE model
(Abowd ., 1999; Finkelstein ., 2016; Card ., 2023), such methods rule out any
complementarities between workers and certain locations, unlike predicted
average counterfactual wages given by our factor-model-based
method.313131Bonhomme . (2019) note that typical empirical tests in the
literature that claim to not detect match effects in fact have no power under
a variety of models that allow for worker-firm complementarities (see Card .
(2013) for an example deploying such tests). As mentioned in the discussion of
their results, Card . (2013) do find some evidence for the existence of match
effects and suggest that characterizing them further is worth future study.

To assess how our method performs relative to a TWFE regression at predicting
counterfactual mean wages for workers across provinces, we construct a panel
dataset that fits within our setup of interest as follows. First, we restrict
our attention to workers and the 174,470 unique within-Veneto firms at which
they worked between 1991 and 2001 to diminish the impact of long-run secular
trends on our results. To avoid understating the degree to which wages are
determined by workers’ provinces of employment by ignoring within-province
firm heterogeneity (Card ., 2023),323232Borrowing an example from Card .
(2023), if a worker moves from a low-wage province to a high-wage province,
they are likely moving between firms that pay relatively similar wages. As
such, the origin firm in the low-wage province must pay well relative to the
other firms in the low-wage province, while the destination firm must pay
poorly relative to the other firms in the high-wage province. Thus, if we only
let outcomes correspond to the average wages workers could earn in different
provinces, we would wrongly assume that because workers’ wages didn’t change
much when moving across places, the places must not have large effects on
wages. we cluster the between 8,000 and 30,000 firms located in each province
into three types using the $k$-means-based procedure proposed in Bonhomme .
(2019); details can be found in Appendix E.1. We then let $t$ index the pair
of provinces and firm types within provinces, meaning $T=21$, and define the
observed outcome $Y_{it}$ as the natural logarithm of worker $i$’s average
weekly wage earned while working at all firms belonging to province-and-firm-
type pair $t$ between 1991 and 2001. Outcomes $Y_{it}^{*}$ can then be defined
similarly, but correspond to potentially counterfactual average weekly wages
earned at province-and-firm-type pairs $t$ that worker $i$ never actually
worked at during the relevant time period.

After defining the $t$ dimension of our panel dataset, we restrict our sample
to the 327,125 unique workers who moved between province-and-firm-type pairs
at least once between 1991 and 2001, in keeping with the literature on matched
employer-employee data (Abowd ., 1999; Card ., 2013). Figure 1(b) illustrates
the complexity of the outcome missingness pattern generated by this dataset.
We note that no embedded block missingness pattern exists for many target
cohorts and outcomes, precluding the use of methods that rely on the existence
of reference cohorts, as discussed in Section 2.2. In contrast, the observed
outcome overlap graph $\mathcal{G}_{1}$ (which is illustrated in Figure 1(b))
is connected, so Assumption 4.2 holds and $\Pi(\Gamma)$ and $\mu_{c}$ are
identified so long as $r=1$. As such, we focus on comparing our method to
imputations based on TWFE regressions in our empirical evaluation below. We
note that assuming $r=1$ rules out the existence of multidimensional
unobserved confounders still allows for complementarity between worker and
province-by-firm-type effects in determining outcomes, unlike the TWFE model.

###  5.2 Semi-Synthetic Simulation Study of Estimator Performance

To evaluate the performance of our method at predicting counterfactual outcome
means, we conduct a semi-synthetic simulation study based on the dataset whose
construction we described in Section 5.1. First, we choose a target cohort
$c^{*}$ of workers that move between province-by-firm-types at least twice (so
that at least three of their outcomes are observed) and “mask” one of their
outcomes $t^{*}$, treating it as if it were unobserved. Next, we resample with
replacement from each cohort in this masked dataset 100 times, constructing
100 synthetic datasets drawn from the same distribution and with the same
cohort sizes as our masked dataset but without the target outcome $t^{*}$
observed for the units in the target cohort $c^{*}$. For each masked synthetic
dataset, we compute $\hat{\mu}_{c^{*}t^{*}}$ using our method, as well as a
TWFE-based analog. Finally, we evaluating the accuracy of each estimator by
computing the absolute bias, standard error, and root mean squared error over
the 100 estimates constructed using each estimator, treating the actual sample
mean in our original, unmasked dataset of the target outcome $t^{*}$ for the
units in the target cohort $c^{*}$ as the ground truth. To give a
representative sense of each estimator’s accuracy, we repeat this bootstrapped
estimator evaluation procedure across every observed outcome $t^{*}$ for the
15 largest cohorts $c^{*}$ of workers that work for firms of at least three
different province-by-firm types.

Figure 5.1: Plots of CDFs of error metric ratios between our estimator and the
TWFE-based estimator across masked target parameters. The vertical, black,
dashed line denotes the threshold above which our estimator performs worse
than the TWFE estimator with respect to a given error metric, while the
horizontal dotted lines correspond to the shares of target parameters for
which our estimator performs better than the TWFE estimator on the error
metric corresponding to its color, as labeled in the plot’s legend.

We illustrate the distributions of our estimation error metrics across target
cohorts $c^{*}$ and outcomes $t^{*}$ in Figure 5.1. In particular, Figure 5.1
plots the cumulative distribution functions (CDF) of error metric ratios
between our estimator and the TWFE-based estimator across masked target
parameters. For over 75% of target cohorts and outcomes, our factor-model-
based method attains smaller root mean squared error than the TWFE-based
estimator. This frequent better performance is due to our estimator having
lower bias than the TWFE-based estimator for over 80% of target cohorts and
outcomes. The price of using our estimator over a TWFE-based estimator tends
to be slightly higher variance, as indicated by the fact that the TWFE-based
estimator had lower variance for over 90% of the target cohorts and outcomes
we considered. However, the increased variance of our estimator is not enough
to outweigh our estimator’s smaller bias relative to a TWFE-based estimator
when comparing the root mean squared errors of the two methods, as discussed
above. In Figure E.2, we provide scatter plots illustrating the estimators’
absolute error metric values and error metric values relative to the values of
the target cohort outcome means across target cohorts and outcomes. Overall,
this semi-synthetic simulation study provides more evidence that accounting
for complementarities between workers and firms using a method like ours can
yield more accurate estimates of average counterfactual match outcomes in
bipartite matching settings.

##  6 Conclusion

In this paper, we develop a new approach for identifying and estimating
average counterfactual outcomes with short panel data applicable in settings
like event studies and studies of bipartite match outcomes. Relying only on an
eigendecomposition of a new spectral operator, our method produces consistent,
asymptotically normal estimates of means of outcomes generated by a latent
factor model under general outcome missingness patterns as only the cross-
sectional dimension of a panel grows large. Importantly, our procedure
accommodates more general missingness patterns than other approaches for short
panels, and it automatically stitches together different pieces of identifying
information each used in isolation by existing methods. Through our simulation
study based on the Veneto Worker Histories dataset, we also demonstrate the
practicality of our approach in real-world “short” panel data settings.

## References

  * Abadie . (2010) abadie2010syntheticAbadie, A., Diamond, A. Hainmueller, J.  2010\.  Synthetic control methods for comparative case studies: Estimating the effect of California’s tobacco control program Synthetic control methods for comparative case studies: Estimating the effect of california’s tobacco control program.  Journal of the American Statistical Association105490493–505.
  * Abowd . (2002) abowd2002computingAbowd, JM., Creecy, RH. Kramarz, F.  2002\.  Computing person and firm effects using linked longitudinal employer-employee data Computing person and firm effects using linked longitudinal employer-employee data .  Center for Economic Studies, US Census Bureau.
  * Abowd . (1999) abowd1999highAbowd, JM., Kramarz, F. Margolis, DN.  1999\.  High wage workers and high wage firms High wage workers and high wage firms.  Econometrica672251–333.
  * Agarwal, Dahleh . (2023) agarwal2021causalAgarwal, A., Dahleh, M., Shah, D. Shen, D.  2023\.  Causal matrix completion Causal matrix completion.  The Thirty Sixth Annual Conference on Learning Theory The thirty sixth annual conference on learning theory ( 3821–3826).
  * Agarwal, Shah Shen (2023) agarwal2020syntheticAgarwal, A., Shah, D. Shen, D.  2023\.  Synthetic Interventions Synthetic interventions.  arXiv preprint arXiv:2006.07691.
  * Ahn Horenstein (2013) ahn2013eigenvalueAhn, SC. Horenstein, AR.  2013\.  Eigenvalue ratio test for the number of factors Eigenvalue ratio test for the number of factors.  Econometrica8131203–1227.
  * Ahn . (2013) ahn2013panelAhn, SC., Lee, YH. Schmidt, P.  2013\.  Panel data models with multiple time-varying individual effects Panel data models with multiple time-varying individual effects.  Journal of Econometrics17411–14.
  * Anderson (1963) anderson1963asymptoticAnderson, TW.  1963\.  Asymptotic theory for principal component analysis Asymptotic theory for principal component analysis.  The Annals of Mathematical Statistics341122–148.
  * Anderson (2009) anderson2009introductionAnderson, TW.  2009\.  An Introduction to Multivariate Statistical Analysis, 3rd Edition An introduction to multivariate statistical analysis, 3rd edition.  Wiley-Interscience.
  * Anderson Rubin (1956) anderson1956statisticalAnderson, TW. Rubin, H.  1956\.  Statistical Inference in Factor Analysis Statistical inference in factor analysis.  Proceedings of the Third Berkeley Symposium on Mathematical Statistics and Probability: Held at the Statistical Laboratory, University of California, December, 1954, July and August, 1955 Proceedings of the third berkeley symposium on mathematical statistics and probability: Held at the statistical laboratory, university of california, december, 1954, july and august, 1955 ( 1, 111).
  * Andrews . (2008) andrews2008highAndrews, MJ., Gill, L., Schank, T. Upward, R.  2008\.  High wage workers and low wage firms: negative assortative matching or limited mobility bias? High wage workers and low wage firms: negative assortative matching or limited mobility bias?  Journal of the Royal Statistical Society: Series A (Statistics in Society)1713673–697.
  * Angrist Pischke (2009) angrist2009mostlyAngrist, JD. Pischke, JS.  2009\.  Mostly harmless econometrics: An empiricist’s companion Mostly harmless econometrics: An empiricist’s companion.  Princeton university press.
  * Arkhangelsky . (2021) arkhangelsky2021syntheticArkhangelsky, D., Athey, S., Hirshberg, DA., Imbens, GW. Wager, S.  2021\.  Synthetic difference-in-differences Synthetic difference-in-differences.  American Economic Review111124088–4118.
  * Arkhangelsky Hirshberg (2023) arkhangelsky2023largesampleArkhangelsky, D. Hirshberg, D.  2023\.  Large-Sample Properties of the Synthetic Control Method under Selection on Unobservables Large-sample properties of the synthetic control method under selection on unobservables.  arXiv preprint arXiv:2311.13575.
  * Ashenfelter Card (1985) ashenfelter1985usingAshenfelter, O. Card, D.  1985\.  Using the longitudinal structure of earnings to estimate the effect of training programs Using the longitudinal structure of earnings to estimate the effect of training programs.  The Review of Economics and Statistics674648–660.
  * Ater . () ater2023canAter, I., Shany, A., Ross, B., Turkel, E. Vasserman, S.  .  Can usage-based pricing reduce traffic congestion? Can usage-based pricing reduce traffic congestion?  Work in progress.
  * Athey . (2021) athey2021matrixAthey, S., Bayati, M., Doudchenko, N., Imbens, GW. Khosravi, K.  2021\.  Matrix completion methods for causal panel data models Matrix completion methods for causal panel data models.  Journal of the American Statistical Association1165361716–1730.
  * J. Bai (2003) bai2003inferentialBai, J.  2003\.  Inferential theory for factor models of large dimensions Inferential theory for factor models of large dimensions.  Econometrica711135–171.
  * J. Bai Ng (2002) bai2002determiningBai, J. Ng, S.  2002\.  Determining the number of factors in approximate factor models Determining the number of factors in approximate factor models.  Econometrica701191–221.
  * J. Bai Ng (2013) bai2013principalBai, J. Ng, S.  2013\.  Principal components estimation and identification of static factors Principal components estimation and identification of static factors.  Journal of Econometrics176118–29.
  * J. Bai Ng (2021) bai2021matrixBai, J. Ng, S.  2021\.  Matrix completion, counterfactuals, and factor analysis of missing data Matrix completion, counterfactuals, and factor analysis of missing data.  Journal of the American Statistical Association1165361746–1763.
  * J. Bai Ng (2023) bai2021approximateBai, J. Ng, S.  2023\.  Approximate factor models with weaker loadings Approximate factor models with weaker loadings.  Journal of Econometrics.
  * Z. Bai Silverstein (2010) bai2010spectralBai, Z. Silverstein, JW.  2010\.  Spectral analysis of large dimensional random matrices Spectral analysis of large dimensional random matrices ( 20).  Springer.
  * Belloni . (2017) belloni2017programBelloni, A., Chernozhukov, V., Fernandez-Val, I. Hansen, C.  2017\.  Program evaluation and causal inference with high-dimensional data Program evaluation and causal inference with high-dimensional data.  Econometrica851233–298.
  * Ben-Michael . (2022) ben2022syntheticBen-Michael, E., Feller, A. Rothstein, J.  2022\.  Synthetic controls with staggered adoption Synthetic controls with staggered adoption.  Journal of the Royal Statistical Society Series B842351–381.
  * Bertrand . (2004) bertrand2004muchBertrand, M., Duflo, E. Mullainathan, S.  2004\.  How much should we trust differences-in-differences estimates? How much should we trust differences-in-differences estimates?  The Quarterly Journal of Economics1191249–275.
  * Bertsimas . (2017) bertsimas2017certifiablyBertsimas, D., Copenhaver, MS. Mazumder, R.  2017\.  Certifiably optimal low rank factor analysis Certifiably optimal low rank factor analysis.  The Journal of Machine Learning Research181907–959.
  * Bing . (2022) bing2022inferenceBing, X., Bunea, F. Wegkamp, M.  2022\.  Inference in latent factor regression with clusterable features Inference in latent factor regression with clusterable features.  Bernoulli282997–1020.
  * Bonhomme (2020) bonhomme2020econometricBonhomme, S.  2020\.  Econometric analysis of bipartite networks Econometric analysis of bipartite networks.  The econometric analysis of network data The econometric analysis of network data ( 83–121).  Elsevier.
  * Bonhomme . (2023) bonhomme2023muchBonhomme, S., Holzheu, K., Lamadon, T., Manresa, E., Mogstad, M. Setzler, B.  2023\.  How much should we trust estimates of firm effects and worker sorting? How much should we trust estimates of firm effects and worker sorting?  Journal of Labor Economics412291–322.
  * Bonhomme . (2019) bonhomme2019distributionalBonhomme, S., Lamadon, T. Manresa, E.  2019\.  A distributional framework for matched employer employee data A distributional framework for matched employer employee data.  Econometrica873699–739.
  * Borusyak . (2023) borusyak2021revisitingBorusyak, K., Jaravel, X. Spiess, J.  2023\.  Revisiting event study designs: Robust and efficient estimation Revisiting event study designs: Robust and efficient estimation.  Review of Economic Studies.
  * Brown Butts (2022) brown2022generalizedBrown, N. Butts, K.  2022\.  Generalized Imputation Estimators for Factor Models Generalized imputation estimators for factor models.
  * Brown . (2023) brown2023differenceBrown, N., Butts, K. Westerlund, J.  2023\.  Difference-in-Differences via Common Correlated Effects Difference-in-differences via common correlated effects.  arXiv preprint arXiv:2301.11358.
  * Callaway Karami (2023) callaway2023treatmentCallaway, B. Karami, S.  2023\.  Treatment effects in interactive fixed effects models with a small number of time periods Treatment effects in interactive fixed effects models with a small number of time periods.  Journal of Econometrics2331184–208.
  * Callaway Sant’Anna (2021) callaway2021differenceCallaway, B. Sant’Anna, PH.  2021\.  Difference-in-differences with multiple time periods Difference-in-differences with multiple time periods.  Journal of Econometrics2252200–230.
  * Callaway Tsyawo (2023) callaway2023treatment2Callaway, B. Tsyawo, ES.  2023\.  Treatment Effects in Staggered Adoption Designs with Non-Parallel Trends Treatment effects in staggered adoption designs with non-parallel trends.  arXiv preprint arXiv:2308.02899.
  * Card . (2013) card2013workplaceCard, D., Heining, J. Kline, P.  2013\.  Workplace heterogeneity and the rise of West German wage inequality Workplace heterogeneity and the rise of west german wage inequality.  The Quarterly Journal of Economics1283967–1015.
  * Card . (2023) card2023locationCard, D., Rothstein, J. Yi, M.  2023\.  Location, location, location Location, location, location .  National Bureau of Economic Research.
  * Chamberlain Rothschild (1983) chamberlain1983arbitrageChamberlain, G. Rothschild, M.  1983\.  Arbitrage, Factor Structure, and Mean-Variance Analysis on Large Asset Markets Arbitrage, factor structure, and mean-variance analysis on large asset markets.  Econometrica5151281–1304.
  * Chernozhukov . (2013) chernozhukov2013inferenceChernozhukov, V., Fernández-Val, I. Melly, B.  2013\.  Inference on counterfactual distributions Inference on counterfactual distributions.  Econometrica8162205–2268.
  * Chernozhukov . (2023) chernozhukov2021inferenceChernozhukov, V., Hansen, C., Liao, Y. Zhu, Y.  2023\.  Inference for low-rank models Inference for low-rank models.  The Annals of statistics5131309–1330.
  * Chetty Hendren (2018) chetty2018impactsChetty, R. Hendren, N.  2018\.  The impacts of neighborhoods on intergenerational mobility II: County-level estimates The impacts of neighborhoods on intergenerational mobility ii: County-level estimates.  The Quarterly Journal of Economics13331163–1228.
  * Choi . (2023) choi2023inferenceChoi, J., Kwon, H. Liao, Y.  2023\.  Inference for Low-rank Completion without Sample Splitting with Application to Treatment Effect Estimation Inference for low-rank completion without sample splitting with application to treatment effect estimation.  arXiv preprint arXiv:2307.16370.
  * Choi Yuan (2023) choi2023matrixChoi, J. Yuan, M.  2023\.  Matrix Completion When Missing Is Not at Random and Its Applications in Causal Panel Data Models Matrix completion when missing is not at random and its applications in causal panel data models.  arXiv preprint arXiv:2308.02364.
  * Connor Korajczyk (1986) connor1986performanceConnor, G. Korajczyk, RA.  1986\.  Performance measurement with the arbitrage pricing theory: A new framework for analysis Performance measurement with the arbitrage pricing theory: A new framework for analysis.  Journal of financial economics153373–394.
  * Deaner (2018) deaner2018proxyDeaner, B.  2018\.  Proxy controls and panel data Proxy controls and panel data.  arXiv preprint arXiv:1810.00283.
  * De Chaisemartin d’Haultfoeuille (2020) de2020twoDe Chaisemartin, C. d’Haultfoeuille, X.  2020\.  Two-way fixed effects estimators with heterogeneous treatment effects Two-way fixed effects estimators with heterogeneous treatment effects.  American Economic Review11092964–2996.
  * Duan . (2023) duan2023targetDuan, J., Pelger, M. Xiong, R.  2023\.  Target PCA: Transfer learning large dimensional panel data Target pca: Transfer learning large dimensional panel data.  Journal of Econometrics105521.
  * Fan Liao (2022) fan2022learningFan, J. Liao, Y.  2022\.  Learning latent factors from diversified projections and its applications to over-estimated and weak factors Learning latent factors from diversified projections and its applications to over-estimated and weak factors.  Journal of the American Statistical Association117538909–924.
  * Farias . (2021) farias2021learningFarias, V., Li, A. Peng, T.  2021\.  Learning treatment effects in panels with general intervention patterns Learning treatment effects in panels with general intervention patterns.  Advances in Neural Information Processing Systems3414001–14013.
  * Fernández-Val . (2021) fernandez2021lowFernández-Val, I., Freeman, H. Weidner, M.  2021\.  Low-rank approximations of nonseparable panel models Low-rank approximations of nonseparable panel models.  The Econometrics Journal242C40–C77.
  * Finkelstein . (2016) finkelstein2016sourcesFinkelstein, A., Gentzkow, M. Williams, H.  2016\.  Sources of geographic variation in health care: Evidence from patient migration Sources of geographic variation in health care: Evidence from patient migration.  The Quarterly Journal of Economics13141681–1726.
  * Fortin . (2022) fortin2022eigenvalueFortin, AP., Gagliardini, P. Scaillet, O.  2022\.  Eigenvalue tests for the number of latent factors in short panels Eigenvalue tests for the number of latent factors in short panels.  arXiv preprint arXiv:2210.16042.
  * Fortin . (2023) fortin2023latentFortin, AP., Gagliardini, P. Scaillet, O.  2023\.  Latent Factor Analysis in short panels Latent factor analysis in short panels.  arXiv preprint arXiv:2306.14004.
  * Freyberger (2018) freyberger2018nonFreyberger, J.  2018\.  Non-parametric panel data models with interactive fixed effects Non-parametric panel data models with interactive fixed effects.  The Review of Economic Studies8531824–1851.
  * Gagliardini . (2019) gagliardini2019diagnosticGagliardini, P., Ossola, E. Scaillet, O.  2019\.  A diagnostic criterion for approximate factor structure A diagnostic criterion for approximate factor structure.  Journal of Econometrics2122503–521.
  * Ghanem . (2022) ghanem2022selectionGhanem, D., Sant’Anna, PH. Wüthrich, K.  2022\.  Selection and parallel trends Selection and parallel trends.  arXiv preprint arXiv:2203.09001.
  * Gobillon Magnac (2016) gobillon2016regionalGobillon, L. Magnac, T.  2016\.  Regional policy evaluation: Interactive fixed effects and synthetic controls Regional policy evaluation: Interactive fixed effects and synthetic controls.  Review of Economics and Statistics983535–551.
  * Harding . (2022) harding2022estimationHarding, M., Lamarche, C. Muris, C.  2022\.  Estimation of a Factor-Augmented Linear Model with Applications Using Student Achievement Data Estimation of a factor-augmented linear model with applications using student achievement data.  arXiv preprint arXiv:2203.03051.
  * Heckman Scheinkman (1987) heckman1987importanceHeckman, J. Scheinkman, J.  1987\.  The importance of bundling in a Gorman-Lancaster model of earnings The importance of bundling in a gorman-lancaster model of earnings.  The Review of Economic Studies542243–255.
  * Hull (2018) hull2018estimatingHull, P.  2018\.  Estimating treatment effects in mover designs Estimating treatment effects in mover designs.  arXiv preprint arXiv:1804.06721.
  * Imbens . (2021) imbens2021controllingImbens, GW., Kallus, N. Mao, X.  2021\.  Controlling for unmeasured confounding in panel data using minimal bridge functions: From two-way fixed effects to factor models Controlling for unmeasured confounding in panel data using minimal bridge functions: From two-way fixed effects to factor models.  arXiv preprint arXiv:2108.03849.
  * Jackson . (2014) jackson2014teacherJackson, CK., Rockoff, JE. Staiger, DO.  2014\.  Teacher effects and teacher-related policies Teacher effects and teacher-related policies.  Annual Review of Economics61801–825.
  * Jochmans Weidner (2019) jochmans2019fixedJochmans, K. Weidner, M.  2019\.  Fixed-effect regressions on network data Fixed-effect regressions on network data.  Econometrica8751543–1560.
  * Kato (1949) kato1949convergenceKato, T.  1949\.  On the convergence of the perturbation method. I On the convergence of the perturbation method. i.  Progress of Theoretical Physics44514–523.
  * Kato (1980) kato2013perturbationKato, T.  1980\.  Perturbation theory for linear operators Perturbation theory for linear operators ( 132).  Springer Science & Business Media.
  * Khamaru Mazumder (2019) khamaru2019computationKhamaru, K. Mazumder, R.  2019\.  Computation of the maximum likelihood estimator in low-rank factor analysis Computation of the maximum likelihood estimator in low-rank factor analysis.  Mathematical Programming176279–310.
  * P. Lancaster Farahat (1972) lancaster1972normsLancaster, P. Farahat, HK.  1972\.  Norms on direct sums and tensor products Norms on direct sums and tensor products.  mathematics of computation26118401–414.
  * T. Lancaster (2000) lancaster2000incidentalLancaster, T.  2000\.  The incidental parameter problem since 1948 The incidental parameter problem since 1948.  Journal of econometrics952391–413.
  * Lei (2020) lei2020unifiedLei, L.  2020\.  Unified $\ell_{2\rightarrow\infty}$ Eigenspace Perturbation Theory for Symmetric Random Matrices Unified $\ell_{2\rightarrow\infty}$ eigenspace perturbation theory for symmetric random matrices.  arXiv preprint arXiv:1909.04798.
  * Madansky (1964) madansky1964instrumentalMadansky, A.  1964\.  Instrumental variables in factor analysis Instrumental variables in factor analysis.  Psychometrika292105–113.
  * Mao . (2021) mao2021estimatingMao, X., Sarkar, P. Chakrabarti, D.  2021\.  Estimating mixed memberships with sharp eigenvector deviations Estimating mixed memberships with sharp eigenvector deviations.  Journal of the American Statistical Association1165361928–1940.
  * Miao . (2018) miao2018identifyingMiao, W., Geng, Z. Tchetgen Tchetgen, EJ.  2018\.  Identifying causal effects with proxy variables of an unmeasured confounder Identifying causal effects with proxy variables of an unmeasured confounder.  Biometrika1054987–993.
  * Montiel Olea Plagborg-Møller (2019) montiel2019simultaneousMontiel Olea, JL. Plagborg-Møller, M.  2019\.  Simultaneous confidence bands: Theory, implementation, and an application to SVARs Simultaneous confidence bands: Theory, implementation, and an application to SVARs.  Journal of Applied Econometrics3411–17.
  * Neyman Scott (1948) neyman1948consistentNeyman, J. Scott, EL.  1948\.  Consistent estimates based on partially consistent observations Consistent estimates based on partially consistent observations.  Econometrica: Journal of the Econometric Society1–32.
  * Oliveira (2010) oliveira2010concentrationOliveira, RI.  2010\.  Concentration of the adjacency matrix and of the Laplacian in random graphs with independent edges Concentration of the adjacency matrix and of the laplacian in random graphs with independent edges.  arXiv preprint arXiv:0911.0600.
  * Onatski (2009) onatski2009testingOnatski, A.  2009\.  Testing hypotheses about the number of factors in large factor models Testing hypotheses about the number of factors in large factor models.  Econometrica7751447–1479.
  * Pesaran (2006) pesaran2006estimationPesaran, MH.  2006\.  Estimation and inference in large heterogeneous panels with a multifactor error structure Estimation and inference in large heterogeneous panels with a multifactor error structure.  Econometrica744967–1012.
  * Robertson Sarafidis (2015) robertson2015ivRobertson, D. Sarafidis, V.  2015\.  IV estimation of panels with factor residuals Iv estimation of panels with factor residuals.  Journal of Econometrics1852526–541.
  * Rubin (1981) rubin1981bayesianRubin, DB.  1981\.  The bayesian bootstrap The bayesian bootstrap.  The Annals of Statistics130–134.
  * Simons (2023) simons2023inferenceSimons, JR.  2023\.  Inference on eigenvectors of non-symmetric matrices Inference on eigenvectors of non-symmetric matrices.  arXiv preprint arXiv:2303.18233.
  * Jg. Sun (1991) sun1991perturbationSun, Jg.  1991\.  Perturbation expansions for invariant subspaces Perturbation expansions for invariant subspaces.  Linear Algebra and its Applications15385–97.
  * L. Sun Abraham (2021) sun2021estimatingSun, L. Abraham, S.  2021\.  Estimating dynamic treatment effects in event studies with heterogeneous treatment effects Estimating dynamic treatment effects in event studies with heterogeneous treatment effects.  Journal of Econometrics2252175–199.
  * van der Vaart (2000) van2000asymptoticvan der Vaart, A.  2000\.  Asymptotic statistics Asymptotic statistics ( 3).  Cambridge university press.
  * van der Vaart Wellner (1996) vaart1996empiricalvan der Vaart, A. Wellner, JA.  1996\.  Empirical processes Empirical processes.  Weak Convergence and Empirical Processes: With Applications to Statistics. Weak convergence and empirical processes: With applications to statistics.  Springer.
  * Vershynin (2018) vershynin2018highVershynin, R.  2018\.  High-dimensional probability: An introduction with applications in data science High-dimensional probability: An introduction with applications in data science ( 47).  Cambridge university press.
  * Westerlund . (2019) westerlund2019cceWesterlund, J., Petrova, Y. Norkute, M.  2019\.  CCE in fixed-T panels Cce in fixed-t panels.  Journal of Applied Econometrics345746–761.
  * Xiong Pelger (2023) xiong2023largeXiong, R. Pelger, M.  2023\.  Large dimensional latent factor modeling with missing observations and applications to causal inference Large dimensional latent factor modeling with missing observations and applications to causal inference.  Journal of Econometrics2331271–301.
  * Xu (2017) xu2017generalizedXu, Y.  2017\.  Generalized synthetic control method: Causal inference with interactive fixed effects models Generalized synthetic control method: Causal inference with interactive fixed effects models.  Political Analysis25157–76.
  * Yu . (2015) yu2015usefulYu, Y., Wang, T. Samworth, RJ.  2015\.  A useful variant of the Davis–Kahan theorem for statisticians A useful variant of the davis–kahan theorem for statisticians.  Biometrika1022315–323.

##  Appendix A Proofs of Results in the Main Text

###  A.1 Proof of Lemma 4.1

By Assumption 4.2, there must exist at least one other cohort $c^{\prime}$
that is a neighbor of cohort $c$ in the observed outcome overlap graph
$\mathcal{G}_{r}$. Further, since $E_{c^{\prime}}$ is a square matrix, the
column space of $\Gamma^{\prime}E_{c}E_{c^{\prime}}$ must be a subset of the
column space of $\Gamma^{\prime}E_{c}$, which, along with (4.1), implies that

$\displaystyle\mathrm{rank}(\underbrace{E_{c}\Gamma}_{T\times r})$
$\displaystyle=\mathrm{rank}(\Gamma^{\prime}E_{c})$ ($T>r$, rank is min. of
row rank and column rank) (A.1)
$\displaystyle\geq\mathrm{rank}(\Gamma^{\prime}E_{c}E_{c^{\prime}})$ (A.2)
$\displaystyle=r$ $\displaystyle\text{($c$ and $c^{\prime}$ are neighbors, so
\eqref{cond:observed_outcome_overlap_graph_edge_def} holds)}.$ (A.3)

Since $E_{c}\Gamma$ is $T\times r$ with $T>r$, $\mathrm{rank}(E_{c}\Gamma)\leq
r$ must also hold, completing the proof.

###  A.2 Proof of Theorem 4.2

First, as in Footnote 15, we show that for any cohort $c$, the columns of
$\Gamma$ lie in the null space of $E_{c}-\Pi(E_{c}\Gamma)$, in which case
$\mathrm{col}(\Gamma)\subseteq\mathrm{null}(A)$, where $\mathrm{col}(\Gamma)$
denotes the column space of $\Gamma$. Since $E_{c}^{2}=E_{c}$ and
$E_{c}\Pi(E_{c}\Gamma)=\Pi(E_{c}\Gamma)E_{c}=\Pi(E_{c}\Gamma)$, we have that

$(E_{c}-\Pi(E_{c}\Gamma))\Gamma=(E_{c}^{2}-E_{c}\Pi(E_{c}\Gamma)E_{c})\Gamma=E_{c}\cdot(I-\Pi(E_{c}\Gamma))E_{c}\Gamma=\operatorname{\mathbf{0}}_{T\times
r},$ (A.4)

as required.

Next, we show that $\mathrm{null}(A)\subseteq\mathrm{col}(\Gamma)$ also holds.
Define the matrix

$\tilde{A}\coloneqq\begin{bmatrix}E_{1}-\Pi(E_{1}\Gamma)\\\ \vdots\\\
E_{C}-\Pi(E_{C}\Gamma)\end{bmatrix},$ (A.5)

and note that since

$(E_{c}-\Pi(E_{c}\Gamma))^{2}=E_{c}^{2}-E_{c}\Pi(E_{c}\Gamma)-\Pi(E_{c}\Gamma)E_{c}-\Pi(E_{c}\Gamma)^{2}=E_{c}-\Pi(E_{c}\Gamma)$
(A.6)

and $E_{c}-\Pi(E_{c}\Gamma)$ is symmetric, $A=\tilde{A}^{\prime}\tilde{A}$,
implying that $\mathrm{null}(A)=\mathrm{null}(\tilde{A})$.333333It is
straightforward to check that if $v\in\mathrm{null}(\tilde{A})$, then
$v\in\mathrm{null}(A)=\mathrm{null}(\tilde{A}^{\prime}\tilde{A})$. To check
the other direction of inclusion, for any $v\in\mathrm{null}(A)$, we have that
$\operatorname{\mathbf{0}}=Av=\tilde{A}^{\prime}\tilde{A}v$. Further, we have
that
$0=v^{\prime}\operatorname{\mathbf{0}}=v^{\prime}\tilde{A}^{\prime}\tilde{A}v=(\tilde{A}v)^{\prime}(\tilde{A}v)$,
so it must be that $\tilde{A}v=\operatorname{\mathbf{0}}$, as required.

Consider any $v\in\mathrm{null}(\tilde{A})$; by definition, it must be that
for any cohort $c$, $(E_{c}-\Pi(E_{c}\Gamma))v=\operatorname{\mathbf{0}}_{T}$,
which in turn implies that $E_{c}v\in\mathrm{col}(E_{c}\Gamma)$. Thus, there
exists a vector $\omega_{c}\in\operatorname{\mathbb{R}}^{r}$ such that
$E_{c}v=E_{c}\Gamma\omega_{c}$. Then, for any edge $(c_{1},c_{2})$ in the
observed outcome overlap graph $\mathcal{G}_{r}$, since diagonal matrices
commute,

$E_{c_{1}}E_{c_{2}}v=E_{c_{2}}E_{c_{1}}v=E_{c_{2}}E_{c_{1}}\Gamma\omega_{c_{1}}=E_{c_{1}}E_{c_{2}}\Gamma\omega_{c_{1}},\quad
E_{c_{1}}E_{c_{2}}v=E_{c_{1}}E_{c_{2}}\Gamma\omega_{c_{2}},$ (A.7)

implying that
$E_{c_{1}}E_{c_{2}}\Gamma(\omega_{c_{1}}-\omega_{c_{2}})=\operatorname{\mathbf{0}}_{T}$.
Since (4.1) implies that the null space of $E_{c_{1}}E_{c_{2}}\Gamma$ is
trivial and $\omega_{c_{1}}-\omega_{c_{2}}$ lies in it,
$\omega_{c_{1}}=\omega_{c_{2}}$ must also hold. Since $\mathcal{G}_{r}$ is
connected, there exists a path in $\mathcal{G}_{r}$ between any two cohorts,
so $\omega_{1}=\dotsc=\omega_{C}=\omega$ for some common vector
$\omega\in\operatorname{\mathbb{R}}^{r}$. For the $t$-th entry $v_{t}$ of $v$,
since every outcome is observed for at least one cohort, there exists at least
one cohort $c_{t}$ for whom $t\in\mathcal{T}_{c_{t}}$. Letting $e_{t}$ denote
the $t$-th standard basis vector, we have that

$v_{t}=e_{t}^{\prime}E_{c_{t}}\Gamma\omega=\gamma_{t}^{\prime}\omega,$ (A.8)

which, stacking across $t=1,\dotsc,T$, implies that $v=\Gamma\omega$, i.e.
$v\in\mathrm{col}(\Gamma)$. Putting everything together, we have that
$\mathrm{null}(A)=\mathrm{null}(\tilde{A})\subseteq\mathrm{col}(\Gamma)$, as
required.

###  A.3 Proof of Corollary 4.3

First, we show that Assumption 4.2 is sufficient to ensure that
$E_{c^{*}}\tilde{\Gamma}$ is has full rank $r$ (since $T>r$ by assumption),
and thus the $r\times r$ matrix
$\tilde{\Gamma}^{\prime}E_{c^{*}}\tilde{\Gamma}=\left(E_{c^{*}}\tilde{\Gamma}\right)^{\prime}\cdot
E_{c^{*}}\tilde{\Gamma}$ is invertible. Since $\tilde{\Gamma}$ is a basis for
column space of $\Gamma$ by Theorem 4.2, there exists some invertible basis
matrix $Q\in\operatorname{\mathbb{R}}^{r\times r}$ for which
$\tilde{\Gamma}Q^{\prime}=\Gamma$, so
$\tilde{\Gamma}=\Gamma\left(Q^{\prime}\right)^{-1}$. Then

$\displaystyle\mathrm{rank}(E_{c^{*}}\tilde{\Gamma})$
$\displaystyle=\mathrm{rank}\left(E_{c^{*}}\Gamma\left(Q^{\prime}\right)^{-1}\right)$
(A.9) $\displaystyle=\mathrm{rank}(\underbrace{E_{c^{*}}\Gamma}_{T\times r})$
($\left(Q^{\prime}\right)^{-1}$ is full-rank) (A.10) $\displaystyle=r.$ (Lemma
4.1) (A.11)

Next, since for full rank matrices $M$ with more rows than columns,
$M^{+}=(M^{\prime}M)^{-1}M^{\prime}$, we have that

$\displaystyle\tilde{\Gamma}\left(E_{c^{*}}\tilde{\Gamma}\right)^{+}\operatorname{\mathbb{E}}[E_{c^{*}}Y_{i}~{}|~{}C_{i}=c^{*}]$
$\displaystyle=\Gamma Q^{\prime}\cdot(Q\Gamma^{\prime}E_{c^{*}}\Gamma
Q^{\prime})^{-1}Q\Gamma^{\prime}E_{c^{*}}\cdot\operatorname{\mathbb{E}}[E_{c^{*}}Y_{i}~{}|~{}C_{i}=c^{*}]$
(A.12)
$\displaystyle=\Gamma\cdot(\Gamma^{\prime}E_{c^{*}}\Gamma)^{-1}\Gamma^{\prime}E_{c^{*}}\Gamma\operatorname{\mathbb{E}}[\lambda_{i}~{}|~{}C_{i}=c^{*}]$
(A.13)
$\displaystyle=\Gamma\cdot\operatorname{\mathbb{E}}[\lambda_{i}~{}|~{}C_{i}=c^{*}]$
(A.14) $\displaystyle=\mu_{c^{*}},$ (A.15)

as required.

###  A.4 Proof of Theorem 4.4

First, we show that
$\lVert\hat{A}-A\rVert_{\mathrm{op}}=O_{p}\left(N^{-1/2}\right)$:

$\displaystyle\lVert\hat{A}-A\rVert_{\mathrm{op}}$
$\displaystyle\leq\lVert\hat{A}-A\rVert_{F}$ (Equivalence of matrix norms)
(A.16) $\displaystyle=\lVert\mathrm{vec}(\hat{A}-A)\rVert_{2}$ (A.17)
$\displaystyle\leq\sum_{c=1}^{C}\lVert\mathrm{vec}(\Pi(\hat{\Gamma}_{c})-\Pi(E_{c}\Gamma))\rVert_{2}.$
(Definitions of $\hat{A}$ and $A$) (A.18)

By Assumption 4.3,
$\mathrm{vec}(\Pi(\hat{\Gamma}_{c})-\Pi(E_{c}\Gamma))=O_{p}\left(N^{-1/2}\right)$.
Since $T=O(1)$, it must be that $C=O(1)$ and hence
$\lVert\hat{A}-A\rVert_{\mathrm{op}}=O_{p}\left(N^{-1/2}\right)$ as well.

Next, let $\lambda_{j}$ be the $j$th smallest eigenvalue of the population APM
$A$, and denote a generic eigendecomposition of $A$ as follows:

$A=\sum_{j=1}^{T}\lambda_{j}\Pi(u_{j}),$ (A.19)

where $u_{j}$ is an eigenvector corresponding to the $j$th eigenvalue
$\lambda_{j}$, and $\Pi(u_{j})$ is the projection onto the span of
$u_{j}$.343434Even though the matrix of eigenvectors whose columns are $u_{j}$
is required to be orthonormal, its columns are only unique up to signs.
Furthermore, if $A$ has repeated eigenvalues, then any of the eigenvectors
corresponding to those repeated eigenvalues are interchangeable. We note that
all of $A$’s eigenvalues must be non-
negative.353535$E_{c}-\Pi(E_{c}\Gamma)=E_{c}[I-\Pi(E_{c}\Gamma)]E_{c}$, so
since $E_{c}$ and $I-\Pi(E_{c}\Gamma)$ are both projection matrices, their
product must be positive semidefinite. The sum of positive semidefinite
matrices must also be positive semidefinite, so then $A$ must be positive
semidefinite as well by (2.5).

Given this notation, we now apply Theorem B.1 with $M=A$, $\hat{M}=\hat{A}$ as
defined in Section 3, $s=0$, and $r=r$. We also note that
$\lambda_{j}=\lambda_{j}(A)$ using the notation defined in Appendix B. Since
Theorem 4.2 implies that

$0=\lambda_{1}=\dotsc=\lambda_{r}<\lambda_{r+1},$ (A.20)

we have that $\Delta(A)=4^{-1}\lambda_{r+1}>0$, satisfying (B.4). Further,
Theorem 4.2 implies that

$\Pi(\Gamma)=\Pi(U_{1:r}(A)),$ (A.21)

and we define $\hat{\Gamma}$ in Section 3 such that

$\Pi(\hat{\Gamma})=\Pi(U_{1:r}(\hat{A})).$ (A.22)

Then so long as

$\lVert\hat{A}-A\rVert_{\mathrm{op}}\leq 4^{-1}\lambda_{r+1},$ (A.23)

Theorem B.1 implies that

$\displaystyle\Bigg{\lVert}\Pi(\hat{\Gamma})-\Pi(\Gamma)-\sum_{j=1}^{r}\sum_{k=r+1}^{T}\lambda_{k}^{-1}\left[\Pi(u_{j})(\hat{A}-A)\Pi(u_{k})+\Pi(u_{k})(\hat{A}-A)\Pi(u_{j})\right]\Bigg{\rVert}_{\mathrm{op}}$
(A.24)
$\displaystyle\leq\frac{32}{\pi\lambda_{r+1}^{2}}\lVert\hat{A}-A\rVert_{\mathrm{op}}^{2}.$

The fact that $\lVert\hat{A}-A\rVert_{\mathrm{op}}=O_{p}\left(N^{-1/2}\right)$
implies (A.23) holds with probability approaching one as $N\rightarrow\infty$,
so by (A.24) and the definitions of $\hat{A}$ and $A$,

$\displaystyle\Pi(\hat{\Gamma})-\Pi(\Gamma)$ (A.25)
$\displaystyle=\sum_{j=r+1}^{T}\frac{1}{\lambda_{j}}\sum_{k=1}^{r}\bigg{\\{}\Pi(u_{j})(\hat{A}-A)\Pi(u_{k})+\Pi(u_{k})(\hat{A}-A)\Pi(u_{j})\bigg{\\}}+o_{p}\left(N^{-1/2}\right).$
(A.26)
$\displaystyle=\left(\sum_{j=r+1}^{T}\frac{1}{\lambda_{j}}\Pi(u_{j})\right)(\hat{A}-A)\left(\sum_{k=1}^{r}\Pi(u_{k})\right)$
(A.27)
$\displaystyle\phantom{=}+\left(\sum_{k=1}^{r}\Pi(u_{k})\right)(\hat{A}-A)\left(\sum_{j=r+1}^{T}\frac{1}{\lambda_{j}}\Pi(u_{j})\right)+o_{p}\left(N^{-1/2}\right).$
(A.28)

Next, since $\lambda_{1}=\dotsc=\lambda_{r}=0$, it must be that

$\sum_{j=r+1}^{T}\frac{1}{\lambda_{j}}\Pi(u_{j})=\sum_{j=1}^{T}\frac{\mathbbm{1}\\{j>r\\}}{\lambda_{j}}\Pi(u_{j})=A^{+},$
(A.29)

and since $\mathrm{rank}(\Gamma)=r$ and
$\mathrm{null}(A)=\mathrm{col}(\Gamma)$,

$\sum_{k=1}^{r}\Pi(u_{k})=\Pi(\Gamma).$ (A.30)

Substituting these two simplifications back into (A.28) and applying the
definitions of $\hat{A}$ and $A$, we have

$\displaystyle\Pi(\hat{\Gamma})-\Pi(\Gamma)$ (A.31)
$\displaystyle=\sum_{c=1}^{C}\left\\{A^{+}[\Pi(\hat{\Gamma}_{c})-\Pi(E_{c}\Gamma)]\Pi(\Gamma)+\Pi(\Gamma)[\Pi(\hat{\Gamma}_{c})-\Pi(E_{c}\Gamma)]A^{+}\right\\}+o_{p}\left(N^{-1/2}\right).$

Next, we note that the so-called mixed product property of the Kronecker
product implies that, for any matrices $A\in\operatorname{\mathbb{R}}^{d\times
m}$, $B\in\operatorname{\mathbb{R}}^{m\times p}$, and
$C\in\operatorname{\mathbb{R}}^{p\times q}$,

$\mathrm{vec}(ABC)=(C^{\prime}\otimes A)\mathrm{vec}(B).$ (A.32)

Applying this fact to the vectorization of $\sqrt{N}$ times both sides of
(A.31), since any projection matrix $\Pi(\cdot)$ is symmetric, we can write
the vectorization of our scaled and centered statistic
$\sqrt{N}(\Pi(\hat{\Gamma})-\Pi(\Gamma))$ as follows:

$\displaystyle\sqrt{N}\mathrm{vec}(\Pi(\hat{\Gamma})-\Pi(\Gamma))$ (A.33)
$\displaystyle=\left(A^{+}\otimes\Pi(\Gamma)+\Pi(\Gamma)\otimes
A^{+}\right)\left(\sum_{c=1}^{C}\sqrt{N}\mathrm{vec}\left(\Pi(\hat{\Gamma}_{c^{*}})-\Pi(E_{c^{*}}\Gamma)\right)\right)+o_{p}\left(1\right).$

Combining (A.33) and (4.4) yields (4.5), as required.

By the triangle inequality and the product property of Kronecker product
spectra (P. Lancaster Farahat, 1972),

$\left\lVert H\right\rVert_{\mathrm{op}}\leq\left\lVert
A^{+}\otimes\Pi(\Gamma)\right\rVert_{\mathrm{op}}+\left\lVert\Pi(\Gamma)\otimes
A^{+}\right\rVert_{\mathrm{op}}=2\left\lVert
A^{+}\right\rVert_{\mathrm{op}}\left\lVert\Pi(\Gamma)\right\rVert_{\mathrm{op}}=\frac{2}{\lambda_{r+1}}<\infty.$

The zero mean and bounded squared norm properties of
$H\sum_{c=1}^{C}\phi_{c}(C_{i},Y_{i})$ follow from Assumption 4.3 and the fact
that the operator norm of $H$ is bounded.

###  A.5 Proof of Theorem 4.5

First, we show that with probability approaching one as $N\rightarrow\infty$,
$\hat{\Gamma}^{\prime}E_{c^{*}}\hat{\Gamma}$ is invertible. To do so, we note
that by Theorem 4.4, $\Pi(\hat{\Gamma})-\Pi(\Gamma)=o_{p}\left(1\right)$, so
there must be a random, $r\times r$ invertible basis matrix $\hat{Q}$ such
that $\hat{\Gamma}\hat{Q}^{\prime}=\Gamma+o_{p}\left(1\right)$; if not, then
for all random, $r\times r$ invertible basis matrices $\hat{Q}$,
$\Pi(\hat{\Gamma})=\Pi(\hat{\Gamma}\hat{Q}^{\prime})\neq\Pi(\Gamma)+o_{p}\left(1\right)$,
which would be a violation of Theorem 4.4. Since Lemma 4.1 implies that that
$\Gamma^{\prime}E_{c^{*}}\Gamma$ is invertible and
$\hat{\Gamma}^{\prime}E_{c^{*}}\hat{\Gamma}=(\hat{Q})^{-1}\Gamma
E_{c^{*}}\Gamma(\hat{Q}^{\prime})^{-1}+o_{p}\left(1\right)$, it must be that
$\hat{\Gamma}^{\prime}E_{c^{*}}\hat{\Gamma}$ is invertible with probability
approaching one as $N\rightarrow\infty$. As such, in the remainder of the
proof below, we condition on the event that
$\hat{\Gamma}^{\prime}E_{c^{*}}\hat{\Gamma}$ is invertible.

Next, let
$B_{c^{*}}\coloneqq\Gamma(\Gamma^{\prime}E_{c^{*}}\Gamma)^{-1}\Gamma^{\prime}E_{c^{*}}$
and
$\hat{B}_{c^{*}}\coloneqq\hat{\Gamma}(\hat{\Gamma}^{\prime}E_{c^{*}}\hat{\Gamma})^{-1}\hat{\Gamma}^{\prime}E_{c^{*}}$,
and note that

$\displaystyle\hat{\mu}_{c^{*}}-\mu_{c^{*}}$
$\displaystyle=\hat{B}_{c^{*}}\hat{\operatorname{\mathbb{E}}}_{N}[E_{c^{*}}Y_{i}~{}|~{}C_{i}=c^{*}]-B_{c^{*}}\operatorname{\mathbb{E}}[E_{c^{*}}Y_{i}~{}|~{}C_{i}=c^{*}]$
(A.34)
$\displaystyle=(\hat{B}_{c^{*}}-B_{c^{*}})E_{c^{*}}\mu_{c^{*}}+B_{c^{*}}(\hat{\operatorname{\mathbb{E}}}_{N}-\operatorname{\mathbb{E}})[E_{c^{*}}Y_{i}~{}|~{}C_{i}=c^{*}]$
(A.35)
$\displaystyle\phantom{=}+(\hat{B}_{c^{*}}-B_{c^{*}})(\hat{\operatorname{\mathbb{E}}}_{N}-\operatorname{\mathbb{E}})[E_{c^{*}}Y_{i}~{}|~{}C_{i}=c^{*}].$
(A.36)

To approximate the expression $\hat{B}_{c^{*}}-B_{c^{*}}$, we will rely on the
following lemma, which we prove in Appendix C.1:

######  Lemma A.1.

For any two $T\times r$ matrices $\hat{\Gamma}$ and $\Gamma$ such that
$E_{c^{*}}\hat{\Gamma}$ and $E_{c^{*}}\Gamma$ are full-rank, let
$\hat{\Delta}_{\Gamma}\coloneqq\Pi(\hat{\Gamma})-\Pi(\Gamma)$. If
$\lVert\hat{\Delta}_{\Gamma}\rVert_{\mathrm{op}}\leq 1/\left(2\lVert
R_{c^{*}}\rVert_{\mathrm{op}}\right)$, then the following first-order
approximation holds:

$\displaystyle\left\lVert(\hat{B}_{c^{*}}-B_{c^{*}})-\left(2E_{c^{*}}\hat{\Delta}_{\Gamma}R_{c^{*}}^{\prime}+(I+R_{c^{*}})\hat{\Delta}_{\Gamma}(I-R_{c^{*}})^{\prime}\right)E_{c^{*}}\right\rVert_{\mathrm{op}}\leq
b\lVert\hat{\Delta}_{\Gamma}\rVert_{\mathrm{op}}^{2},$ (A.37)

where $b\coloneqq 12\left(1+\left\lVert
R_{c^{*}}\right\rVert_{\mathrm{op}}\right)^{3}<\infty$.

Since
$\lVert\hat{\Delta}_{\Gamma}\rVert_{\mathrm{op}}=O_{p}\left(N^{-1/2}\right)$,
the probability that $\lVert\hat{\Delta}_{\Gamma}\rVert_{\mathrm{op}}\leq
1/\left(2\lVert R_{c^{*}}\rVert_{\mathrm{op}}\right)$ converges to one. For
the rest of the proof we also condition on this event so that Lemma A.1 can be
applied. Letting

$\hat{F}_{c^{*}}\coloneqq\left[2E_{c^{*}}\hat{\Delta}_{\Gamma}R_{c^{*}}^{\prime}+(I+R_{c^{*}})\hat{\Delta}_{\Gamma}(I-R_{c^{*}})^{\prime}\right]E_{c^{*}}.$
(A.38)

By Lemma A.1,

$\displaystyle\lVert\hat{B}_{c^{*}}-B_{c^{*}}\rVert_{\mathrm{op}}\leq\left[2\left\lVert
R_{c^{*}}\right\rVert_{\mathrm{op}}+(1+\left\lVert
R_{c^{*}}\right\rVert_{\mathrm{op}})^{2}\right]\lVert\hat{\Delta}_{\Gamma}\rVert_{\mathrm{op}}+b\lVert\hat{\Delta}_{\Gamma}\rVert_{\mathrm{op}}^{2}\leq
b\left(\lVert\hat{\Delta}_{\Gamma}\rVert_{\mathrm{op}}+\lVert\hat{\Delta}_{\Gamma}\rVert_{\mathrm{op}}^{2}\right).$
(A.39)

Using the fact that $\hat{F}_{c^{*}}E_{c^{*}}=\hat{F}_{c^{*}}$, we can further
expand (A.36) as follows:

$\displaystyle\hat{\mu}_{c^{*}}-\mu_{c^{*}}$ (A.40)
$\displaystyle=\hat{F}_{c^{*}}\mu_{c^{*}}+B_{c^{*}}(\hat{\operatorname{\mathbb{E}}}_{N}-\operatorname{\mathbb{E}})[E_{c^{*}}Y_{i}~{}|~{}C_{i}=c^{*}]$
$\displaystyle\phantom{=}+(\hat{B}_{c^{*}}-B_{c^{*}}-\hat{F}_{c^{*}})E_{c^{*}}\mu_{c^{*}}+(\hat{B}_{c^{*}}-B_{c^{*}})(\hat{\operatorname{\mathbb{E}}}_{N}-\operatorname{\mathbb{E}})[E_{c^{*}}Y_{i}~{}|~{}C_{i}=c^{*}]$

Rearranging the display above and taking the norm, we have that, again by
Lemma A.1,

$\displaystyle\left\lVert\left[\hat{\mu}_{c^{*}}-\mu_{c^{*}}\right]-\hat{F}_{c^{*}}\mu_{c^{*}}-B_{c^{*}}(\hat{\operatorname{\mathbb{E}}}_{N}-\operatorname{\mathbb{E}})[E_{c^{*}}Y_{i}~{}|~{}C_{i}=c^{*}]\right\rVert_{2}$
(A.41)
$\displaystyle\leq\lVert\hat{B}_{c^{*}}-B_{c^{*}}-\hat{F}_{c^{*}}\rVert_{\mathrm{op}}\left\lVert
E_{c^{*}}\mu_{c^{*}}\right\rVert_{2}$
$\displaystyle\phantom{=}+\lVert\hat{B}_{c^{*}}-B_{c^{*}}\rVert_{\mathrm{op}}\lVert(\hat{\operatorname{\mathbb{E}}}_{N}-\operatorname{\mathbb{E}})[E_{c^{*}}Y_{i}~{}|~{}C_{i}=c^{*}]\rVert_{2}$
$\displaystyle\leq
b\lVert\hat{\Delta}_{\Gamma}\rVert_{\mathrm{op}}^{2}\left\lVert
E_{c^{*}}\mu_{c^{*}}\right\rVert_{2}+b\left(\lVert\hat{\Delta}_{\Gamma}\rVert_{\mathrm{op}}+\lVert\hat{\Delta}_{\Gamma}\rVert_{\mathrm{op}}^{2}\right)\lVert(\hat{\operatorname{\mathbb{E}}}_{N}-\operatorname{\mathbb{E}})[E_{c^{*}}Y_{i}~{}|~{}C_{i}=c^{*}]\rVert_{2}.$

Under Assumption 4.3, Theorem 4.4 implies that
$\lVert\hat{\Delta}_{\Gamma}\rVert_{\mathrm{op}}=O_{p}\left(N^{-1/2}\right)$,
and under Assumption 4.4,

$\lVert(\hat{\operatorname{\mathbb{E}}}_{N}-\operatorname{\mathbb{E}})[E_{c^{*}}Y_{i}~{}|~{}C_{i}=c^{*}]\rVert_{2}\leq\lVert(\hat{\operatorname{\mathbb{E}}}_{N}-\operatorname{\mathbb{E}})[Y_{i}^{*}~{}|~{}C_{i}=c^{*}]\rVert_{2}=O_{p}\left(N^{-1/2}\right).$
(A.42)

Therefore, by (A.41),

$\displaystyle\sqrt{N}(\hat{\mu}_{c^{*}}-\mu_{c^{*}})$
$\displaystyle=\sqrt{N}\hat{F}_{c^{*}}\mu_{c^{*}}+\sqrt{N}B_{c^{*}}(\hat{\operatorname{\mathbb{E}}}_{N}-\operatorname{\mathbb{E}})[E_{c^{*}}Y_{i}~{}|~{}C_{i}=c^{*}]+o_{p}\left(1\right).$
(A.43)

Applying (A.32) to the first term in the expansion (A.43), we have that

$\displaystyle\sqrt{N}\hat{F}_{c^{*}}\mu_{c^{*}}$ (A.44)
$\displaystyle=\left((2\mu_{c^{*}}^{\prime}E_{c^{*}}R_{c^{*}})\otimes
E_{c^{*}}+(\mu_{c^{*}}^{\prime}E_{c^{*}}(I-R_{c^{*}}))\otimes(I+R_{c^{*}})\right)\cdot\sqrt{N}\mathrm{vec}(\Pi(\hat{\Gamma})-\Pi(\Gamma))$
(A.45)
$\displaystyle=M_{c^{*}}\cdot\sqrt{N}\mathrm{vec}(\Pi(\hat{\Gamma})-\Pi(\Gamma))$
(A.46)
$\displaystyle=\sqrt{N}\hat{\operatorname{\mathbb{E}}}_{N}\left[M_{c^{*}}H\sum_{c=1}^{C}\phi_{c}(C_{i},Y_{i})\right]+o_{p}\left(1\right).$
(A.47)

The following lemma, which we prove in Appendix C.2, will help us complete the
proof:

######  Lemma A.2.

Let $V_{1},V_{2},\dotsc$ be a sequence of i.i.d. random variables such that
$\operatorname{\mathbb{E}}[\lVert V_{1}\rVert_{2}^{2}]<\infty$, and let
$D_{i}$ be a binary random variable such that
$\operatorname{\mathbb{P}}(D_{i}=1)>0$. Then

$\displaystyle\sqrt{N}\hat{\operatorname{\mathbb{E}}}_{N}[V_{i}~{}|~{}D_{i}=1]-\operatorname{\mathbb{E}}[V_{i}~{}|~{}D_{i}=1]=\sqrt{N}\hat{\operatorname{\mathbb{E}}}_{N}\left[\frac{D_{i}}{\operatorname{\mathbb{P}}(D_{i}=1)}\left(V_{i}-\operatorname{\mathbb{E}}[V_{i}~{}|~{}D_{i}=1]\right)\right]+o_{p}\left(1\right).$
(A.48)

By Lemma A.2, we have that

$(\hat{\operatorname{\mathbb{E}}}_{N}-\operatorname{\mathbb{E}})[E_{c^{*}}Y_{i}~{}|~{}C_{i}=c^{*}]=\hat{\operatorname{\mathbb{E}}}_{N}\left[\frac{\mathbbm{1}\\{C_{i}=c^{*}\\}}{\operatorname{\mathbb{P}}(C_{i}=c^{*})}\left(E_{c^{*}}Y_{i}-E_{c^{*}}\mu_{c^{*}}\right)\right].$
(A.49)

Substituting the expressions in the two displays above back into (A.43) yields
the asymptotically linear expansion in (4.7). Further, under (4.4) and (4.4),
$\operatorname{\mathbb{E}}[\psi_{c^{*}}(C_{i},Y_{i})]=\operatorname{\mathbf{0}}_{T}$
and

$\operatorname{\mathbb{E}}[\left\lVert\psi_{c^{*}i}\right\rVert_{2}^{2}]\leq\left\lVert
M_{c^{*}}\right\rVert_{\mathrm{op}}^{2}\left\lVert
H\right\rVert_{\mathrm{op}}^{2}\sum_{c=1}^{C}\operatorname{\mathbb{E}}[\left\lVert\phi_{c}(C_{i},Y_{i})\right\rVert_{2}^{2}]+2\operatorname{\mathbb{E}}\left[\lVert
Y_{i}^{*}\rVert_{2}^{2}~{}\middle|~{}C_{i}=c^{*}\right]<\infty.$ (A.50)

##  Appendix B An Exact, First-Order Expansion of the Eigenspace Operator on
Symmetric Matrices

In this section, we derive an exact, first-order expansion of the operator
mapping a symmetric matrix into the projection matrix onto the space spanned
by some subset of its eigenvectors. Our expansion is based on Kato’s integral,
which characterizes the projection matrix onto the space spanned by some
subset of a matrix’s eigenvectors as a contour integral of that matrix’s
resolvent Kato (1949). We characterize the approximation error in our
expansion up to exact constants, unlike the asymptotic expansions given in
e.g. Kato (1980) and Jg. Sun (1991) (which is applied in Simons (2023) to
construct hypothesis tests concerning eigenspaces). As such, our result below
may be of independent interest.

To describe our expansion, we first introduce some notation. Let
$\mathbb{S}^{d}$ denote the set of real-valued, $d$-dimensional symmetric
matrices, and for any $M\in\mathbb{S}^{d}$, let $\lambda_{j}(M)$ be the $j$th
eigenvalue of $A$, where

$\lambda_{0}(M)\coloneqq-\infty<\lambda_{1}(M)\leq\cdots\leq\lambda_{d}(M)<\lambda_{d+1}(M)\coloneqq\infty,$
(B.1)

and we denote a generic eigendecomposition of $M$ as follows:

$M=U(M)\Lambda(M)U(M)^{\prime}=\sum_{j=1}^{d}\lambda_{j}(M)\Pi(u_{j}(M)),$
(B.2)

where $U(M)$ is any matrix with $j$th column $u_{j}(M)$ being an eigenvector
corresponding to the $j$th eigenvalue $\lambda_{j}(M)$, and $\Pi(u_{j}(M))$ is
the projection onto the span of $u_{j}(M)$.373737Even though $U(M)$ is
required to be orthonormal, the columns of $U(M)$ are only unique up to signs
and permutations any eigenvectors corresponding to repeated eigenvalues. For
notational convenience, for any integers $1\leq j,k\leq d$, we let

$U_{j:k}(M)\coloneqq\begin{bmatrix}u_{j}(M)&\cdots&u_{k}(M)\end{bmatrix}$
(B.3)

denote the matrix whose columns are the eigenvectors corresponding to the
$j$th through $k$th eigenvalues.

Armed with this notation, we state our main result. In particular, we
construct an exact bound on the error incurred by a first-order expansion of
the difference between the projection matrices onto eigenspaces of two
matrices $M$ and $\hat{M}$ in terms of the magnitude of the difference between
$M$ and $\hat{M}$:

######  Theorem B.1.

Consider any integers $s,r$ such that $1\leq s+1\leq s+r\leq d$ and any real-
valued, $d$-dimensional symmetric matrix $M$ satisfying the following eigen-
gap condition:

$\Delta(M)>0,\quad\Delta(M)\coloneqq
4^{-1}\min\left\\{\lambda_{s+1}(M)-\lambda_{s}(M),\lambda_{s+r+1}(M)-\lambda_{s+r}(M)\right\\},$
(B.4)

and define the following neighborhood of $M$:

$\mathcal{B}(M)\coloneqq\left\\{\hat{M}\in\mathbb{S}^{d}~{}\operatorname{:}~{}\lVert\hat{M}-M\rVert_{\mathrm{op}}\leq\Delta(M)\right\\}.$
(B.5)

Then for any $\hat{M}\in\mathcal{B}(M)$, the following first-order
approximation holds:

$\displaystyle\Bigg{\lVert}\left[\Pi(U_{(s+1):(s+r)}(\hat{M}))-\Pi(U_{(s+1):(s+r)}(M))\right]$
(B.6)
$\displaystyle\phantom{\big{\lVert}}-\sum_{j=s+1}^{s+r}\sum_{k\not\in[s+1,s+r]}\frac{1}{\lambda_{k}(M)-\lambda_{j}(M)}\Big{[}\Pi(u_{j}(M))(\hat{M}-M)\Pi(u_{k}(M))$
$\displaystyle\phantom{\big{\lVert}-\sum_{j=s+1}^{s+r}\sum_{k\not\in[s+1,s+r]}\frac{1}{\lambda_{k}(M)-\lambda_{j}(M)}\Big{[}}+\Pi(u_{k}(M))(\hat{M}-M)\Pi(u_{j}(M))\Big{]}\Bigg{\rVert}_{\text{op}}$
$\displaystyle\leq\frac{2}{\pi\Delta(M)^{2}}\lVert\hat{M}-M\rVert_{\text{op}}^{2}.$

###### Proof.

Consider any $\hat{M}\in\mathcal{B}(M)$. By Weyl’s inequality (see e.g.
Theorem 4.5.3 in Vershynin (2018)),

$\max_{j}\lvert\lambda_{j}(M)-\lambda_{j}(\hat{M})\rvert\leq\lVert
M-\hat{M}\rVert_{\text{op}}\leq\Delta(M).$ (B.7)

Let

$a(M)\coloneqq\begin{cases}\frac{\lambda_{s+1}(M)+\lambda_{s}(M)}{2},&s>0\\\
\lambda_{1}(M)-2\Delta(M),&s=0\end{cases},\quad
b(M)\coloneqq\begin{cases}\frac{\lambda_{s+r+1}(M)+\lambda_{s+r}(M)}{2},&s+r<d\\\
\lambda_{d}(M)+2\Delta(M),&s+r=d\end{cases}.$ (B.8)