The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
100,948 characters
Low-Rank Approximations of Nonseparable Panel Models
\begin{abstract}
We provide estimation methods for nonseparable panel models based on low-rank factor structure approximations. The factor structures are estimated by matrix-completion methods to deal with the computational challenges of principal component analysis in the presence of missing data. We show that the resulting estimators are consistent in large panels, but suffer from approximation and shrinkage biases. We correct these biases using matching and difference-in-differences approaches. Numerical examples and an empirical application to the effect of election day registration on voter turnout in the U.S. illustrate the properties and usefulness of our methods.
\keywords{Nonseparable Panel, Low-Rank Approximations, Matrix Completion, Debias, Two-Way Matching, Election Day Registration}
\end{abstract}
\section{Introduction}
\setcounter{equation}{0}
Nonseparable models are useful to capture multidimensional unobserved heterogeneity, which is an important feature of economic data. The presence of this heterogeneity makes the effect of covariates on the outcome of interest different for each unit due to factors that are unobservable or unavailable to the researcher. In the absence of further restrictions, a different data generating process essentially operates for each unit, which creates identification and estimation challenges. One way to deal with these challenges is the use of panel data, where each unit is observed on multiple occasions.
In this paper, we develop an approach to estimate nonseparable models from panel data based on homogeneity restrictions and low-rank factor approximations. Whilst homogeneity restrictions have been used previously in this context, the application of low-rank factor approximations is more novel.
The nonseparable model that we consider includes observed discrete covariates or treatments, multidimensional unobserved individual and time effects, and idiosyncratic errors. We construct the effects of interest as averages or quantiles of potential outcomes constructed from the model by exogenously manipulating the value of the treatments. These effects are generally not identified from the observed data because the treatment assignment is usually determined by the unobserved individual and time effects. Following the previous panel literature, we impose cross-section and time-series homogeneity restrictions to identify the effects of interest, see, e.g. \citet{Chamberlain82}, \citet{Manski1987}, \citet{Honore1992}, \citet{evdokimov2010identification}, \citet{GrahamPowell2012}, \citet{HoderleinWhite2012} and \citet{CFHN13}.
The estimation of the nonseparable model is challenging due to the presence of the multidimensional unobserved individual and time effects. We cannot just exclude these effects because they are endogenous, i.e., related to the treatments. We deal with this problem by approximating their effect with a low-rank factor structure. This approach can be interpreted as a series or sieve approximation on the unobservables. We characterize the error of this approximation in terms of the functional singular value decomposition of the expectation of the outcome conditional on the treatment and unobserved effects. For smooth conditional expectation functions, the mean squared error of the approximation error vanishes with the rank of the factor structure at a polynomial rate.
We develop an estimator of the low-rank factor approximation in the case where the covariate of interest is binary. This is an empirically relevant case as it covers the treatment effect model for panel data. We also show how to extend the model to include additive controls and fixed effects. Here, we rely on the analogy between the estimation of treatment effects and the matrix completion problem previously noted by \citet{athey17} and \citet{amjad18}. Thus, given that the principal components program is combinatorially hard in the presence of missing data, we consider the convex relaxation of this program that replaces a constraint in the rank of a matrix by a constraint in its nuclear norm, following \citet{srebro03} and \citet{fazel03}. The resulting estimator is the matrix-completion estimator.
The main theoretical result of the paper is to show that the matrix-completion estimator is consistent under asymptotic sequences where the two dimensions of the panel grow to infinity at the same rate. This result does not follow from the existing matrix completion literature that assumes that the matrix to complete has low-rank. In our case, the underlying matrix of interest can have full rank, but we impose
appropriate smoothness assumptions on the data generating process that guarantee that the singular values of the matrix
form a rapidly decreasing sequence. This allows a low-rank approximation, and it also implies a bound on the nuclear norm of the
matrix. Our consistency proof for the matrix completion estimator therefore crucially relies on the bound of the nuclear norm, but does
not impose any low-rank conditions. Our proof strategy also avoids the high-level \emph{restricted strong convexity} assumption
(see e.g.\ \citet{negahban2012restricted}). We instead provide interpretable conditions on the underlying process of the observable and unobservable variables directly.
The matrix-completion estimator is consistent, but can be biased in small samples. This bias comes from two different sources: approximation bias due to the low-rank factor structure approximation and shrinkage bias due to the nuclear norm regularization of the principal component analysis program \citet{cai10,ma11,bai19joe}. We propose matching approaches to debias the estimator. For each treatment level, the simplest approach consists of finding the observation in the other treatment level that is the closest in terms of the estimated factor structure. We also propose a two-way matching procedure that combines matching with a differences-in-differences approach. The two-way procedure is related to several recent proposals such as the matching approach of \citet{imai2019use} to estimate causal effects from panel data and the blind regression of \citet{li17b} for matrix completion. The difference with these proposals is in the information used to match the observations. \citet{imai2019use} use the treatment variable and \citet{li17b} the outcome, whereas we use the estimated factor structure. In this sense, the estimation of the factor structure can be seen as a preliminary de-noising step of the data \citet{chatterjee2015}. \citet{amjad18} proposed a similar debiasing procedure based on the estimated factor structure, but they rely on synthetic control methods instead of matching. In contemporaneous and independent work, \citet{chlz20} have developed an alternative rotation-debiasing method that can be applied to make inference on heterogenous treatment effects in low-rank models. This method consists of the application of iterative least squares to the left and right singular vectors of the matrix-completion estimator.
We illustrate our methods with an empirical application to the effect of election day registration (EDR) on voter turnout and numerical simulations. We estimate average and quantile effects using a state-level panel dataset on the 24 U.S. presidential elections between 1920 and 2012 collected by \citet{xu17}. We find that, after controlling for possible non-random adoption, EDR has a positive effect, especially at the bottom of the voter turnout distribution. Our methods uncover stronger effects than standard difference-in-differences methods that rely on restrictive parallel trend assumptions. The simulation results show that our theoretical results provide a good representation of the behavior of the estimators in small samples.
The rest of the paper is organized as follows. Section \ref{sec:model} describes the model and effects of interest. Section \ref{sec:est} introduces the low-rank factor approximation and derives the properties of its matrix-completion estimator. The matching methods to debias the matrix-completion estimator are discussed in Section \ref{sec:debias}. Section \ref{sec:examples} reports the results of the numerical examples.
All the proofs of the theoretical results are gathered in the Appendix.
\section{Model and Effects of Interest}\label{sec:model}
\setcounter{equation}{0}
Throughout this paper we consider the following nonseparable and nonparametric panel data model:
\begin{assumption}[Model]\label{ass:model}
\begin{align}
Y_{it} &= g(\boldsymbol{X}_{it}, \boldsymbol{A}_i, \boldsymbol{B}_t, \boldsymbol{U}_{it} ) ,
\qquad
i \in \mathbb{N} = \{1,\ldots,N\} ,
\;
t \in \mathbb{T} = \{1,\ldots,T\} ,
\label{model}
\end{align}
where $i$ and $t$ index individual units and time periods, respectively; $Y_{it}$ is an observed outcome or response variable with support $\mathbb{Y} \subseteq \mathbb{R}$; $g$ is an unknown function; $\boldsymbol{X}_{it}$ is a vector of observed covariates or treatments with finite support $\mathbb{X}$; $\boldsymbol{A}_i$
and $\boldsymbol{B}_t$ are vectors of individual and time unobserved effects, possibly correlated with $\boldsymbol{X}_{it}$, with supports $\mathbb{A} \subseteq \mathbb{R}^{d_a}$ and $\mathbb{B} \subseteq \mathbb{R}^{d_b}$, respectively; and $\boldsymbol{U}_{it}$ is a vector of unobserved error terms of unspecified dimension,
for which we assume that
\begin{align}
\boldsymbol{U}_{it} \overset{d}{=} \boldsymbol{U}_{js} \mid \boldsymbol{X}^{NT}, \boldsymbol{A}^N, \boldsymbol{B}^T,
\qquad
\text{for all } i,j \in \mathbb{N}, \; t,s \in \mathbb{T},
\label{ass:Main}
\end{align}
and
\begin{equation}\label{ass:Main2}
\boldsymbol{U}_{it} \perp\!\!\!\perp (\boldsymbol{X}^{NT}, \boldsymbol{A}^N, \boldsymbol{B}^T) \mid \boldsymbol{A}_i, \boldsymbol{B}_t, \qquad
\text{for all } i \in \mathbb{N}, \; t \in \mathbb{T},
\end{equation}
where $\boldsymbol{X}^{NT} = \{ \boldsymbol{X}_{it}: i \in \mathbb{N}, t \in \mathbb{T}\}$, $\boldsymbol{A}^N = \{ \boldsymbol{A}_i: i \in \mathbb{N}\}$, $\boldsymbol{B}^T = \{ \boldsymbol{B}_t: t \in \mathbb{T}\}$, and $\perp\!\!\!\perp$ denotes stochastic independence. We also assume that, for all $i \in \mathbb{N}, \; t \in \mathbb{T}$, the support of $(\boldsymbol{X}_{it}, \boldsymbol{A}_i, \boldsymbol{B}_t)$ is
equal to the Cartesian product $\mathbb{X} \times \mathbb{A} \times \mathbb{B}$, and that $\mathrm{E} \, Y_{it}^2 < \infty$.
\end{assumption}
This model can be motivated from a purely statistical perspective as a latent variable model using the Aldous-Hoover representation for exchangeable random matrices, e.g. \citet{xu14}, \citet{chatterjee2015}, \citet{or15}, and \citet{li17}.\footnote{In the Aldous-Hoover representation, $\boldsymbol{A}_i$, $\boldsymbol{B}_t$ and $\boldsymbol{U}_{it}$ are independent uniform random variables.} We motivate it instead as a structural model where the unobserved effects $\boldsymbol{A}_i$ and $\boldsymbol{B}_t$ are associated with individual heterogeneity and aggregate shocks, respectively. Additional exogenous covariates can be incorporated in the usual way by carrying out the analysis conditional on them.
We focus on discrete covariates but, from a theoretical perspective, the extension to continuous covariates is
straightforward by using appropriate smoothing methods --- it is, however, not clear to us whether that extension would be practically useful with realistic sample
sizes. We therefore think that it would complicate our presentation without much benefit.
The main restriction imposed by Assumption \ref{ass:model} is the unit and time homogeneity in \eqref{ass:Main}. A sufficient condition for unit homogeneity is that the observations are identically distributed across $i$, which is a common sampling assumption for panel data. Time homogeneity has also been commonly used in panel data models \citep{Chamberlain82,Manski1987,Honore1992,evdokimov2010identification,GrahamPowell2012,HoderleinWhite2012,CFHN13}. It implies that time is randomly assigned, conditional on covariates and unobserved effects. The additional restrictions in \eqref{ass:Main2} are exogeneity conditions on $(\boldsymbol{X}^{NT}, \boldsymbol{A}^N, \boldsymbol{B}^T)$ with respect to $\boldsymbol{U}_{it}$, conditional on $\boldsymbol{A}_i$ and $\boldsymbol{B}_t$. The most substantive is the exogeneity of $\boldsymbol{X}_{it}$. Given \eqref{ass:Main}, this is a mild condition as time homogeneity already imposes that any relationship between $\boldsymbol{U}_{it}$ and $\boldsymbol{X}_{it}$ can only be unit and time-invariant. Taken together, \eqref{ass:Main} and \eqref{ass:Main2} impose that
\begin{align}\label{ass:Main3}
\boldsymbol{U}_{it} \mid \boldsymbol{A}_i, \boldsymbol{B}_t \overset{d}{=} \boldsymbol{U}_{js} \mid \boldsymbol{A}_j, \boldsymbol{B}_s,
\qquad
\text{for all } i,j \in \mathbb{N}, \; t,s \in \mathbb{T}.
\end{align}
The product support condition guarantees overlap in the support of the unobserved effects for all values of the treatments. This condition is similar to the overlap condition used in cross section treatment effect models under unconfoundedness or selection on observables. Thus, together with \eqref{ass:Main2}, it implies that $P_{it}(x) := \Pr \left( X_{it} = x \mid \boldsymbol{A}^N, \boldsymbol{B}^T \right)>0$, a.s., for all $i \in \mathbb{N}$, $t \in \mathbb{T}$ and $x \in \mathbb{X}$, where $P_{it}(x)$ is the analog of the propensity score in our setting. This condition is plausible in many applications. For example,
in our empirical application in Section~\ref{sec:Election}, $X_{it}= \mathbbm{1}\{ t \geq \tau_i\} $, where $\tau_i$
is the date of the law change in state $i$. In that case, if we consider $\tau_i$ to be a random variable with sufficiently large support conditional on the unobserved effects, then the condition $ P_{it}(x)>0$, a.s., is satisfied.
The model considered is similar to the static model in \citet{CFHN13}, but there are three important differences. First, the structural function $g$ has time effects as arguments and therefore allows the relationship between $Y_{it}$ and $\boldsymbol{X}_{it}$ to vary over time in an unrestricted fashion even under \eqref{ass:Main}. For example, it can include location and scale time effects. Second, \citet{CFHN13} impose that $Y_{it}$ and $\boldsymbol{X}_{it}$ are identically distributed across $i$, which is stronger than the unit homogeneity in \eqref{ass:Main3}. Thus, unit homogeneity does not restrict the treatment assignment process. Third, they analyze short panels, whereas we rely on large $T$ for identification. Our model also encompasses the nonseparable model with time effects in \citet{freyberger18}, where in our notation $Y_{it} = g_t(\boldsymbol{X}_{it}, \boldsymbol{A}_i^\mathrm{\scriptscriptstyle{T}}\boldsymbol{B}_t + \boldsymbol{U}_{it} )$.\footnote{Note that our model allows for $g$ to depend on $t$ because the dimension of $\boldsymbol{B}_t$ is unspecified.}
We provide more examples of models covered by Assumption \ref{ass:model} below.
The structural function $g$ is generally not identified, but can be used to construct interesting effects. Let $Y_{it}(\boldsymbol{x}) := g(\boldsymbol{x}, \boldsymbol{A}_i, \boldsymbol{B}_t , \boldsymbol{U}_{it}(\boldsymbol{x}) ) $ be the potential outcome for individual $i$ at time $t$ obtained by setting exogenously $\boldsymbol{X}_{it} =\boldsymbol{x} \in \mathbb{X}$, where
\begin{equation}\label{eq:rs}
\boldsymbol{U}_{it}(\boldsymbol{x}) \overset{d}{=} \boldsymbol{U}_{it} \mid \boldsymbol{A}^N, \boldsymbol{B}^T.
\end{equation}
Here we impose rank similarity as the distribution of $\boldsymbol{U}_{it}(\boldsymbol{x})$ conditional on $\boldsymbol{A}^N$ and $\boldsymbol{B}^T$ does not change with $\boldsymbol{x}$. The main effects of interest are the average structural functions (ASFs)
\begin{align}\label{eq:asf}
\mu_{t}(\boldsymbol{x}) := \frac 1 {N} \sum_{i=1}^N \mathrm{E} \left[ Y_{it}(\boldsymbol{x}) \mid \boldsymbol{A}^N, \boldsymbol{B}^T \right], \ \ \mu(\boldsymbol{x}) := \frac 1 {T} \sum_{t=1}^T \mu_{t}(\boldsymbol{x}),
\end{align}
and the conditional average structural functions (CASFs)
\begin{align}\label{eq:casf}
\mu_{t}(\boldsymbol{x} \mid \mathbb{X}_0) &:= \frac 1 {N_t(\mathbb{X}_0)} \sum_{i=1}^N \mathbbm{1}\{\boldsymbol{X}_{it} \in \mathbb{X}_0 \} \mathrm{E} \left[ Y_{it}(\boldsymbol{x}) \mid \boldsymbol{A}^N, \boldsymbol{B}^T \right], \ \ N_t(\mathbb{X}_0) = \sum_{i=1}^N \mathbbm{1}\{\boldsymbol{X}_{it} \in \mathbb{X}_0 \}, \notag \\
\mu(\boldsymbol{x} \mid \mathbb{X}_0) &:= \frac 1 {n(\mathbb{X}_0)} \sum_{t=1}^T N_t(\mathbb{X}_0) \mu_{t}(\boldsymbol{x} \mid \mathbb{X}_0), \ \ n(\mathbb{X}_0) = \sum_{t=1}^T N_t(\mathbb{X}_0),
\end{align}
where $\mathbb{X}_0 \subseteq \mathbb{X}$, provided that $n(\mathbb{X}_0) > 0$. The ASFs and CASFs correspond to averages of the potential outcome $Y_{it}(\boldsymbol{x})$ at a given time period or aggregated over the observed time periods. In both cases the average is over the cross sectional units in the observed sample or finite population. Infinite-population versions of the effects can be obtained by taking probability limits as $N \to \infty$. If $\boldsymbol{X}_{it}$ includes only a binary treatment, the ASFs and CASFs can be used to form treatment effects. For example, $\mu(1) - \mu(0)$ is the time-aggregated average treatment effect and $\mu_t(1 \,|\, \{1\}) - \mu_t(0 \,|\, \{1\})$ is the average treatment effect on the treated at time $t$. Distribution structural functions (DSFs) can be constructed analogously replacing $Y_{it}(\boldsymbol{x})$ by $\mathbbm{1} \{Y_{it}(\boldsymbol{x}) \leq y\}$ in \eqref{eq:asf} and \eqref{eq:casf} for $y \in \mathbb{Y}$. Quantile effects can then be formed by taking left-inverses of the DSFs and taking differences. For example, the $\tau$-quantile treatment effect at time $t$ is $q_{t,\tau}(1) - q_{t,\tau}(0)$, where
$$
q_{t,\tau}(\boldsymbol{x}) = \inf\left\{y \in \mathbb{Y}: \frac 1 {N} \sum_{i=1}^N \mathrm{E} \left[ \mathbbm{1}\{Y_{it}(\boldsymbol{x}) \leq y\} \mid \boldsymbol{A}^N, \boldsymbol{B}^T \right] \geq \tau \right\}.
$$
We provide some examples of data generating processes that satisfy Assumption \ref{ass:model}. The purpose is to show that Assumption \ref{ass:model} covers a great variety of models commonly used in empirical analysis. Our estimation methods are generic in that we do not need to specify the data generating process, besides of satisfying Assumption \ref{ass:model}. Of course, using more information about the data generating process would lead to more efficient estimators, but at the cost of robustness to model misspecification.
\begin{example}[Linear factor model]\label{ex:lfm} \textnormal{Consider the linear panel model with factor structure in the error terms:
$$
Y_{it}(\boldsymbol{x}) = \boldsymbol{x}^\mathrm{\scriptscriptstyle{T}} \boldsymbol{\beta} + \boldsymbol{\lambda}_i^\mathrm{\scriptscriptstyle{T}} \boldsymbol{f}_t + \sigma_i(\boldsymbol{x}) \sigma_t(\boldsymbol{x}) U_{it}(\boldsymbol{x}), \ \ U_{it}(\boldsymbol{x}) \mid \boldsymbol{X}^{NT}, \boldsymbol{A}^N, \boldsymbol{B}^T \sim i.i.d. \ F_{U},
$$
where $U_{it}(\boldsymbol{x})$ is a zero mean random variable with marginal distribution $F_{U}$, which does not depend on $\boldsymbol{x}$.
This is special case of Assumption \ref{ass:model} with $Y_{it} = Y_{it}(\boldsymbol{X}_{it})$, $\boldsymbol{A}_i = \left(\boldsymbol{\lambda}_i, \{\sigma_i(\boldsymbol{x}) : \boldsymbol{x} \in \mathbb{X}\}\right),$ $\boldsymbol{B}_t = \left(\boldsymbol{f}_t, \{\sigma_t(\boldsymbol{x}) : \boldsymbol{x} \in \mathbb{X}\}\right)$, and $\boldsymbol{U}_{it} = U_{it}(\boldsymbol{X}_{it})$. The average effect of changing the covariate from $\boldsymbol{x}_0$ to $\boldsymbol{x}_1$ at $t$ is
$$
\mu_t(\boldsymbol{x}_1) - \mu_t(\boldsymbol{x}_0) = \mu_t(\boldsymbol{x}_1 \mid \{\boldsymbol{x}_1\}) - \mu_t(\boldsymbol{x}_0 \mid \{\boldsymbol{x}_1\}) = (\boldsymbol{x}_1 - \boldsymbol{x}_0)^\mathrm{\scriptscriptstyle{T}} \boldsymbol{\beta}.
$$
A version of this model was considered by \citet{KimOka2014} to analyze the effect of unilateral divorce laws on divorce rates in the U.S. This model encompasses the standard difference-in-differences model, $Y_{it}(\boldsymbol{x}) = \boldsymbol{x}^\mathrm{\scriptscriptstyle{T}} \boldsymbol{\beta} + \lambda_i + f_t + \sigma_i(\boldsymbol{x}) \sigma_t(\boldsymbol{x}) U_{it}(\boldsymbol{x})$, by setting $\boldsymbol{\lambda}_i = (\lambda_i,1)^{\mathrm{\scriptscriptstyle{T}}}$ and $\boldsymbol{f}_t = (1,f_t)^{\mathrm{\scriptscriptstyle{T}}}$.}
\end{example}
\begin{example}[Binary response model] \textnormal{Assume that the potential outcome $Y_{it}(\boldsymbol{x})$ is binary and generated by
\begin{align*}
Y_{it}(\boldsymbol{x}) &= \mathbbm{1}\{ m(\boldsymbol{x},\boldsymbol{A}_i,\boldsymbol{B}_t) \geq U_{it}(\boldsymbol{x}) \},
\quad
U_{it}(\boldsymbol{x}) \mid \boldsymbol{X}^{NT},\boldsymbol{A}^N,\boldsymbol{B}^T \sim i.i.d. \, {\cal U}(0,1),
\end{align*}
for some unknown function $m$. Here, assuming that $U_{it}(\boldsymbol{x})$ is uniform is a normalization, since $m$ can be arbitrary. This latent index model with unobserved effects is a special case of Assumption \ref{ass:model} with $Y_{it} = Y_{it}(\boldsymbol{X}_{it})$ and $\boldsymbol{U}_{it} = U_{it}(\boldsymbol{X}_{it})$. The ASFs at $\boldsymbol{x}$ and $t$ is
$$
\mu_t(\boldsymbol{x}) = \frac 1 {N} \sum_{i=1}^N m(\boldsymbol{x},\boldsymbol{A}_i,\boldsymbol{B}_t).
$$
Similar latent index models for count or censored responses are also covered by Assumption \ref{ass:model}. }
\end{example}
\begin{example}[Treatment effect factor model]\label{ex:factor} \textnormal{Assume that $\boldsymbol{X}_{it}$ contains only a binary treatment indicator, i.e., $\mathbb{X} = \{0,1\}$. The potential outcomes are generated by the linear factor model
$$
Y_{it}(\boldsymbol{x}) = \boldsymbol{\lambda}_i(\boldsymbol{x})^\mathrm{\scriptscriptstyle{T}} \boldsymbol{f}_t(\boldsymbol{x}) + \sigma_i(\boldsymbol{x}) \sigma_t(\boldsymbol{x}) U_{it}(\boldsymbol{x}), \ \ U_{it}(\boldsymbol{x}) \mid \boldsymbol{X}^{NT}, \boldsymbol{A}^N, \boldsymbol{B}^T \sim i.i.d. \ F_{U}, \ \boldsymbol{x} \in \mathbb{X},
$$
where $U_{it}(\boldsymbol{x})$ is a zero mean random variable with marginal distribution $F_{U}$, which does not depend on $\boldsymbol{x}$.
This is special case of Assumption \ref{ass:model} with $
Y_{it} = Y_{it}(\boldsymbol{X}_{it})
$, $\boldsymbol{A}_i = \left(\{\boldsymbol{\lambda}_i(\boldsymbol{x}), \sigma_i(\boldsymbol{x}) : \boldsymbol{x} \in \mathbb{X}\}\right)$, $\boldsymbol{B}_t = \left(\{\boldsymbol{f}_t(\boldsymbol{x}), \sigma_t(\boldsymbol{x}) : \boldsymbol{x} \in \mathbb{X}\}\right)$, and $\boldsymbol{U}_{it} = U_{it}(\boldsymbol{X}_{it})$. The average treatment effect at $t$ is
$$
\mu_t(1) - \mu_t(0)= \frac 1 {N} \sum_{i=1}^N [ \boldsymbol{\lambda}_i(1)^\mathrm{\scriptscriptstyle{T}} \boldsymbol{f}_t(1) - \boldsymbol{\lambda}_i(0)^\mathrm{\scriptscriptstyle{T}} \boldsymbol{f}_t(0)],
$$
and the average effect on the treated at $t$ is
$$
\mu_t(1 \mid \{1\}) - \mu_t(0 \mid \{1\})= \frac 1 {N_t(1)} \sum_{i=1}^N \mathbbm{1}\{\boldsymbol{X}_{it} = 1 \} [ \boldsymbol{\lambda}_i(1)^\mathrm{\scriptscriptstyle{T}} \boldsymbol{f}_t(1) - \boldsymbol{\lambda}_i(0)^\mathrm{\scriptscriptstyle{T}} \boldsymbol{f}_t(0)],
$$
provided that $N_t(1) = \sum_{i=1}^N \mathbbm{1}\{\boldsymbol{X}_{it} = 1 \} > 0$. Versions of this model have been considered by \citet{hsiao12}, \citet{gobillon16}, \citet{athey17}, \citet{li17}, \citet{xu17}, \citet{li18}, \citet{bai19matrix}, \citet{xiong19}, and \citet{chan20}.
Example \ref{ex:lfm} is a special case with $\boldsymbol{\lambda}_i(\boldsymbol{x})^\mathrm{\scriptscriptstyle{T}} \boldsymbol{f}_t(\boldsymbol{x}) = \boldsymbol{x}^{\mathrm{\scriptscriptstyle{T}}}\boldsymbol{\beta} + \boldsymbol{\lambda}_i^\mathrm{\scriptscriptstyle{T}} \boldsymbol{f}_t$. }
\end{example}
Throughout this paper we use standard panel data notation, with the two panel dimensions being denoted by units $i$ and time $t$. However,
one could also consider pseudo-panel or network applications of our results, where the two panel dimensions are denoted by $i$ and $j$, and
$Y_{ij}$ could, for example, be wage of worker $i$ in firm $j$, consumption of member $i$ in household $j$, a friendship indicator between individuals $i$ and $j$, or the volume of trade from country $i$ to country $j$.
The existing literature on two-way heterogeneity in network models usually either makes stronger parametric assumptions than we impose here
(e.g.~\citet{graham2017econometric}, \citet{dzemski2019empirical}, \citet{chen2020nonlinear}, \citet{zeleneev2020identification})
or uses stochastic blockmodels or graphon models, which typically ignore the effect of covariates
(e.g.~\citet{holland1983stochastic}, \citet{wolfe2013nonparametric}, \citet{gao2015rate}, \citet{auerbach2019identification}). Our methods of estimating non-parametric models
with two-way heterogeneity may therefore also be of interest in a network context.
\section{Estimation via Factor Structure Approximation}\label{sec:est}
\setcounter{equation}{0}
A natural starting point to estimate the effects in \eqref{eq:asf} and \eqref{eq:casf} is to use empirical analogs. This amounts to replacing $ \mathrm{E} \left[ Y_{it}(\boldsymbol{x}) \mid \boldsymbol{A}^N, \boldsymbol{B}^T \right]$ by an estimator. There are two complications with this approach. First, the potential outcome $Y_{it}(\boldsymbol{x})$ is not observable. We deal with this complication by noting that
\begin{multline*}
\mathrm{E} \left[ Y_{it}(\boldsymbol{x}) \mid \boldsymbol{A}^N, \boldsymbol{B}^T \right] = \mathrm{E} \left[ g(\boldsymbol{x}, \boldsymbol{A}_i, \boldsymbol{B}_t , \boldsymbol{U}_{it}(\boldsymbol{x}) ) \mid \boldsymbol{A}^N, \boldsymbol{B}^T \right] \\ = \mathrm{E} \left[ g(\boldsymbol{x}, \boldsymbol{A}_i, \boldsymbol{B}_t , \boldsymbol{U}_{it} ) \mid \boldsymbol{A}^N, \boldsymbol{B}^T \right] = \mathrm{E} \left[ g(\boldsymbol{x}, \boldsymbol{A}_i, \boldsymbol{B}_t , \boldsymbol{U}_{it} ) \mid \boldsymbol{X}_{it} = \boldsymbol{x}, \boldsymbol{A}_i, \boldsymbol{B}_t \right] \\ = \mathrm{E} \left[ Y_{it} \mid \boldsymbol{X}_{it} = \boldsymbol{x}, \boldsymbol{A}_i, \boldsymbol{B}_t \right],
\end{multline*}
under the rank similarity in \eqref{eq:rs} and Assumption \ref{ass:model}. Hence, we can write the expectation of the potential outcome as an expectation of the observed outcome.
The second complication is that $\boldsymbol{A}_i$ and $\boldsymbol{B}_t$ are not observable, so that we cannot directly estimate $\mathrm{E} \left[ Y_{it} \mid \boldsymbol{X}_{it} = \boldsymbol{x}, \boldsymbol{A}_i, \boldsymbol{B}_t \right]$. To deal with this complication, we start by noticing that
\begin{align}\label{eq:cef}
\begin{split}
\mathrm{E} \left[ Y_{it} \mid \boldsymbol{X}_{it} = \boldsymbol{x}, \boldsymbol{A}_i = \boldsymbol{a}, \boldsymbol{B}_t = \boldsymbol{b} \right] &= \mathrm{E} \left[ g(\boldsymbol{x}, \boldsymbol{a}, \boldsymbol{b} , \boldsymbol{U}_{it} ) \mid \boldsymbol{A}_i = \boldsymbol{a}, \boldsymbol{B}_t = \boldsymbol{b} \right] \\
&=: m(\boldsymbol{x}, \boldsymbol{a}, \boldsymbol{b}),
\end{split}
\end{align}
where the function $m$ does not vary with $i$ and $t$, by implication \eqref{ass:Main3} of Assumption \ref{ass:model}. We show next how this function can be approximated and estimated using a low-rank factor structure.
\subsection{Low-rank factor structure approximation}
For ease of exposition, we assume in the rest of the paper that the covariate vector $\boldsymbol{X}_{it}$ includes only a binary treatment and $\mathbb{X} = \{0,1\}$. Accordingly, we denote the covariate and its values by $X_{it}$ and $x$ instead of $\boldsymbol{X}_{it}$ and $\boldsymbol{x}$. In what follows, $x$ denotes a generic element of $\mathbb{X}$ and all the assumptions and results hold for all $x \in \mathbb{X}_1 \subseteq \mathbb{X}$, where $\mathbb{X}_1 = \mathbb{X}$ if we are interested in the entire population, $\mathbb{X}_1 = \{0\}$ if we are interested in the treated subpopulation, and $\mathbb{X}_1 = \{1\}$ if we are interested in the untreated subpopulation.
The approximation that we propose is based on the singular value decomposition of the function $(\boldsymbol{a}, \boldsymbol{b}) \mapsto m(x,\boldsymbol{a},\boldsymbol{b}) $ for each $x \in \mathbb{X}$. We make two assumptions on this decomposition. The first assumption is a sampling condition on the unobserved effects that will be useful to define a norm for the eigenfunctions.
\begin{assumption}[Sampling of $\boldsymbol{A}_i$ and $\boldsymbol{B}_t$]~
\label{ass:SamplingAB}
(i) $\boldsymbol{A}_i$ is independent and identically distributed across $i \in \mathbb{N}$ with distribution $F_{\boldsymbol{A}}$, (ii) $\boldsymbol{B}_t$ is independent and identically distributed over $t \in \mathbb{T}$ with distribution $F_{\boldsymbol{B}}$, and (iii) $\boldsymbol{A}_i$ and $\boldsymbol{B}_t$ are independent for all $i,t$.
\end{assumption}
For simplicity we consider the case where both $\boldsymbol{A}_i$ and $\boldsymbol{B}_t$ are independently distributed across $i$ and over $t$,
but since we consider asymptotic sequences where both $N$ and $T$ become large one could also allow for
appropriate weak dependence across both $i$ and $t$. Formalizing this weak dependence would complicate both the assumption
and the proof of the following results, which is why we decided to stick to independence in our presentation here.
The next assumption is a regularity condition on the function $ m(x,\boldsymbol{a},\boldsymbol{b}) $.
\begin{assumption}[Smoothness of $(\boldsymbol{a},\boldsymbol{b}) \mapsto m(x,\boldsymbol{a}, \boldsymbol{b})$]
\label{ass:Smoothness}
The function $(\boldsymbol{a},\boldsymbol{b}) \mapsto m(x,\boldsymbol{a}, \boldsymbol{b})$ admits a singular value decomposition
\begin{align*}
m(x,\boldsymbol{a}, \boldsymbol{b}) = \sum_{j=1}^\infty \, s_j(x) \, u_j(x,\boldsymbol{a}) \, v_j(x,\boldsymbol{b}) ,
\end{align*}
under the $L_2(F_{\boldsymbol{A}} \times F_{\boldsymbol{B}})$ norm, where the eigenfunctions $u_j(x,\boldsymbol{a})$ and $v_j(x,\boldsymbol{b})$ are orthonormal, i.e.,
\begin{align*}
\mathbb{E} \, u_j(x,\boldsymbol{A}_i)^2 &= 1 ,
&
\mathbb{E} \, u_j(x,\boldsymbol{A}_i) u_k(x,\boldsymbol{A}_i) &=0,\\
\mathbb{E} \, v_j(x,\boldsymbol{B}_t)^2 &= 1,
&
\mathbb{E} \, v_j(x,\boldsymbol{B}_t) v_k(x,\boldsymbol{B}_t) &=0,
& j \neq k&\in \{1,2,3\ldots\},
\end{align*}
and the singular values $s_1(x) \geq s_2(x) \geq s_3(x) \geq \ldots \geq 0$ satisfy
\begin{align*}
\sum_{j=1}^\infty s_j(x) &< \infty.
\end{align*}
\end{assumption}
There is a large literature on singular value decompositions of functions, which shows that, under appropriate conditions, the singular values satisfy $s_j(x) \lesssim j^{-\alpha}$,\footnote{
Here, $s_j(x) \lesssim j^{-\alpha}$ means that there exists a constant $c>0$ such that $s_j(x) \leq c\, j^{-\alpha}$, for all $j$.
} where the decay coefficient $\alpha$
depends on the dimensions of the arguments $\boldsymbol{a}$, $\boldsymbol{b}$, and on the smoothness of $(\boldsymbol{a},\boldsymbol{b}) \mapsto m(x,\boldsymbol{a}, \boldsymbol{b}) $. For sufficiently smooth functions,
$\alpha>1$ and therefore $ \sum_{j=1}^\infty s_j(x) < \infty$. For example, if $(\boldsymbol{a}, \boldsymbol{b}) \mapsto m(x,\boldsymbol{a},\boldsymbol{b})$ is continuously differentiable up to order $s$ and $\mathbb{A}$ and $\mathbb{B}$ are compact, then
$$
s_j(x) \lesssim j^{- \frac{s}{d_a \wedge d_b} },
$$
by Theorem 3.3 of \citet{gh13}, where $d_a \wedge d_b$ is the minimum of $d_a$ and $d_b$. This implies that $ \sum_{j=1}^\infty s_j(x) < \infty$ if $s > d_a \wedge d_b $.
Assumption~\ref{ass:Smoothness} is therefore a high-level
smoothness assumption on $(\boldsymbol{a},\boldsymbol{b}) \mapsto m(x,\boldsymbol{a}, \boldsymbol{b})$,
very similar to the Assumption~2.2. in \citet{menzel2018bootstrap}, where an analogous condition on the
singular values is imposed, with the same aim of controlling the behaviour of a function of unobserved two-dimensional heterogeneity.
The formulation of this smoothness assumption is convenient for our purposes, because it immediately leads to a low-rank approximation of
$m(x,\boldsymbol{a},\boldsymbol{b})$.
The low-rank approximation truncates the singular value decomposition to the first $R$ elements,
\begin{equation}\label{eq:approx}
m(x,\boldsymbol{a},\boldsymbol{b}) = \sum_{j=1}^{\infty} \, \underbrace{s_j(x)^{1/2} u_j(x,\boldsymbol{a})}_{=:\phi_{j}(x,\boldsymbol{a}) }
\, \underbrace{s_j(x)^{1/2} v_j(x,\boldsymbol{b}) }_{=: \psi_{j}(x,\boldsymbol{b})} = \sum_{j=1}^{R} \phi_{j}(x,\boldsymbol{a}) \psi_{j}(x,\boldsymbol{b}) + \zeta_R(x,\boldsymbol{a},\boldsymbol{b}).
\end{equation}
The first term is the approximation and the second term is the approximation error. Under Assumption~\ref{ass:Smoothness},
$$
\mathrm{E} \ \zeta_R(x,\boldsymbol{A}_i,\boldsymbol{B}_t)^2 \to 0 \ \ \text{ as } \ \ R \to \infty.
$$
In other words, the approximation error can be made negligible by increasing the truncation point $R$. For example, if $s_j(x) \lesssim j^{-\alpha}$ with $\alpha > 1$, then
\begin{multline*}
\mathrm{E} \ \zeta_R(x,\boldsymbol{A}_i,\boldsymbol{B}_t)^2 = \mathrm{E}\left[ \sum_{j=R+1}^\infty \, s_j(x) \, u_j(x,\boldsymbol{A}_i) \, v_j(x,\boldsymbol{B}_t) \right]^2
\\ = \sum_{j,k =R+1}^\infty s_j(x)s_k(x) \mathrm{E} \left[ u_j(x,\boldsymbol{A}_i)u_k(x,\boldsymbol{A}_i) \right] \, \mathrm{E} \left[ v_j(x,\boldsymbol{B}_t)v_k(x,\boldsymbol{B}_t) \right] \\ = \sum_{j =R+1}^\infty s_j(x)^2 \lesssim \sum_{j=R+1}^\infty j^{-2\alpha} \leq \int_{R}^{\infty} j^{-2\alpha} \mathrm{d} j \lesssim R^{1 - 2 \alpha },
\end{multline*}
by Assumptions \ref{ass:SamplingAB} and \ref{ass:Smoothness}.
Hence, $\zeta_R(x,\boldsymbol{A}_i,\boldsymbol{B}_t)$ converges in mean square to zero at a polynomial rate with $R$.
Combining \eqref{eq:cef} and \eqref{eq:approx}, we obtain the approximate factor model
\begin{equation}\label{eq:afm}
Y_{it} = \boldsymbol{\lambda}_{i}(X_{it})^\mathrm{\scriptscriptstyle{T}} \boldsymbol{f}_{t}(X_{it}) + \zeta_R(X_{it},\boldsymbol{A}_i,\boldsymbol{B}_t) + E_{it}, \ \ E_{it} := Y_{it} - \mathrm{E} \left[ Y_{it} \mid X_{it}, \boldsymbol{A}_i, \boldsymbol{B}_t \right],
\end{equation}
where $\boldsymbol{\lambda}_{i}(x) = [\phi_{1}(x,\boldsymbol{A}_i), \ldots, \phi_{R}(x,\boldsymbol{A}_i)]^\mathrm{\scriptscriptstyle{T}}$, $\boldsymbol{f}_{t}(x) = [\psi_{1}(x,\boldsymbol{B}_t), \ldots, \psi_{R}(x,\boldsymbol{B}_t)]^\mathrm{\scriptscriptstyle{T}}$, and the composite error $\nu_{it} := \zeta_R(X_{it},\boldsymbol{A}_i,\boldsymbol{B}_t) + E_{it}$ contains the approximation error, $\zeta_R(X_{it},\boldsymbol{A}_i,\boldsymbol{B}_t) $, and the conditional expectation error, $E_{it}$.
The factor structure can be seen as a series or sieve approximation to the function $(\boldsymbol{a},\boldsymbol{b}) \mapsto m(x,\boldsymbol{a},\boldsymbol{b})$ with basis functions $\{\phi_j(x,\boldsymbol{a}) \psi_j(x,\boldsymbol{b})\}_{j=1}^{\infty}$ if we let $R=R_{N,T}$ to grow with $N$ and $T$ such that $\zeta_{R}(x,\boldsymbol{a},\boldsymbol{b})$ vanishes as $N,T \to \infty$. The factor structure approximation is exact in some cases for fixed $R$. For instance, in Example \ref{ex:factor}
\[
m(x,\boldsymbol{A}_i,\boldsymbol{B}_t) = \boldsymbol{\lambda}_i(x)^\mathrm{\scriptscriptstyle{T}} \boldsymbol{f}_t(x),
\]
so that $\zeta_R(x,\boldsymbol{A}_i,\boldsymbol{B}_t) = 0$, a.s., if $R$ is greater or equal to the number of factors.
In the model \eqref{eq:afm} the factor structure changes with the treatment level. In other words, we have a different pure factor model for each $x \in \mathbb{X}$, that is
$$
Y_{it} = \boldsymbol{\lambda}_{i}(x)^\mathrm{\scriptscriptstyle{T}} \boldsymbol{f}_{t}(x) + \nu_{it} \text{ if } X_{it} = x.
$$
This observation leads to our first estimation strategy where the data is partitioned by the treatment level and separate factors and factor loadings are estimated in each element of the partition by solving the least squares program
\begin{equation}\label{eq:pca}
\min_{\{\boldsymbol{\lambda}_{i}\}_{i=1}^N, \{\boldsymbol{f}_{t}\}_{t=1}^T} \frac{1}{2} \sum_{i=1}^N \sum_{t=1}^T D_{it}(x) \left(Y_{it} - \boldsymbol{\lambda}_{i}^\mathrm{\scriptscriptstyle{T}} \boldsymbol{f}_{t}\right)^2,
\end{equation}
where $D_{it}(x) := \mathbbm{1}\{X_{it} = x\} $.
Unfortunately, we cannot solve this problem using standard principal component analysis due to the presence of missing data, that is, each observational unit $(i,t)$ is not available at all treatment levels. In the next section, we apply matrix completion methods to deal with this problem.
\subsection{Estimation by matrix completion methods}
We start by expressing the program \eqref{eq:pca} in matrix form. Let $\boldsymbol{\Gamma}^R(x) = \boldsymbol{\lambda}^N(x) \boldsymbol{f}^T(x)^\mathrm{\scriptscriptstyle{T}}$, where $\boldsymbol{\lambda}^N(x) = [\boldsymbol{\lambda}_1(x), \ldots, \boldsymbol{\lambda}_N(x)]^\mathrm{\scriptscriptstyle{T}}$, a $N \times R$ matrix of factor loadings, and $\boldsymbol{f}^T(x) = [\boldsymbol{f}_1(x), \ldots, \boldsymbol{f}_T(x)]^\mathrm{\scriptscriptstyle{T}}$, a $T \times R$ matrix of factors. The least squares estimator of $\boldsymbol{\Gamma}^R(x)$ is the $N \times T$ matrix $\boldsymbol{\Gamma} $ with typical element $\Gamma_{it}$ that solves
\begin{equation}\label{eq:mc}
\min_{\{\boldsymbol{\Gamma} \in \mathbb{R}^{N \times T}: \operatorname*{rank}(\boldsymbol{\Gamma}) \leq R \}} \frac{1}{2} \sum_{i=1}^N \sum_{t=1}^T D_{it}(x) \left(Y_{it} - \Gamma_{it} \right)^2.
\end{equation}
Let $\boldsymbol{Y}(x)$ be a $N \times T$ matrix whose $(i,t)$ element is $Y_{it}$ if $X_{it} = x$ and is missing otherwise. The previous program is closely related to the problem of completing the missing entries of $\boldsymbol{Y}(x)$ using a low rank approximation matrix $\boldsymbol{\Gamma}^R(x)$ \citet{rennie05,candes09,candes10}. This connection was previously noticed by \citet{athey17} and \citet{amjad18} in the context of treatment effects models. The solution is the $N \times T$ matrix of rank $R$ whose entries are the closest in the mean squared error sense to the corresponding entries of $\boldsymbol{Y}(x)$.
The previous program is combinatorially hard because of the constraint in the rank of the matrix \citet{srebro03}.
Following \citet{fazel03} we consider the convex relaxation of this program.
Let $\|{\bf M}\|_\infty$ be the spectral norm of a $\mathbb{R}^{N \times T}$-matrix ${\bf M}$, and define the nuclear norm (also called trace norm)
of $\boldsymbol{\Gamma}$ as the corresponding dual norm
$\|\boldsymbol{\Gamma}\|_1 := \max_{\left\{ {\bf M} \in \mathbb{R}^{N \times T} \, : \, \|{\bf M}\|_\infty \leq 1 \right\}} {\rm Tr}\left( {\bf M}' \boldsymbol{\Gamma} \right)$.
This nuclear norm can equivalently be defined as the sum of the singular values of $\boldsymbol{\Gamma}$.
Using this norm we can write the convex relaxation of the program \eqref{eq:mc} as follows,
$$
\min_{\{\boldsymbol{\Gamma} \in \mathbb{R}^{N \times T}: \|\boldsymbol{\Gamma}\|_1 \leq R_1 \}} \frac{1}{2}\sum_{i=1}^N \sum_{t=1}^T D_{it}(x) \left(Y_{it} - \Gamma_{it} \right)^2,
$$
where $R_1$ is a positive constant such that $R = f(R_1)$, where $f$ is an increasing function. Hence, $\zeta_R(x,\boldsymbol{A}_i,\boldsymbol{B}_t)$ vanishes in mean square as $R_1 \to \infty$. We replace the rank constraint, $\operatorname*{rank}(\boldsymbol{\Gamma}) \leq R$, by a constraint on the nuclear norm of the matrix, $\|\boldsymbol{\Gamma}\|_1 \leq R_1$, i.e. we replace a constraint in the number of nonzero singular values by a constraint in the sum of singular values. This program is convex in $\boldsymbol{\Gamma}$ and can be reformulated in Lagrange form as
\begin{equation}\label{eq:nn}
\min_{\{\boldsymbol{\Gamma} \in \mathbb{R}^{N \times T}\}} \frac{1}{2}\sum_{i=1}^N \sum_{t=1}^T D_{it}(x) \left(Y_{it} - \Gamma_{it} \right)^2 + \rho(R_1) \|\boldsymbol{\Gamma}\|_1,
\end{equation}
where $\rho(R_1) \geq 0$ is a regularization parameter, which is a one-to-one increasing function of $R_1$. There exist efficient algorithms to solve this program \citet{mazumder10}.
Let $\widehat \boldsymbol{\Gamma}(x)$ be a solution to \eqref{eq:nn} with typical element $\widehat \Gamma_{it}(x)$. Then, we can form estimators of the ASF and CASF as
$$
\widehat \mu_{t}(x) = \frac 1 {N} \sum_{i=1}^N \left[ D_{it}(x) Y_{it} + \{1-D_{it}(x)\} \widehat \Gamma_{it}(x) \right],
$$
and
$$
\widehat \mu_{t}(x \mid \{x_0\}) = \frac {\sum_{i=1}^N D_{it}(x_0) \left[ D_{it}(x)Y_{it} + \{1- D_{it}(x) \} \widehat \Gamma_{it}(x)\right]} {\sum_{i=1}^N D_{it}(x_0)}.
$$
In the next section, we provide conditions under which these estimators are consistent using asymptotic sequences where $N, T \to \infty$. These estimators, however, might display shrinkage biases in finite samples due to the nuclear norm regularization \citet{cai10,ma11,bai19joe}. We propose two matching procedures to debias the estimator in Section \ref{sec:debias}.
\subsection{Consistency of Matrix Completion Estimator}\label{sec:mc}
Let $\boldsymbol{\Gamma}^{\infty}(x)$ be the $N \times T$ matrix with typical element $ \Gamma^{\infty}_{it}(x) = m(x, \boldsymbol{A}_i, \boldsymbol{B}_t)$ and $\boldsymbol{E}(x)$ be the $N \times T$ matrix with typical element
\begin{align}
E_{it}(x) &:= \left\{
\begin{array}{ll}
E_{it} = Y_{it} - \Gamma^{\infty}_{it}(x) & \text{if $X_{it}=x$,}
\\
0 & \text{otherwise.}
\end{array}
\right.
\label{DefE}
\end{align}
Note that $\boldsymbol{\Gamma}^{\infty}(x) = \lim_{R\to \infty} \boldsymbol{\Gamma}^{R}(x)$ a.s.
Furthermore, we introduce the notation $\mathbb{D}(x) = \{ (i,t) \in \mathbb{N} \times \mathbb{T} \, : \, X_{it} = x \}$,
and $n(x)= |\mathbb{D}(x)|$ for the number of observations with $X_{it} = x$.
Recall that
\begin{align}
\widehat \boldsymbol{\Gamma}(x)
&\in
\operatorname*{argmin}_{ \boldsymbol{\Gamma} \in \mathbb{R}^{N \times T}}
Q_{NT}(\boldsymbol{\Gamma},\rho,x) ,
&
Q_{NT}( \boldsymbol{\Gamma},\rho,x) &=
\frac 1 2 \sum_{(i,t) \in \mathbb{D}(x) }
\left(Y_{it} - \Gamma_{it} \right)^2 + \rho \| \boldsymbol{\Gamma} \|_1,
\label{DefEstimator}
\end{align}
where $\rho := \rho(R_1)$.
Here, if the $ \operatorname*{argmin}$ over $ \boldsymbol{\Gamma} \in \mathbb{R}^{N \times T} $ is not unique, then we can choose $ \widehat \boldsymbol{\Gamma}(x)$
arbitrarily from the set of minimizers --- our results are not affected by that, we only require that
$ Q_{NT}( \widehat \boldsymbol{\Gamma}(x),\rho,x) \leq Q_{NT}(\boldsymbol{\Gamma},\rho,x) $,
for all $ \boldsymbol{\Gamma} \in \mathbb{R}^{N \times T}$. We want to show that $ \widehat \boldsymbol{\Gamma}(x)$ converges to $\boldsymbol{\Gamma}^{\infty}(x)$ as $N,T \to \infty$ in some sense such that $ \widehat \mu(x) - \mu(x) = o_P(1)$. For that we require additional assumptions.
\begin{assumption}[Error Moments]~
\label{ass:SamplingE}
Conditional on $\boldsymbol{X}^{NT}$, $\boldsymbol{A}^N$ and $\boldsymbol{B}^T$, $E_{it}(x)$ is independent across $ (i,t) \in \mathbb{D}(x) $,
and there exists a constant $b < \infty$ that does not depend on $i$, $t$, $N$, $T$, such that
$$
\mathrm{E}\left[E_{it}(x)^4 \mid \boldsymbol{A}^N, \boldsymbol{B}^T, \boldsymbol{X}^{NT} \right] \leq b.
$$
Furthermore, we assume that $n(x)^{-1} \sum_{(i,t) \in \mathbb{D}(x) }
\Gamma^{\infty}_{it}(x)^2 = O_P(1)$.
\end{assumption}
For the purpose of showing Lemma~\ref{lemma:consistency1} and Theorem~\ref{th:MAIN1}
we could alternatively replace Assumption~\ref{ass:SamplingE} by the two high-level conditions:
\begin{align*}
\frac 2 {n(x)} \sum_{(i,t) \in \mathbb{D}(x) } \Gamma^{\infty}_{it}(x) E_{it} &= o_P(1) ,
&
\| \boldsymbol{E}(x) \|_\infty &= O_P\left( \sqrt{N+T} \right),
\end{align*}
where again $\|\cdot\|_\infty$ denotes the spectral norm.
The first of those conditions is implied by Assumption~\ref{ass:SamplingE}
through application of the weak law of large numbers,
while the second follows, for example, by the spectral norm inequality in \citet{latala05}.
In principle, we could still derive those high-level conditions if we allowed for appropriate weak dependence of $E_{it}(x)$ across $i$
and over $t$, but we again focus on the independent case for simplicity of presentation.
We first provide a consistency result for the entries of $\widehat \boldsymbol{\Gamma}(x)$ that correspond to the observed values of $\boldsymbol{Y}(x)$.
\begin{lemma}
\label{lemma:consistency1}
Let the Assumptions~\ref{ass:SamplingAB}, \ref{ass:Smoothness} and \ref{ass:SamplingE} hold,
and assume that $\rho = \rho_{NT}$ is chosen such that $\rho_{NT} / \sqrt{N+T} \rightarrow \infty$ and $\rho_{NT} \sqrt{NT} / n(x) \rightarrow 0$ as $N,T \rightarrow \infty$.
Then,
\begin{align*}
\frac 1 { n(x)} \sum_{(i,t) \in \mathbb{D}(x) }
\left[ \widehat \Gamma_{it}(x) - \Gamma^{\infty}_{it}(x) \right]^2
= o_P(1).
\end{align*}
\end{lemma}
A necessary condition for the existence of the sequence $\rho = \rho_{NT}$ in Lemma~\ref{lemma:consistency1} is
$n(x) / \sqrt{(N+T)NT} \rightarrow \infty$, that is, the fraction $n(x) / (NT)$ of observations with $X_{it} = x$ can converge to zero,
but not too fast. Apart from that, Lemma~\ref{lemma:consistency1} does not restrict the assignment process that determines $\boldsymbol{X}^{NT}$.
Notice also that Lemma~\ref{lemma:consistency1} does not require Assumption \ref{ass:model} because $\Gamma^{\infty}(x)$ is a reduced-form parameter.
Applying the Cauchy-Schwarz inequality
$$\left( \frac 1 { n(x)} \sum_{(i,t) \in \mathbb{D}(x) } a_{it} \right)^2 \leq \frac 1 { n(x)} \sum_{(i,t) \in \mathbb{D}(x) } a_{it}^2$$
with $a_{it} = \widehat \Gamma_{it}(x) - \Gamma^{\infty}_{it}(x)$,
Lemma~\ref{lemma:consistency1}
guarantees that
$$\frac 1 { n(x)} \sum_{(i,t) \in \mathbb{D}(x) }
\left[ \widehat \Gamma_{it}(x) - \Gamma^{\infty}_{it}(x) \right]
= o_P(1).$$
Nevertheless, Lemma~\ref{lemma:consistency1}
is not directly useful to show the consistency of the estimators of the ASF, because it only guarantees $L_2$-consistency
of $ \widehat \boldsymbol{\Gamma}(x)$ over the set of entries $(i,t)$ for which $X_{it} = x$. Those are exactly the observations
for which an unbiased estimator of $ \Gamma^{\infty}_{it}(x) = m(x,\boldsymbol{A}_i, \boldsymbol{B}_t)$ is already available, namely $Y_{it}$.
The consistency result we would like to obtain is
\begin{align}
\frac 1 { NT} \sum_{i=1}^N \sum_{t=1}^T
\left[ \widehat \Gamma_{it}(x) - \Gamma^{\infty}_{it}(x) \right]^2 = o_P(1) ,
\label{ConsistencyStrong}
\end{align}
but such a result will certainly require stronger assumptions on $\boldsymbol{X}^{NT}$ than we have imposed so far.
The existing literature on matrix completion relies on the concept
of \textit{restricted strong convexity} to derive \eqref{ConsistencyStrong}. This approach shows that
under certain conditions on a $\mathbb{R}^{N \times T}$-matrix $\boldsymbol{M}$ with entries $M_{it}$, and
on $\boldsymbol{X}^{NT}$ (which determines the set $\mathbb{D}(x)$), there exists a constant $c>0$ such that with high probability
\begin{align*}
\frac 1 { NT} \sum_{i=1}^N \sum_{t=1}^T M_{it}^2
\leq \frac c { n(x)} \sum_{(i,t) \in \mathbb{D}(x) } M_{it}^2 .
\end{align*}
See Theorem 1 in \citet{negahban2012restricted}, Lemma 12 in \citet{klopp2014noisy}, and Lemma 3 in \citet{athey17}.
Thus, if $M_{it} = \widehat \Gamma_{it}(x) - \Gamma^{\infty}_{it}(x)$ and $\boldsymbol{X}^{NT}$ satisfy restricted strong convexity, then \eqref{ConsistencyStrong} would follow from Lemma~\ref{lemma:consistency1}.
We pursue a different strategy than the existing matrix completion literature
to show that
$$
\widehat \mu(x) := \frac 1 {T} \sum_{t=1}^T \widehat \mu_{t}(x)
= \frac 1 {NT} \sum_{i=1}^N \sum_{t=1}^T D_{it}(x) \, Y_{it} + \frac 1 {NT} \sum_{i=1}^N \sum_{t=1}^T [1-D_{it}(x) ] \, \widehat \Gamma_{it}(x)
$$
is a consistent estimator of $(NT)^{-1} \sum_{i=1}^N \sum_{t=1}^T \Gamma^{\infty}_{it}$,
which under Assumption~\ref{ass:model} is equal to $\mu(x)$ defined in \eqref{eq:asf}.
We believe that our approach is simpler in the setting of this paper where
$ \Gamma^{\infty}_{it}(x)$ is not necessarily of low-rank.
In particular, we do not aim to show \eqref{ConsistencyStrong}, but instead we derive consistency of $\widehat \mu(x)$ directly.
However, the following theorem still requires additional assumptions on the assignment process that determines $\boldsymbol{X}^{NT}$,
in the same way that
additional conditions on $\boldsymbol{X}^{NT}$ are required to verify restricted strong convexity.
For simplicity, we focus on consistency of $ \widehat \mu(x)$ in the main text,
but results for more general weighted averages of the form $ (NT)^{-1} \sum_{i=1}^N \sum_{t=1}^T W_{it}(x) \, \Gamma^{\infty}_{it}(x)$,
with known weights $W_{it}(x) \in \mathbb{R}$,
are presented in the appendix. For example, in the case of the treatment effects on the treated that we consider in the empirical application of Section \ref{sec:Election}, $W_{it}(x) = n(1)^{-1} X_{it}$.
\begin{theorem}
\label{th:MAIN1}
Let the Assumptions~\ref{ass:model}, \ref{ass:SamplingAB}, \ref{ass:Smoothness} and \ref{ass:SamplingE} hold.
Consider $N,T \rightarrow \infty$ at the same rate,
and let $\rho = \rho_{NT}$ be chosen such that $\rho_{NT} / \sqrt{N+T} \rightarrow \infty$
and $\rho_{NT} / \sqrt{NT} \rightarrow 0$.
Let
$ P_{it}(x) = \Pr \left( X_{it} = x \mid \boldsymbol{A}^N, \boldsymbol{B}^T \right)$,
and assume
that $ (NT)^{-1} \sum_{i=1}^N \sum_{t=1}^T P_{it}^{-1}(x) = O_P(1)$.
Let $\boldsymbol{G}(x)$ be the $N \times T$ matrix with entries $ G_{it}(x) = P_{it}^{-1}(x) (D_{it}(x) - P_{it}(x)) $,
and assume that $ \| \boldsymbol{G}(x) \|_\infty = O_P(\sqrt{N+T})$, and
\begin{align}
\frac 1 {NT} \sum_{i=1}^N \sum_{t=1}^T P_{it}^{-1}(x) \, G_{it}(x) &= o_P(1) ,
&
\frac {1} {NT} \sum_{i=1}^N \sum_{t=1}^T \Gamma^{\infty}_{it}(x) \, G_{it}(x) &= o_P(1) .
\label{ConditionsGit}
\end{align}
Then,
\begin{align*}
\widehat \mu(x) = \mu(x) + o_P(1).
\end{align*}
\end{theorem}
\medskip
To interpret the conditions in Theorem~\ref{th:MAIN1}, notice that
due to the definitions $D_{it}(x) = \mathbbm{1}\{X_{it} = x\} $ and $ P_{it}(x) = \Pr \left( X_{it} = x \mid \boldsymbol{A}^N, \boldsymbol{B}^T \right)$,
$\mathrm{E} \left[ G_{it}(x) \mid \boldsymbol{A}^N, \boldsymbol{B}^T \right] = 0$ by construction, and $G_{it}(x) $ therefore plays a role very similar
to the error term $E_{it}(x)$. In particular, the conditions in \eqref{ConditionsGit} can be verified by a weak law of large numbers, as long as $ P_{it}^{-1}(x)$ is not too large, and $G_{it}(x)$ is not too strongly correlated across $i$ and over $t$.
Regarding the condition on the spectral norm $ \| \boldsymbol{G}(x) \|_\infty = O_P(\sqrt{N+T})$, there are many results in the random-matrix
theory literature that show this rate for mean-zero random matrices $\boldsymbol{G}(x) $, see, for example, \citet{Geman1980}, \citet{Silverstein1989}, \citet{BaiSilvYin1988},
\citet{BaiKrishYin1988}. In particular, if $G_{it}(x)$ is independent across both $i$ and $t$, then this rate result follows
from the very elegant spectral norm inequality in \citet{latala05}, see the proof of Lemma~\ref{lemma:consistency1} in the appendix,
where we apply that inequality to $E_{it}(x)$. However, that simple argument would require $X_{it} $ to be independently distributed across $i$ and $t$,
conditional on $\boldsymbol{A}^N$, $\boldsymbol{B}^T$.
More generally, we expect $ \| \boldsymbol{G}(x) \|_\infty = O_P(\sqrt{N+T})$ to hold whenever
the matrix entries $ G_{it}(x) $ have zero mean, sufficiently bounded moments,
and weak correlation across both $i$ and $t$, see Section S.2
of the supplementary material of \citet{MoonWeidner2017} for details.
We have thus shown that consistent estimates for ASFs can be obtained via the matrix completion estimator
even if the estimand $ \Gamma^{\infty}_{it}(x) = m(x, \boldsymbol{A}_i, \boldsymbol{B}_t)$ itself is not of low rank.
This is the main technical result of this paper. However, inference on $\mu(x)$ based on $\widehat \mu(x) $ can be problematic, because $\widehat \mu(x) $ is subject to both
low-rank approximation and shrinkage biases.
The low-rank approximation bias is due to the approximation error $ \zeta_R(x,\boldsymbol{a},\boldsymbol{b})$
in the decomposition of $m(x,\boldsymbol{a},\boldsymbol{b})$ in equation \eqref{eq:approx}.
The shrinkage bias comes from bias in
$\widehat \boldsymbol{\Gamma}(x)$ due to the presence of the nuclear norm penalization in the objective function of \eqref{DefEstimator}. To isolate this bias, consider a simple case where $Y_{it}(x)$ follows a deterministic pure factor model
$$
Y_{it}(x) = \Gamma_{it}(x) = \sum_{j=1}^R s_j(x) u_j(x,\boldsymbol{A}_i) v_j(x,\boldsymbol{B}_i).
$$
Then, the matrix completion estimator of $ \Gamma_{it}(x)$ in \eqref{DefEstimator} yields
$$
\widehat \Gamma_{it}(x) = \sum_{j=1}^R [s_j(x) - \rho]_{+} u_j(x,\boldsymbol{A}_i) v_j(x,\boldsymbol{B}_i)
$$
where $[z]_{+} = \max(z,0)$. Compared to $\boldsymbol{\Gamma}(x)$, $\widehat \boldsymbol{\Gamma}(x)$ has the same eigenvectors but the singular values are shrunk toward zero. This argument carries over to the case where $Y_{it}(x)$ follows an approximate factor structure \citet{cai10,ma11,bai19joe}. Because of these biases, we explore alternative estimates for $\mu(x)$ in Section~\ref{sec:debias}.
\subsection{Covariates and fixed effects}
As we mentioned in Section \ref{sec:model}, exogenous covariates can be incorporated by conditioning on their values. This method can produce very noisy estimators in small samples unless the covariates take only on few values. Here we consider a semiparametric version of the model that imposes additivity in the effect of the exogenous covariates, which may be continuous, discrete or mixed. It also allows for additive unobserved individual and time effects that might vary across the covariate level $x$. These effects can be subsumed in the factor structure, but are usually considered separately in empirical analysis as the estimators perform better without regularizing them \citet{athey17}.
Let $\boldsymbol{C}_{it}$ be a $d_c$-vector of covariates, $\boldsymbol{\alpha}(x) = (\alpha_1(x),\ldots, \alpha_N(x))$ be a $N$-vector of individual effects and $\boldsymbol{\delta}(x) = (\delta_1(x), \ldots, \delta_T(x))$ be a $T$-vector of time effects. Then, we can replace the program \eqref{eq:nn} by
\begin{equation*}
\min_{\{\boldsymbol{\beta} \in \mathbb{R}^{d_c}, \boldsymbol{\alpha} \in \mathbb{R}^N, \boldsymbol{\delta} \in \mathbb{R}^T, \boldsymbol{\Gamma} \in \mathbb{R}^{N \times T}\}} \sum_{i=1}^N \sum_{t=1}^T \mathbbm{1}\{X_{it} = x\} \left(Y_{it} - \boldsymbol{C}_{it}^\mathrm{\scriptscriptstyle{T}} \boldsymbol{\beta} - \alpha_i - \delta_t - \Gamma_{it} \right)^2 + \rho(R_1) \|\boldsymbol{\Gamma}\|_1,
\end{equation*}
\citet{chernozhukov18}, \citet{moon18} and \citet{beyhum19} provide algorithms to solve this program.
Let $\widehat \boldsymbol{\beta}(x)$, $\widehat \boldsymbol{\alpha}(x) = (\widehat \alpha_1(x),\ldots, \widehat \alpha_N(x))$, $\widehat \boldsymbol{\delta}(x) = (\widehat \delta_1(x), \ldots, \widehat \delta_T(x))$, and $\widehat \boldsymbol{\Gamma}(x)$ be the solution of the previous program. We can form estimators of the ASF and CASF as
$$
\widehat \mu_{t}(x) = \frac 1 {N} \sum_{i=1}^N \left[ \mathbbm{1}\{X_{it} = x\} Y_{it} + \mathbbm{1}\{X_{it} \neq x\} \left\{ \boldsymbol{C}_{it}^\mathrm{\scriptscriptstyle{T}} \widehat \boldsymbol{\beta}(x) + \widehat \alpha_i(x) + \widehat \delta_t(x) + \widehat \Gamma_{it}(x) \right\} \right],
$$
and
\begin{multline*}
\widehat \mu_{t}(x \mid \{x_0\}) = \\ \frac {\sum_{i=1}^N \left[ \mathbbm{1}\{X_{it} = x_0 = x\} Y_{it} + \mathbbm{1}\{X_{it} = x_0 \neq x\} \left\{ \boldsymbol{C}_{it}^\mathrm{\scriptscriptstyle{T}} \widehat \boldsymbol{\beta}(x) + \widehat \alpha_i(x) + \widehat \delta_t(x) + \widehat \Gamma_{it}(x) \right\} \right]} {\sum_{i=1}^N \mathbbm{1}\{X_{it} = x_0 \}}.
\end{multline*}
\section{Debiasing Using Matching Methods}\label{sec:debias}
\setcounter{equation}{0}
The matrix completion estimator of the ASF is generally biased. As we explained in Section \ref{sec:mc}, the bias comes from two sources: low-rank approximation bias and shrinkage bias. One could attempt to correct the shrinkage bias by shifting the singular values of $\widehat \boldsymbol{\Gamma}(x)$ upwards. However, inference results on the ASFs based on matrix completion are generally very difficult to obtain even if $ \boldsymbol{\Gamma}^{\infty}(x)$ is truly low rank. In our setting, the presence of the additional low-rank approximation bias makes this even more challenging. We instead discuss alternative estimators and show that they have significantly lower biases than the matrix completion estimators in the numerical simulations of Section \ref{sec:MC}.
To construct the estimators of $ \boldsymbol{\Gamma}^{\infty}(x)$, we start by extracting the factor structure of $\widehat \boldsymbol{\Gamma}(x)$ in \eqref{DefEstimator}. Let $ \widehat \boldsymbol{\lambda}_i(x)$ and $\widehat \boldsymbol{f}_t(x)$ be the $R \times 1$ vectors that satisfy
$$
\widehat \Gamma_{it}(x) = \widehat \boldsymbol{\lambda}_i(x)^\mathrm{\scriptscriptstyle{T}} \, \widehat \boldsymbol{f}_t(x) ,
$$
subject to the usual normalizations that $T^{-1} \sum_{t=1}^T \widehat \boldsymbol{f}_t(x) \, \widehat \boldsymbol{f}_t(x)^\mathrm{\scriptscriptstyle{T}} $ is the identity
matrix of size $R$ and
$N^{-1} \sum_{i=1}^N \widehat \boldsymbol{\lambda}_i(x) \, \widehat \boldsymbol{\lambda}_i(x)^\mathrm{\scriptscriptstyle{T}} $ is a diagonal matrix. Next, we apply a matching procedure to this factor structure. In its simplest version, we estimate each entry $ \boldsymbol{\Gamma}_{it}^{\infty}(x)$ such that $X_{it} \neq x$, by matching with the observation with $X_{js} = x$ that is the nearest neighbor in terms of the vectors $\widehat \boldsymbol{\lambda}_i(x)$ and $\widehat \boldsymbol{f}_t(x)$. In particular, $ \breve \Gamma_{it}(x) = Y_{i^{**}(i,t,x),t^{**}(i,t,x)}$ where $i^{**}(i,t,x) \in \mathbb{N}$ and $t^{**}(i,t,x) \in \mathbb{T}$ are a solution to the program
\begin{eqnarray*}
&\min_{j \in \mathbb{N}, s \in \mathbb{T}}& \left\| \widehat \boldsymbol{\lambda}_i(x) - \widehat \boldsymbol{\lambda}_j(x) \right\|^2 + \left\| \widehat \boldsymbol{f}_t(x) - \widehat \boldsymbol{f}_s(x) \right\|^2 \\
&\text{s.t.} & X_{js} =x.
\end{eqnarray*}
We also consider a two-way matching procedure that combines matching with a difference-in-differences approach. It consists of two steps:
\begin{itemize}
\item[(i)] For all $x \in \mathbb{X}$ and $(i,t) \in \mathbb{N} \times \mathbb{T}$ such that $X_{it} \neq x$, find the matches $i^*(i,t,x) \in \mathbb{N}$ and $t^*(i,t,x) \in \mathbb{T}$ that solve the program
\begin{eqnarray*}
&\min_{j \in \mathbb{N}, s \in \mathbb{T}}& \left\| \widehat \boldsymbol{\lambda}_i(x) - \widehat \boldsymbol{\lambda}_j(x) \right\|^2 + \left\| \widehat \boldsymbol{f}_t(x) - \widehat \boldsymbol{f}_s(x) \right\|^2 \\
&\text{s.t.} & X_{is} = X_{jt} = X_{js} =x.
\end{eqnarray*}
\item[(ii)] Estimate $ \Gamma_{it}(x)$ by
$$
\widetilde \Gamma_{it}(x) = Y_{i,t^*(i,t,x)} + Y_{i^*(i,t,x),t} - Y_{i^*(i,t,x),t^*(i,t,x)}.
$$
\end{itemize}
In other words, we find the match $(j,s)$ with $X_{js} = x$ that not only is the closest to $(i,t)$ in terms of the estimated factor structure, but also corresponds to a unit $j$ with $X_{jt} = x$ and a time period $s$ with $X_{is} = x$. Then, we estimate the counterfactual $ \Gamma_{it}(x)$ as a linear combination of $Y_{jt}$, $Y_{is}$ and $Y_{js}$.
The additional difference-in-differences step in the two-way procedure is useful to reduce bias. To see this, we can compare $ \widetilde \Gamma_{it}(x)$ with the simple matching estimator $\breve \Gamma_{it}(x)$.
Thus, abstracting from the estimation error in the factors and loadings,
\begin{multline*}
\mathrm{E}[\breve \Gamma_{it}(x) - \Gamma_{it}(x) \mid \boldsymbol{A}^N, \boldsymbol{B}^T, \boldsymbol{X}^{NT}] = m(x,\boldsymbol{A}_{i^{**}(i,t,x)},\boldsymbol{B}_{t^{**}(i,t,x)}) - m(x,\boldsymbol{A}_i,\boldsymbol{B}_t) \\
= \mathcal{O}_P(\|\boldsymbol{A}_{i^{**}(i,t,x)} - \boldsymbol{A}_i\| + \|\boldsymbol{B}_{t^{**}(i,t,x)} - \boldsymbol{B}_t\|),
\end{multline*}
by a first-order Taylor expansion of $(\boldsymbol{a}_{i}, \boldsymbol{b}_{t}) \mapsto m(x,\boldsymbol{a}_{i},\boldsymbol{b}_{t})$ around $(\boldsymbol{A}_i,\boldsymbol{B}_t)$; whereas
\begin{multline*}
\mathrm{E}[\widetilde \Gamma_{it}(x) - \Gamma_{it}(x) \mid \boldsymbol{A}^N, \boldsymbol{B}^T, \boldsymbol{X}^{NT} ] = m(x,\boldsymbol{A}_{i^{*}(i,t,x)},\boldsymbol{B}_{t^{*}(i,t,x)}) - m(x,\boldsymbol{A}_i,\boldsymbol{B}_t) \\
= \mathcal{O}_P(\|\boldsymbol{A}_{i^{*}(i,t,x)} - \boldsymbol{A}_i\|^2 + \|\boldsymbol{B}_{t^{*}(i,t,x)} - \boldsymbol{B}_t\|^2),
\end{multline*}
by a second-order Taylor expansion of $(\boldsymbol{a}_{i}, \boldsymbol{b}_{t}) \mapsto m(x,\boldsymbol{a}_{i},\boldsymbol{b}_{t})$ around $(\boldsymbol{A}_i,\boldsymbol{B}_t)$. The two-way matching removes the leading term of the Taylor expansion, reducing the bias of the matching by one order of magnitude because $i^{**}(i,t,x) \neq i$ or $t^{**}(i,t,x) \neq t$. On the other hand, $\|\boldsymbol{A}_{i^{*}(i,t,x)} - \boldsymbol{A}_i\| \geq \|\boldsymbol{A}_{i^{**}(i,t,x)} - \boldsymbol{A}_i\|$ and $\|\boldsymbol{B}_{t^{*}(i,t,x)} - \boldsymbol{B}_t\| \geq \|\boldsymbol{B}_{t^{**}(i,t,x)} - \boldsymbol{B}_t\|$ a.s. because the two-way procedure imposes the additional restrictions $X_{is} = X_{jt} =x$. Whether the first or second order bias dominates would generally be determined by the proportion of observations with $X_{js} =x$ and the distributions of $\boldsymbol{A}_i$ and $\boldsymbol{B}_t$. We provide a numerical comparison of the biases of the matching estimators in Section~\ref{sec:MC}.
We develop the theory for a debiased estimator that allows for multiple matches and estimated factors and loadings. Multiple matches are expected to reduce dispersion at the cost of increasing bias.
Let $\boldsymbol{\lambda}_i = \boldsymbol{\lambda}(x,\boldsymbol{A}_i)$ and $\boldsymbol{f}_t = \boldsymbol{f}(x,\boldsymbol{B}_t)$ be the transformations of $\boldsymbol{A}_i$
and $\boldsymbol{B}_t$
that are consistently estimated by $\widehat {\boldsymbol{\lambda}}_i$ and $\widehat {\boldsymbol{f}}_t$.\footnote{The matching method discussed here is also applicable to settings where the matching is based on variables other than the estimated factor structure. These include for example cross section and time series averages of the observable variables. See the appendix for a more general
treatment.}
We define
\begin{align*}
\mathbb{N}_i &= \Big\{ j \in \mathbb{N} \setminus \{i\} \, : \,\Big\| \widehat {\boldsymbol{\lambda}}_i - \widehat {\boldsymbol{\lambda}}_j \Big\| \leq \tau_{NT} \Big\} ,
&
\mathbb{T}_t &= \Big\{ s \in \mathbb{T} \setminus \{t\} \, : \, \Big\| \widehat {\boldsymbol{f}}_t - \widehat {\boldsymbol{f}}_s \Big\| \leq \upsilon_{NT} \Big\} ,
\end{align*}
for some bandwidth parameters $\tau_{NT}>0$ and $\upsilon_{NT} >0$.
The debiased estimator of $\mu(x)$ is then given by
$$
\widetilde \mu(x) = \frac 1 {NT} \sum_{i=1}^N \sum_{t=1}^T \widetilde Y_{it}(x),
$$
with
\begin{align}
\label{DefWidetildeY}
\widetilde Y_{it}(x) = \left\{ \begin{array}{lll}
Y_{it} \hfill \text{if $X_{it} =x$,}
\\[10pt]
\displaystyle \frac 1 {n_{it}} \sum_{j \in \mathbb{N}_i} \sum_{s \in \mathbb{T}_t} \mathbbm{1}\{X_{is}=X_{jt}=X_{js}=x\} ( Y_{is} + Y_{jt} - Y_{js} ) \\
\hfill \text{if $X_{it} \neq x$ and $n_{it}>0$,}
\\[10pt]
\frac 1 { n(x)} \sum_{(j,s) \in \mathbb{D}(x) } Y_{js} \hfill \text{if $n_{it}=0$,}
\end{array}
\right.
\end{align}
where $n_{it} := \sum_{j \in \mathbb{N}_i} \sum_{s \in \mathbb{T}_t} \mathbbm{1}\{X_{is}=X_{jt}=X_{js}=x\}$.
Here, for $X_{it} \neq x$, we construct the counterfactual $\widetilde Y_{it}(x) $
by averaging over all units $(j,s) \in \mathbb{N}_i \times \mathbb{T}_t$ that satisfy the constraint $X_{is}=X_{jt}=X_{js}=x$.
Notice that if $X_{it} \neq x$ and $n_{it}=0$, then we cannot construct a suitable counterfactual by that method. In that case we assign
$ \widetilde Y_{it}(x)$ the average of the observations with $X_{js} = x$ to make sure that $\widetilde \mu(x) $ is always well-defined, but our assumption below guarantees that this rarely happens.
This estimator has similar debiasing properties to the nearest neighbor described above, but it is more tractable theoretically because it varies more smoothly with respect to the factors and loadings.
Indeed, $\widetilde \mu(x)$ can be written as
$$\widetilde \mu(x) = \frac 1 {NT} \sum_{i=1}^N \sum_{t=1}^T \omega_{it} \, Y_{it},$$
where the weights $\omega_{it}$ are functions of
$\widehat {\boldsymbol{\lambda}}_j$ and $\widehat {\boldsymbol{f}}_s$ for all $j \in \mathbb{N}$ and $s \in \mathbb{T}$.
To show that $\widetilde \mu(x)$ is a consistent estimator of $\mu(x)$, we use the following assumption:
\begin{assumption}[Two-way Matching Estimator]
\label{ass:Debias}
There exists a sequence $\xi_{NT}>0$ such that $\xi_{NT} \to 0$ as $N,T \to \infty$, and
\begin{enumerate}[(a)]
\item $ \frac 1 {NT} \sum_{i=1}^N \sum_{t=1}^T \mathbbm{1}\left\{ X_{it} \neq x \, \& \, n_{it}=0 \right\} = O_P \left( \xi_{NT} \right)$.
\item $Y_{it}$ is uniformly bounded over $i,t,N,T$.
\item $Y_{it}$ is independent across both $i$ and $t$, conditional on $\boldsymbol{X}^{NT}$, $\boldsymbol{A}^N$, $\boldsymbol{B}^T$.
\item The function $(\boldsymbol{a}, \boldsymbol{b}) \mapsto m(x, \boldsymbol{a}, \boldsymbol{b})$ is at least twice continuously differentiable with uniformly bounded second derivatives.
\item There exists $c>0$ such that $\left\| \boldsymbol{a}_1 - \boldsymbol{a}_2 \right\| \leq c \left\| \boldsymbol{\lambda}(\boldsymbol{a}_1) - \boldsymbol{\lambda}(\boldsymbol{a}_2) \right\|$ for all $\boldsymbol{a}_1, \boldsymbol{a}_2 \in \mathbb{A}$, and $\left\| \boldsymbol{b}_1 - \boldsymbol{b}_2 \right\| \leq c \left\| \boldsymbol{f}(\boldsymbol{b}_1) - \boldsymbol{f}(\boldsymbol{b}_2) \right\|$ for all $\boldsymbol{b}_1, \boldsymbol{b}_2 \in \mathbb{B}$.
\item
$ \frac 1 {N} \sum_{i=1}^N \left( \left\| \widehat{\boldsymbol{\lambda}}_i - \boldsymbol{\lambda}_i \right\|^2
+ \max_{j \in \mathbb{N}_i} \left\| \widehat{\boldsymbol{\lambda}}_j - \boldsymbol{\lambda}_j \right\|^2 \right)= O_P \left( \xi_{NT} \right) $.
\\
$ \frac 1 {T} \sum_{t=1}^T \left( \left\| \widehat{\boldsymbol{f}}_t - \boldsymbol{f}_t \right\|^2
+ \max_{s \in \mathbb{T}_t} \left\| \widehat{\boldsymbol{f}}_s - \boldsymbol{f}_s \right\|^2 \right)= O_P \left( \xi_{NT} \right) $.
\item
$ \tau_{NT} ^2= O_P \left( \xi_{NT} \right) $ and $\upsilon_{NT}^2 = O_P \left( \xi_{NT} \right) $.
\item $
\frac 1 {NT} \sum_{i=1}^N \sum_{t=1}^T
\mathrm{E} \left[ \omega^2_{it} \,\big| \, \boldsymbol{X}^{NT}, \, \boldsymbol{A}^N, \, \boldsymbol{B}^T \right]
= O_P( NT \, \xi_{NT}^2).
$
\item
Let $\boldsymbol{Y}^{NT}_{-(i,t),-(j,s)}$ be the outcome matrix $\boldsymbol{Y}^{NT}$, but with $Y_{it}$ and $Y_{js}$ replace by zero (or some other non-random number),
and all other outcomes unchanged. We assume
\begin{multline*}
\frac 1 {(NT)^2} \sum_{i,j=1}^N \sum_{t,s=1}^T
\mathbbm{1}\left\{ (i,t) \neq (j,s) \right\}
\mathrm{E} \bigg[ \Big| \omega_{it}\left( \boldsymbol{Y}^{NT}_{-(i,t),-(j,s)} \right) \, \omega_{js}\left( \boldsymbol{Y}^{NT}_{-(i,t),-(j,s)} \right)
\\
- \omega_{it}(\boldsymbol{Y}^{NT}) \, \omega_{js}(\boldsymbol{Y}^{NT}) \Big|
\, \bigg| \, \boldsymbol{X}^{NT}, \, \boldsymbol{A}^N, \, \boldsymbol{B}^T \bigg]
= O_P\left( \xi_{NT}^2 \right) .
\end{multline*}
\end{enumerate}
\end{assumption}
\begin{remark}[Assumption~\ref{ass:Debias}]
\textnormal{
Part (i) guarantees that $X_{it} \neq x$ and $n_{it}=0$ only happens for a small fraction of observations $(i,t)$. We are therefore able to construct proper counterfactuals $ \widetilde Y_{it}(x)$ for most observations. Part (ii) is a boundedness condition that is standard in the matrix completion literature. Part (iii) is an independence condition that is convenient to simplify the derivations but can be generalized to weak correlation across both $i$ and $t$. We use part (iv) to bound the error terms of the Taylor expansions for the bias. Part (v) imposes an injectivity condition. The functions $\boldsymbol{a} \mapsto \boldsymbol{\lambda}(\boldsymbol{a})$ and $ \boldsymbol{b} \mapsto \boldsymbol{f}(\boldsymbol{b})$ need to be such that
$\boldsymbol{A}_i$ and $\boldsymbol{B}_t$ can be uniquely recovered from $\boldsymbol{\lambda}_i = \boldsymbol{\lambda}(\boldsymbol{A}_i)$ and $\boldsymbol{f}_t = \boldsymbol{f}(\boldsymbol{B}_t)$. A necessary condition is that the dimensions of $\boldsymbol{\lambda}_i$ and $\boldsymbol{f}_t$ are greater than or equal to the dimensions of $\boldsymbol{A}_i$ and $\boldsymbol{B}_t$, respectively. This holds in our factor structure approximation when let $R$ grow with the sample size, provided that the dimensions of $\boldsymbol{A}_i$ and $\boldsymbol{B}_t$ are fixed. Part (vi) holds if $ \widehat{\boldsymbol{\lambda}}_i - \boldsymbol{\lambda}_i$
and $\widehat{\boldsymbol{f}}_t - \boldsymbol{f}_t $ are of order $N^{-1/2}$ and $T^{-1/2}$.
We expect this assumption to be satisfied for rates $\xi_{NT} \gg \max(N^{-1},T^{-1})$.
The bandwidth parameters $ \tau_{NT}$ and $\upsilon_{NT}$ should not be chosen too large according to part (vii). For example, if we want to achieve a rate $ \xi_{NT} \ll \max(N^{-1/2},T^{-1/2})$, then
we need $ \tau_{NT} \ll \max(N^{-1/4},T^{-1/4})$ and $\upsilon_{NT} \ll \max(N^{-1/4},T^{-1/4})$.
Part (viii) requires that any given outcome $Y_{it}$ is not chosen too often with too high weight
in the construction of the counterfactuals $ \widetilde Y_{js}(x) $. Finally, part (ix) is a high-level assumption that could be justified by appropriate distributional assumptions
on $X_{it}$, $\boldsymbol{A}_i$, $\boldsymbol{B}_t$, and on the estimators $\widehat {\boldsymbol{\lambda}}_i$ and $\widehat {\boldsymbol{f}}_t$.
We prefer to present it as a high-level assumption,
because formally working out the distributional assumptions is quite cumbersome.
Intuitively, if $n_{it}$ is sufficiently large, then changing $\boldsymbol{Y}^{NT}$ to $ \boldsymbol{Y}^{NT}_{-(i,t),-(j,s)}$
should not change the constructions of the counterfactual $\widehat Y_{it}(x)$ very much.
If that is true for all $(i,t)$, then the weights $ \omega_{it}(\boldsymbol{Y}^{NT}) $
should be very close to the weights $ \omega_{it}\left( \boldsymbol{Y}^{NT}_{-(i,t),-(j,s)} \right) $
and the assumption is satisfied.}
\end{remark}
\begin{theorem}
\label{ConvergenceRateMatching}
Under Assumptions \ref{ass:model} and \ref{ass:Debias},
\begin{align*}
\widetilde \mu(x) - \mu(x) & = O_P \left( \xi_{NT} \right).
\end{align*}
\end{theorem}
As discussed in the above remark, one can achieve rates $ \xi_{NT} \ll \max(N^{-1/2},T^{-1/2})$
for sufficiently regular data generating processes, and if
the bandwidth parameters $ \tau_{NT}$ and $\upsilon_{NT}$ are chosen sufficiently small.
By contrast, the low-rank approximation bias in $\widehat \mu(x)$ will usually prevent us from achieving such a
convergence rate for $\widehat \mu(x)$. This finding is consistent with our Monte Carlo results in Section~\ref{sec:MC}, where
$ \widetilde \mu(x) $ is found to typically have much smaller bias than $\widehat \mu(x)$.
\section{Numerical Examples}\label{sec:examples}
\setcounter{equation}{0}
\subsection{Election day registration and voter turnout}
\label{sec:Election}
We illustrate the methods of the paper with an empirical application to the effect of allowing voter registration during the election day on voter turnout in the U.S. \citet{xu17}. Voting in the U.S. used to require registration prior to the election day in most states. Registration increased the cost of voting and was considered as one possible reason for low turnout rates. In response, some states implemented Election Day Registration (EDR) laws that allowed eligible voters to register on election day when they arrive at the polling stations. These laws were not passed by all the states, and there was variation in the time of adoption across states. Thus, they were enacted by Maine, Minnesota and Wisconsin in 1976; Wyoming, Indiana and New Hampshire in 1994, and Connecticut in 2012.
We use a dataset on the 24 presidential elections for 47 states between 1920 and 2012 collected by \citet{xu17}. It includes state-level information about the turnout rate, $Y_{it}$, measured as the total ballots counted divided by voting-age population in state $i$ at election $t$, and a treatment indicator for EDR, $X_{it}$, that equals one if the state $i$ has an EDR law enacted at election $t$. Following \citet{xu17}, we exclude North Dakota where registration was never needed, and Alaska and Hawaii that were not states until 1959. Since there are only 9 states that are ever treated and the treatment started in the 1976 election, we focus on effects on the treated at the elections between 1976 and 2012. We estimate average treatment effects and quantile treatment effects at multiple quantile indices.
Figure \ref{fig:pretrends} compares the average turnout of states that are ever treated with states that are never treated in elections prior to the first implementation of the EDR laws in 1976. It shows that ever treated states have higher turnout rates on average than never treated states without the EDR treatment. We consider several methods to deal with this likely nonrandom assignment of EDR to estimate the ATTs for each election after 1976. First, we do a naive comparison of means between treated and nontreated states in each election (Dmeans). Second, we consider a difference-in-differences method that uses the nontreated states as controls at each election (DiD). In particular, we estimate the effects from a linear regression with state effects and election effects interacted with a EDR indicator. This method yields the ATT for each election under a parallel trend assumption between treated and nontreated states.\footnote{The DiD model is a special case our model with additive effects. In this case, it imposes that there are only additive state and election effects that are the same for both treatment levels.} Third, we compute our estimator based on matrix completion methods without debiasing (MC) with additive state and election effects and the parameter $\rho$ such that the number of factors is $R = 6$. Fourth, we debias the MC estimates using the two-way matching method with 10 matches (TWM-10). Fifth, we consider the simple matching method with 5 matches (SM-5). We choose the number of matches roughly based on the numerical simulations of Section \ref{sec:MC}.
\begin{figure}
\begin{center}
\includegraphics[height=.8\textwidth,width=.8\textwidth,angle=0]{ptrends.pdf}\caption{Pretrends in turnout rate}\label{fig:pretrends}
\end{center}
\end{figure}
Figure \ref{fig:ATT} reports the estimates of the ATT of EDR at each election. The methods that account for possible nonrandom assignment of the EDR produce lower estimates of the effect than the naive comparison of means between treated and nontreated states. This finding agrees with the pre-EDR differences found in fig.~\ref{fig:pretrends}. MC, TWM-10 and SM-5 estimates are generally larger and more stable across elections than DiD estimates. According to TWM-10, EDR laws increase voter turnout between 5 and 9\% depending on the election. This effect is an economically significant relative to 55\%, the average turnout rate for states without EDR. The estimates of the election-aggregated ATTs are 10.71\%, 0.67\%, 7.35\%, 5.56\%, and 4.87\% for Dmeans, DiD, MC, TWM-10, and SM-3, respectively.
\begin{figure}[ht]
\begin{center}
\includegraphics[height=.8\textwidth,width=.8\textwidth,angle=0]{ATT_with_additive_effects.pdf}
\caption{Average treatment effect on treated}\label{fig:ATT}
\end{center}
\end{figure}
Figure \ref{fig:QTT} plots the estimates of the election-aggregated quantile treatment effect on the treated (QTT) of EDR as a function of the quantile index. We report estimates from four methods: a naive comparison of quantiles between treated and non-treated states (Dquantiles), our estimator based on matrix completion methods without debiasing (MC) with additive state and election effects and the parameter $\rho$ such that the number of factors is $R = 3$, two-way matching with 10 matches (TWM-10), and simple matching with 5 matches (SM-5). The QTT is the difference of the quantiles between the observed turnout for the treated observations and the corresponding potential turnout have they not been treated. The quantiles of the observed turnout are estimated using sample quantiles. The estimates of the quantiles of the potential outcomes are obtained by inverting the corresponding estimates of the distribution, which are obtained by our methods replacing $Y_{it}$ by the indicator $\mathbbm{1}(Y_{it} \leq y)$ and repeating the procedure over a grid of values of $y$ that includes the sample quantiles of observed turnout with indices $\{.10, .11, \ldots, .98\}$.\footnote{We rearrange the estimates of the distribution to guarantee that they are increasing with respect to $y$ \citet{CFG10}.} Here, we find that the effect of EDR is decreasing across the distribution of turnout and ranges between 10 and 0\% according to TWM-10. EDR is therefore more effective at the bottom of the voter turnout distribution. Comparing with the Dquantiles estimates, we find that the sign of the selection bias switches from positive to negative around the middle of the turnout distribution.
\begin{figure}[ht]
\begin{center}
\includegraphics[height=.8\textwidth,width=.8\textwidth,angle=0]{QTT_with_additive_effects.pdf}
\caption{Time-averaged QTT}\label{fig:QTT}
\end{center}
\end{figure}
\subsection{Monte Carlo simulations}
\label{sec:MC}
To evaluate the performance of our methods in a controlled synthetic environment, we generate potential outcomes from an additive linear model where
$$
Y_{it}(x) = x + g(A_i, B_t) + U_{it}(x), \quad x \in \{0,1\}, i \in \{1, \ldots, 30\}, t \in \{1, \ldots, 30\},
$$
${U}_{it}(x) \sim N(0,1/4)$ independently over $i$, $t$ and $x$, $A_i \sim U(0,1)$ independently over $i$, $B_t \sim U(0,1)$ independently over $t$, ${U}_{it}(x)$, $A_j$ and $B_s$ are independent for all $i$, $t$, $j$ and $s$, and $g$ is the Gaussian kernel, i.e.,
\begin{align*}
g(a,b) = \frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{(a-b)^2}{\sigma^2}\right).
\end{align*}
This design is similar to that used in \citet{bordenave2020detection}, with kernel function specification from the numerical simulations in \citet{gh13}.\footnote{We find similar results in a multiplicative model where $Y_{it}(x) = (1 + x) g(A_i, B_t) + U_{it}(x)$. We omit these results for the sake of brevity.} The parameter $\sigma$ controls the decay of the singular values of $g$ and can be calibrated to make sure the singular values decay slowly. Smaller values for $\sigma$ lead to greater dispersion in the kernel function $(a,b) \mapsto g(a,b)$ and a slower singular value decay, hence can be interpreted as a measure of smoothness.\footnote{Smoothness here is specifically related to numerical smoothness, i.e. variability in the function within close neighbourhoods of its arguments.} The assignment of $X_{it}$ that determines what potential outcomes are observed is similar to the election application. In particular, only observations for the first half of the units, $i \in \{1,\ldots, 15\}$, and the second half of the panel, $t \in \{ 15,\ldots,30 \}$, may be treated. For these observations, $X_{it}$ is related to the unobserved effects $(A_i,B_t)$ via $X_{it} = \mathbbm{1} \{g(A_i,B_t)\geq c\}$, where $c$ is a constant calibrated to $\Pr(g(A_i,B_t)\geq c) =.5$.
\begin{table}[ht]
\caption{Results for $\mu(0 \mid \{1\})$ }
\centering
\begin{tabular}{lcccc}
\hline\hline
&Bias & St. Dev. & RMSE \\%&&Bias & St. Dev. & RMSE\\
\hline
Dmeans & 0.59 & 0.02 & 0.59 \\%& & 0.59 & 0.02 & 0.59 \\
DiD & 0.70 & 0.03 & 0.70 \\%& & 0.69 & 0.03 & 0.70 \\
MC & 0.74 & 0.02 & 0.74 \\%& & 0.74 & 0.02 & 0.74 \\
TWM-1 & 0.03 & 0.14 & 0.14 \\%& & 0.02 & 0.13 & 0.14 \\
TWM-5 & 0.03 & 0.11 & 0.12 \\%& & 0.03 & 0.11 & 0.11 \\
TWM-10 & 0.04 & 0.10 & 0.11 \\%& & 0.04 & 0.10 & 0.11 \\
TWM-30 & 0.07 & 0.09 & 0.12 \\%& & 0.06 & 0.09 & 0.11 \\
SM-1 & 0.12 & 0.10 & 0.16 \\%& & 0.12 & 0.11 & 0.16 \\
SM-5 & 0.15 & 0.07 & 0.17 \\%& & 0.15 & 0.07 & 0.17\\
SM-10 & 0.19 & 0.06 & 0.20 \\%& & 0.19 & 0.06 & 0.20\\
SM-30 & 0.31 & 0.05 & 0.31\\% & & 0.31 & 0.05 & 0.31\\
\hline\hline
\multicolumn{4}{l}{\footnotesize{Notes: based on $1,000$ simulations}}
\end{tabular}
\label{table:CASF}
\end{table}
We apply similar methods to Section \ref{sec:Election} to estimate the CASFs $\mu_t(0 \mid \{1\}),$ $t \in \{ 15,\ldots,30 \}$, and $\mu(0 \mid \{1\})$ using the observed variables $X_{it}$ and $Y_{it} = Y_{it}(X_{it})$. Thus, we consider Dmeans, DiD, MC without additive effects and with the parameter $\rho$ such that $R=5$, and multiple versions of TWM and SM with the number of matches equal to $1$, $5$, $10$, and $30$. For each method, we compute the bias, standard deviation and rmse from $1,000$ simulations. Across the simulations, we redraw the values of $U_{it}(x)$ and hold $A_i$, $B_t$ and $X_{it}$ fixed. Table \ref{table:CASF} reports the results for the time-aggregated CASF, $\mu(0 \mid \{1\})$, and Figure \ref{fig:sim_linear} plots the results for the CASF, $\mu_t(0 \mid \{1\})$, as a function of $t$. The results show that Dmeans, DiD and MC are severely biased relative to their standard deviations. All the matching estimators reduce bias and rmse, despite of increasing dispersion. As one would expect, increasing the number of matches reduces the variability of the matching estimators but increases their biases. The number of matches that minimizes the rmse is larger for the TWM than for the SM. Overall, these small-sample findings agree with the asymptotic results of Sections \ref{sec:mc} and \ref{sec:debias}.
\begin{figure}[ht]
\centering
\includegraphics[width=0.32\linewidth, height = 6cm]{generic_biasOt.jpeg}
\includegraphics[width=0.32\linewidth, height = 6cm]{generic_stdevOt.jpeg}
\includegraphics[width=0.32\linewidth, height = 6cm]{generic_RMSEOt.jpeg}
\caption{Results for $t \mapsto \mu_t(0 \,|\, \{1\})$.}
\label{fig:sim_linear}
\end{figure}
\section*{Acknowledgements}
This paper was prepared for the Econometrics Journal Special Session on ``Econometrics of Panel Data'' at the Royal Economic Society 2019 Annual Conference in Warwick University. We thank the editor Jaap Abbring, two anonymous referees, Shuowen Chen, and the participants of this conference and the 25$^{\text{th}}$ International Panel Data Conference for comments. This research was
supported by the Economic and Social Research Council through the ESRC Centre for
Microdata Methods and Practice grant RES-589-28-0001, and by the
European Research Council grants ERC-2014-CoG-646917-ROMIA and
ERC-2018-CoG-819086-PANEDA.
\begin{thebibliography}{}
\bibitem[\protect\citeauthoryear{Amjad, Shah, and Shen}{Amjad
et~al.}{2018}]{amjad18}
Amjad, M., D.~Shah, and D.~Shen (2018).
\newblock Robust synthetic control.
\newblock {\em Journal of Machine Learning Research\/}~{\em 19\/}(22), 1--51.
\bibitem[\protect\citeauthoryear{Athey, Bayati, Doudchenko, Imbens, and
Khosravi}{Athey et~al.}{2017}]{athey17}
Athey, S., M.~Bayati, N.~Doudchenko, G.~Imbens, and K.~Khosravi (2017).
\newblock Matrix completion methods for causal panel data models.
\newblock arXiv eprint arXiv:1710.10251.
\bibitem[\protect\citeauthoryear{Auerbach}{Auerbach}{2019}]{auerbach2019identification}
Auerbach, E. (2019).
\newblock Identification and estimation of a partially linear regression model
using network data.
\newblock arXiv eprint arXiv:1903.09679.
\bibitem[\protect\citeauthoryear{Bai and Ng}{Bai and Ng}{2019a}]{bai19matrix}
Bai, J. and S.~Ng (2019a).
\newblock Matrix completion, counterfactuals, and factor analysis of missing
data.
\newblock arXiv eprint arXiv:1910.06677.
\bibitem[\protect\citeauthoryear{Bai and Ng}{Bai and Ng}{2019b}]{bai19joe}
Bai, J. and S.~Ng (2019b).
\newblock {Rank regularized estimation of approximate factor models}.
\newblock {\em Journal of Econometrics\/}~{\em 212\/}(1), 78--96.
\bibitem[\protect\citeauthoryear{Bai, Silverstein, and Yin}{Bai
et~al.}{1988}]{BaiSilvYin1988}
Bai, Z.~D., J.~W. Silverstein, and Y.~Q. Yin (1988).
\newblock A note on the largest eigenvalue of a large dimensional sample
covariance matrix.
\newblock {\em Journal of multivariate analysis\/}~{\em 26\/}(2), 166--168.
\bibitem[\protect\citeauthoryear{Beyhum and Gautier}{Beyhum and
Gautier}{2019}]{beyhum19}
Beyhum, J. and E.~Gautier (2019).
\newblock Square-root nuclear norm penalized estimator for panel data models
with approximately low-rank unobserved heterogeneity.
\newblock arXiv eprint arXiv:1904.09192.
\bibitem[\protect\citeauthoryear{Bordenave, Coste, and Nadakuditi}{Bordenave
et~al.}{2020}]{bordenave2020detection}
Bordenave, C., S.~Coste, and R.~R. Nadakuditi (2020).
\newblock Detection thresholds in very sparse matrix completion.
\newblock arXiv eprint arXiv:2005.06062.
\bibitem[\protect\citeauthoryear{Cai, Cand{\`e}s, and Shen}{Cai
et~al.}{2010}]{cai10}
Cai, J.-F., E.~J. Cand{\`e}s, and Z.~Shen (2010).
\newblock A singular value thresholding algorithm for matrix completion.
\newblock {\em SIAM Journal on optimization\/}~{\em 20\/}(4), 1956--1982.
\bibitem[\protect\citeauthoryear{Cand{\`e}s and Recht}{Cand{\`e}s and
Recht}{2009}]{candes09}
Cand{\`e}s, E.~J. and B.~Recht (2009).
\newblock Exact matrix completion via convex optimization.
\newblock {\em Foundations of Computational mathematics\/}~{\em 9\/}(6), 717.
\bibitem[\protect\citeauthoryear{{Cand\'es} and {Tao}}{{Cand\'es} and
{Tao}}{2010}]{candes10}
{Cand\'es}, E.~J. and T.~{Tao} (2010).
\newblock The power of convex relaxation: Near-optimal matrix completion.
\newblock {\em IEEE Transactions on Information Theory\/}~{\em 56\/}(5),
2053--2080.
\bibitem[\protect\citeauthoryear{Chamberlain}{Chamberlain}{1982}]{Chamberlain82}
Chamberlain, G. (1982).
\newblock Multivariate regression models for panel data.
\newblock {\em Journal of Econometrics\/}~{\em 18\/}(1), 5--46.
\bibitem[\protect\citeauthoryear{Chan and Kwok}{Chan and Kwok}{2020}]{chan20}
Chan, M.~K. and S.~Kwok (2020, March).
\newblock {The PCDID Approach: Difference-in-Differences when Trends are
Potentially Unparallel and Stochastic}.
\newblock Working Papers 2020-03, University of Sydney, School of Economics.
\bibitem[\protect\citeauthoryear{Chatterjee et~al.}{Chatterjee
et~al.}{2015}]{chatterjee2015}
Chatterjee, S. et~al. (2015).
\newblock Matrix estimation by universal singular value thresholding.
\newblock {\em Annals of Statistics\/}~{\em 43\/}(1), 177--214.
\bibitem[\protect\citeauthoryear{Chen, Fern{\'a}ndez-Val, and Weidner}{Chen
et~al.}{2021}]{chen2020nonlinear}
Chen, M., I.~Fern{\'a}ndez-Val, and M.~Weidner (2021).
\newblock Nonlinear factor models for network and panel data.
\newblock {\em Journal of Econometrics\/}~{\em 220\/}(2), 296--324.
\bibitem[\protect\citeauthoryear{Chernozhukov, Fern{\'a}ndez-Val, and
Galichon}{Chernozhukov et~al.}{2010}]{CFG10}
Chernozhukov, V., I.~Fern{\'a}ndez-Val, and A.~Galichon (2010).
\newblock Quantile and probability curves without crossing.
\newblock {\em Econometrica\/}~{\em 78\/}(3), 1093--1125.
\bibitem[\protect\citeauthoryear{Chernozhukov, Fern{\'a}ndez-Val, Hahn, and
Newey}{Chernozhukov et~al.}{2013}]{CFHN13}
Chernozhukov, V., I.~Fern{\'a}ndez-Val, J.~Hahn, and W.~Newey (2013).
\newblock Average and quantile effects in nonseparable panel models.
\newblock {\em Econometrica\/}~{\em 81\/}(2), 535--580.
\bibitem[\protect\citeauthoryear{Chernozhukov, Hansen, Liao, and
Zhu}{Chernozhukov et~al.}{2018}]{chernozhukov18}
Chernozhukov, V., C.~Hansen, Y.~Liao, and Y.~Zhu (2018).
\newblock Inference for heterogeneous effects using low-rank estimation of
factor slopes.
\newblock arXiv eprint arXiv:1812.08089.
\bibitem[\protect\citeauthoryear{Chernozhukov, Hansen, Liao, and
Zhu}{Chernozhukov et~al.}{2020}]{chlz20}
Chernozhukov, V., C.~Hansen, Y.~Liao, and Y.~Zhu (2020).
\newblock Inference for low-rank models.
\newblock Working Paper.
\bibitem[\protect\citeauthoryear{Dzemski}{Dzemski}{2019}]{dzemski2019empirical}
Dzemski, A. (2019).
\newblock An empirical model of dyadic link formation in a network with
unobserved heterogeneity.
\newblock {\em Review of Economics and Statistics\/}~{\em 101\/}(5), 763--776.
\bibitem[\protect\citeauthoryear{Evdokimov}{Evdokimov}{2010}]{evdokimov2010identification}
Evdokimov, K. (2010).
\newblock Identification and estimation of a nonparametric panel data model
with unobserved heterogeneity.
\newblock Working Paper.
\bibitem[\protect\citeauthoryear{Fazel}{Fazel}{2003}]{fazel03}
Fazel, S.~M. (2003).
\newblock Matrix rank minimization with applications.
\newblock Elec Eng Dept Stanford University 54, 1-130
\bibitem[\protect\citeauthoryear{Freyberger}{Freyberger}{2017}]{freyberger18}
Freyberger, J. (2017, 09).
\newblock {Non-parametric Panel Data Models with Interactive Fixed Effects}.
\newblock {\em The Review of Economic Studies\/}~{\em 85\/}(3), 1824--1851.
\bibitem[\protect\citeauthoryear{Gao, Lu, Zhou, et~al.}{Gao
et~al.}{2015}]{gao2015rate}
Gao, C., Y.~Lu, H.~H. Zhou, et~al. (2015).
\newblock Rate-optimal graphon estimation.
\newblock {\em The Annals of Statistics\/}~{\em 43\/}(6), 2624--2652.
\bibitem[\protect\citeauthoryear{Geman}{Geman}{1980}]{Geman1980}
Geman, S. (1980, April).
\newblock A limit theorem for the norm of random matrices.
\newblock {\em Annals of Probability\/}~{\em 8\/}(2), 252--261.
\bibitem[\protect\citeauthoryear{Gobillon and Magnac}{Gobillon and
Magnac}{2016}]{gobillon16}
Gobillon, L. and T.~Magnac (2016).
\newblock Regional policy evaluation: Interactive fixed effects and synthetic
controls.
\newblock {\em The Review of Economics and Statistics\/}~{\em 98\/}(3),
535--551.
\bibitem[\protect\citeauthoryear{Graham}{Graham}{2017}]{graham2017econometric}
Graham, B.~S. (2017).
\newblock An econometric model of network formation with degree heterogeneity.
\newblock {\em Econometrica\/}~{\em 85\/}(4), 1033--1063.
\bibitem[\protect\citeauthoryear{Graham and Powell}{Graham and
Powell}{2012}]{GrahamPowell2012}
Graham, B.~S. and J.~L. Powell (2012).
\newblock Identification and estimation of average partial effects in
ñirregularî correlated random coefficient panel data models.
\newblock {\em Econometrica\/}~{\em 80\/}(5), 2105--2152.
\bibitem[\protect\citeauthoryear{Griebel and Harbrecht}{Griebel and
Harbrecht}{2013}]{gh13}
Griebel, M. and H.~Harbrecht (2013, 05).
\newblock {Approximation of bi-variate functions: singular value decomposition
versus sparse grids}.
\newblock {\em IMA Journal of Numerical Analysis\/}~{\em 34\/}(1), 28--54.
\bibitem[\protect\citeauthoryear{Hoderlein and White}{Hoderlein and
White}{2012}]{HoderleinWhite2012}
Hoderlein, S. and H.~White (2012).
\newblock Nonparametric identification in nonseparable panel data models with
generalized fixed effects.
\newblock {\em Journal of Econometrics\/}~{\em 168\/}(2), 300--314.
\bibitem[\protect\citeauthoryear{Holland, Laskey, and Leinhardt}{Holland
et~al.}{1983}]{holland1983stochastic}
Holland, P.~W., K.~B. Laskey, and S.~Leinhardt (1983).
\newblock Stochastic blockmodels: First steps.
\newblock {\em Social networks\/}~{\em 5\/}(2), 109--137.
\bibitem[\protect\citeauthoryear{Honor{\'e}}{Honor{\'e}}{1992}]{Honore1992}
Honor{\'e}, B. (1992).
\newblock {Trimmed LAD and least squares estimation of truncated and censored
regression models with fixed effects}.
\newblock {\em Econometrica\/}~{\em 60\/}(3), 533--565.
\bibitem[\protect\citeauthoryear{Hsiao, Steve~Ching, and Ki~Wan}{Hsiao
et~al.}{2012}]{hsiao12}
Hsiao, C., H.~Steve~Ching, and S.~Ki~Wan (2012).
\newblock A panel data approach for program evaluation: Measuring the benefits
of political and economic integration of hong kong with mainland china.
\newblock {\em Journal of Applied Econometrics\/}~{\em 27\/}(5), 705--740.
\bibitem[\protect\citeauthoryear{Imai and Kim}{Imai and
Kim}{2019}]{imai2019use}
Imai, K. and I.~S. Kim (2019).
\newblock On the use of two-way fixed effects regression models for causal
inference with panel data.
\newblock Forthcoming in {\em Political Analysis}.
\bibitem[\protect\citeauthoryear{Kim and Oka}{Kim and Oka}{2014}]{KimOka2014}
Kim, D. and T.~Oka (2014).
\newblock Divorce law reforms and divorce rates in the usa: an interactive
fixed-effects approach.
\newblock {\em Journal of Applied Econometrics\/}~{\em 29\/}(2), 231--245.
\bibitem[\protect\citeauthoryear{Klopp et~al.}{Klopp
et~al.}{2014}]{klopp2014noisy}
Klopp, O. et~al. (2014).
\newblock Noisy low-rank matrix completion with general sampling distribution.
\newblock {\em Bernoulli\/}~{\em 20\/}(1), 282--303.
\bibitem[\protect\citeauthoryear{Lata\l{}a}{Lata\l{}a}{2005}]{latala05}
Lata\l{}a, R. (2005).
\newblock Some estimates of norms of random matrices.
\newblock {\em Proceedings of the American Mathematical Society\/}~{\em
133\/}(5), 1273--1282.
\bibitem[\protect\citeauthoryear{Li}{Li}{2018}]{li18}
Li, K. (2018).
\newblock Inference for factor model based average treatment effects.
\newblock Available at SSRN 3112775.
\bibitem[\protect\citeauthoryear{Li and Bell}{Li and Bell}{2017}]{li17}
Li, K.~T. and D.~R. Bell (2017).
\newblock Estimation of average treatment effects with panel data: Asymptotic
theory and implementation.
\newblock {\em Journal of Econometrics\/}~{\em 197\/}(1), 65 -- 75.
\bibitem[\protect\citeauthoryear{Li, Shah, Song, and Yu}{Li
et~al.}{2017}]{li17b}
Li, Y., D.~Shah, D.~Song, and C.~L. Yu (2017).
\newblock Nearest neighbors for matrix estimation interpreted as blind
regression for latent variable model.
\newblock arXiv eprint arXiv:1705.04867.
\bibitem[\protect\citeauthoryear{Ma, Goldfarb, and Chen}{Ma
et~al.}{2011}]{ma11}
Ma, S., D.~Goldfarb, and L.~Chen (2011).
\newblock Fixed point and bregman iterative methods for matrix rank
minimization.
\newblock {\em Mathematical Programming\/}~{\em 128\/}(1-2), 321--353.
\bibitem[\protect\citeauthoryear{Manski}{Manski}{1987}]{Manski1987}
Manski, C. (1987).
\newblock {Semiparametric analysis of random effects linear models from binary
panel data}.
\newblock {\em Econometrica\/}~{\em 55\/}(2), 357--362.
\bibitem[\protect\citeauthoryear{Mazumder, Hastie, and Tibshirani}{Mazumder
et~al.}{2010}]{mazumder10}
Mazumder, R., T.~Hastie, and R.~Tibshirani (2010).
\newblock Spectral regularization algorithms for learning large incomplete
matrices.
\newblock {\em Journal of Machine Learning Research\/}~{\em 11\/}(80),
2287--2322.
\bibitem[\protect\citeauthoryear{Menzel}{Menzel}{2018}]{menzel2018bootstrap}
Menzel, K. (2018).
\newblock Bootstrap with cluster-dependence in two or more dimensions.
\newblock arXiv preprint arXiv:1703.03043
\bibitem[\protect\citeauthoryear{Moon and Weidner}{Moon and
Weidner}{2017}]{MoonWeidner2017}
Moon, H.~R. and M.~Weidner (2017).
\newblock Dynamic linear panel regression models with interactive fixed
effects.
\newblock {\em Econometric Theory\/}~{\em 33\/}(1), 158--195.
\bibitem[\protect\citeauthoryear{Moon and Weidner}{Moon and
Weidner}{2018}]{moon18}
Moon, H.~R. and M.~Weidner (2018).
\newblock Nuclear norm regularized estimation of panel regression models.
\newblock arXiv eprints arXiv:1810.10987.
\bibitem[\protect\citeauthoryear{Negahban and Wainwright}{Negahban and
Wainwright}{2012}]{negahban2012restricted}
Negahban, S. and M.~J. Wainwright (2012).
\newblock Restricted strong convexity and weighted matrix completion: Optimal
bounds with noise.
\newblock {\em The Journal of Machine Learning Research\/}~{\em 13\/}(1),
1665--1697.
\bibitem[\protect\citeauthoryear{{Orbanz} and {Roy}}{{Orbanz} and
{Roy}}{2015}]{or15}
{Orbanz}, P. and D.~M. {Roy} (2015).
\newblock Bayesian models of graphs, arrays and other exchangeable random
structures.
\newblock {\em IEEE Transactions on Pattern Analysis and Machine
Intelligence\/}~{\em 37\/}(2), 437--461.
\bibitem[\protect\citeauthoryear{Rennie and Srebro}{Rennie and
Srebro}{2005}]{rennie05}
Rennie, J. D.~M. and N.~Srebro (2005).
\newblock Fast maximum margin matrix factorization for collaborative
prediction.
\newblock In {\em Proceedings of the 22nd International Conference on Machine
Learning}, ICML Õ05, New York, NY, USA, pp.\ 713--719. Association for
Computing Machinery.
\bibitem[\protect\citeauthoryear{Silverstein}{Silverstein}{1989}]{Silverstein1989}
Silverstein, J.~W. (1989).
\newblock On the eigenvectors of large dimensional sample covariance matrices.
\newblock {\em Journal of Multivariate Analysis\/}~{\em 30\/}(1), 1--16.
\bibitem[\protect\citeauthoryear{Srebro and Jaakkola}{Srebro and
Jaakkola}{2003}]{srebro03}
Srebro, N. and T.~Jaakkola (2003).
\newblock Weighted low-rank approximations.
\newblock In {\em Proceedings of the 20th International Conference on Machine
Learning (ICML-03)}, pp.\ 720--727.
\bibitem[\protect\citeauthoryear{Wolfe and Olhede}{Wolfe and
Olhede}{2013}]{wolfe2013nonparametric}
Wolfe, P.~J. and S.~C. Olhede (2013).
\newblock Nonparametric graphon estimation.
\newblock arXiv eprint arXiv:1309.5936.
\bibitem[\protect\citeauthoryear{Xiong and Pelger}{Xiong and
Pelger}{2019}]{xiong19}
Xiong, R. and M.~Pelger (2019).
\newblock Large dimensional latent factor modeling with missing observations
and applications to causal inference.
\newblock arXiv eprint arXiv:1910.08273.
\bibitem[\protect\citeauthoryear{Xu, Massouli, and Lelarge}{Xu
et~al.}{2014}]{xu14}
Xu, J., L.~Massouli, and M.~Lelarge (2014).
\newblock Edge label inference in generalized stochastic block models: from
spectral theory to impossibility results.
\newblock arXiv eprint arXiv:1406.6897.
\bibitem[\protect\citeauthoryear{Xu}{Xu}{2017}]{xu17}
Xu, Y. (2017).
\newblock Generalized synthetic control method: Causal inference with
interactive fixed effects models.
\newblock {\em Political Analysis\/}~{\em 25\/}(1), 57--76.
\bibitem[\protect\citeauthoryear{Yin, Bai, and Krishnaiah}{Yin
et~al.}{1988}]{BaiKrishYin1988}
Yin, Y.~Q., Z.~D. Bai, and P.~Krishnaiah (1988).
\newblock On the limit of the largest eigenvalue of the large-dimensional
sample covariance matrix.
\newblock {\em Probability Theory Related Fields\/}~{\em 78}, 509--521.
\bibitem[\protect\citeauthoryear{Zeleneev}{Zeleneev}{2020}]{zeleneev2020identification}
Zeleneev, A. (2020).
\newblock Identification and estimation of network models with nonparametric
unobserved heterogeneity.
\newblock Working Paper.
\end{thebibliography}