Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
100,948 characters · 12 sections · 77 citation commands
Low-Rank Approximations of Nonseparable Panel Models
\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. Chamberlain82, Manski1987, Honore1992, evdokimov2010identification, GrahamPowell2012, HoderleinWhite2012 and 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 athey17 and 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 srebro03 and 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 restricted strong convexity assumption (see e.g.\ 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 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 imai2019use to estimate causal effects from panel data and the blind regression of li17b for matrix completion. The difference with these proposals is in the information used to match the observations. imai2019use use the treatment variable and 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 chatterjee2015. 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, 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 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) describes the model and effects of interest. Section (ref) 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). Section (ref) reports the results of the numerical examples. All the proofs of the theoretical results are gathered in the Appendix.
\setcounter{equation}{0}
Throughout this paper we consider the following nonseparable and nonparametric panel data model:
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. xu14, chatterjee2015, or15, and 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) is the unit and time homogeneity in (ref). 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 Chamberlain82,Manski1987,Honore1992,evdokimov2010identification,GrahamPowell2012,HoderleinWhite2012,CFHN13. It implies that time is randomly assigned, conditional on covariates and unobserved effects. The additional restrictions in (ref) 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 (ref), 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, (ref) and (ref) impose that
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 (ref), 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), $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 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 (ref). For example, it can include location and scale time effects. Second, CFHN13 impose that $Y_{it}$ and $\boldsymbol{X}_{it}$ are identically distributed across $i$, which is stronger than the unit homogeneity in (ref). 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 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) 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
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)
and the conditional average structural functions (CASFs)
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 (ref) and (ref) 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). The purpose is to show that Assumption (ref) 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). 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.
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. graham2017econometric, dzemski2019empirical, chen2020nonlinear, zeleneev2020identification) or uses stochastic blockmodels or graphon models, which typically ignore the effect of covariates (e.g. holland1983stochastic, wolfe2013nonparametric, gao2015rate, auerbach2019identification). Our methods of estimating non-parametric models with two-way heterogeneity may therefore also be of interest in a network context.
\setcounter{equation}{0}
A natural starting point to estimate the effects in (ref) and (ref) 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
under the rank similarity in (ref) and Assumption (ref). 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
where the function $m$ does not vary with $i$ and $t$, by implication (ref) of Assumption (ref). We show next how this function can be approximated and estimated using a low-rank factor structure.
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.
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}) $.
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 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) 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 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,
The first term is the approximation and the second term is the approximation error. Under Assumption (ref), $$ \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
by Assumptions (ref) and (ref). Hence, $\zeta_R(x,\boldsymbol{A}_i,\boldsymbol{B}_t)$ converges in mean square to zero at a polynomial rate with $R$.
Combining (ref) and (ref), we obtain the approximate factor model
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) \[ 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 (ref) 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
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.
We start by expressing the program (ref) 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
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)$ rennie05,candes09,candes10. This connection was previously noticed by athey17 and 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 srebro03. Following 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 (ref) 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
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 mazumder10.
Let $\widehat \boldsymbol{\Gamma}(x)$ be a solution to (ref) 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 cai10,ma11,bai19joe. We propose two matching procedures to debias the estimator in Section (ref).
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
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
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.
For the purpose of showing Lemma (ref) and Theorem (ref) we could alternatively replace Assumption (ref) by the two high-level conditions:
where again $\|\cdot\|_\infty$ denotes the spectral norm. The first of those conditions is implied by Assumption (ref) through application of the weak law of large numbers, while the second follows, for example, by the spectral norm inequality in 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)$.
A necessary condition for the existence of the sequence $\rho = \rho_{NT}$ in Lemma (ref) 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) does not restrict the assignment process that determines $\boldsymbol{X}^{NT}$. Notice also that Lemma (ref) does not require Assumption (ref) 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) 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) 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
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 restricted strong convexity to derive (ref). 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
See Theorem 1 in negahban2012restricted, Lemma 12 in klopp2014noisy, and Lemma 3 in athey17. Thus, if $M_{it} = \widehat \Gamma_{it}(x) - \Gamma^{\infty}_{it}(x)$ and $\boldsymbol{X}^{NT}$ satisfy restricted strong convexity, then (ref) would follow from Lemma (ref).
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) is equal to $\mu(x)$ defined in (ref). 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 (ref), 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), $W_{it}(x) = n(1)^{-1} X_{it}$.
To interpret the conditions in Theorem (ref), 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 (ref) 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, Geman1980, Silverstein1989, BaiSilvYin1988, 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 latala05, see the proof of Lemma (ref) 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 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 (ref). 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 (ref). 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 (ref) 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 cai10,ma11,bai19joe. Because of these biases, we explore alternative estimates for $\mu(x)$ in Section (ref).
As we mentioned in Section (ref), 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 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 (ref) by
chernozhukov18, moon18 and 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
\setcounter{equation}{0}
The matrix completion estimator of the ASF is generally biased. As we explained in Section (ref), 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).
To construct the estimators of $ \boldsymbol{\Gamma}^{\infty}(x)$, we start by extracting the factor structure of $\widehat \boldsymbol{\Gamma}(x)$ in (ref). 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
We also consider a two-way matching procedure that combines matching with a difference-in-differences approach. It consists of two steps:
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,
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
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).
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
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
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:
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), where $ \widetilde \mu(x) $ is found to typically have much smaller bias than $\widehat \mu(x)$.
\setcounter{equation}{0}
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. 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 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 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) 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).
Figure (ref) 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). 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.
Figure (ref) 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$ 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.
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.,
This design is similar to that used in bordenave2020detection, with kernel function specification from the numerical simulations in 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$.
We apply similar methods to Section (ref) 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) reports the results for the time-aggregated CASF, $\mu(0 \mid \{1\})$, and Figure (ref) 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) and (ref).
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.