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.
87,395 characters · 23 sections · 65 citation commands
Linear Regressions with Combined Data
It is often impossible to run the ideal regressions one would like to consider. A common reason for this is that the outcome $Y$ and covariates $X$ of interest are not observed in the same dataset. For instance, in many intergenerational studies (e.g., intergenerational income or wealth mobility), one cannot link parents' and children's outcomes. Even if the outcome and covariates of interest do appear in the same dataset, key control variables are often missing. For instance, when measuring the wage returns to education, one may wish to control for a measure of cognitive skills, but such measure may not be available in the main labor market dataset, even though it appears in another source.
To better explain our contribution in this context, we first detail our setup. We assume that $X$ includes two sets of covariates: “outside” regressors $X_o$, which only appear in a separate dataset from that including the outcome $Y$, and “common” regressors $X_c$, which appear in both datasets. We also consider auxiliary variables, $W_a$, which researchers do not seek to include in the regression but that also appear in both datasets. These types of auxiliary variables are often available to empirical researchers. For instance, if a common variable is a proxy for a variable of interest $X_o$, or a so-called “bad control”, it seems preferable to focus on the regression of $Y$ on $X_o$, without controlling for that variable. We denote the set of common variables, included or not in the regression, by $W$, so that $W=(X_c', W_a')'$.
In this context, empirical researchers have traditionally relied on imputation methods. The most common, which corresponds to two-sample two-stage least squares (TSTSLS), consists in first predicting $X_o$ by $W_a$ in the “$X_o$ dataset”, and then using this prediction in the “$Y$ dataset”. One must recognize, however, that these approaches implicitly rely on exclusion restrictions and can therefore be sensitive to their violation. For instance, imputation based on TSTSLS requires the coefficient of $W_a$ in the infeasible regression of $Y$ on $X$ and $W_a$ to be 0. The goal of our paper is to study identification, estimation and inference on the regression coefficients without such exclusion restrictions.
In the absence of common variables, we obtain sharp bounds on each regression coefficient by applying the Frisch-Waugh-Lovell theorem together with the Cambanis-Simons-Stout inequality. Our contribution here is to show that this approach delivers sharp bounds, which is not obvious when $X_o$ is multivariate. We then extend this result to account for common variables. Sharp bounds still take a simple form in this case. Moreover, by leveraging the variation in $W$, one may be able to identify the sign of regression coefficients, and even obtain point identification in special cases. Importantly, we also show that it is possible to reject the exclusion restriction underlying the imputation based on TSTSLS.
Based on the identification results, we next turn to the estimation of the sharp bounds and inference on the regression coefficients. We propose simple plug-in estimators and establish their asymptotic normality. To do so, we build on results on $L$-statistics and on the statistical optimal transport literature. Our proof relies in particular on a refinement of the convergence rates in fournier2015rate, which we obtain by extending a result of boucheron2015. We also provide simple confidence intervals for the regression coefficients and establish their asymptotic validity. Simulation results indicate that our inference method performs well in finite samples, while being implementable at a modest computational cost.
Finally, we apply our methodology to two different contexts of data combination. In our first application, we revisit racial disparities in U.S. patent approval and show that conclusions about such disparities hinge on the validity of the exclusion restrictions underlying the TSTSLS estimation strategy. When relaxing these restrictions, our bounds are generally wide, pointing to the lack of robustness of the conclusions one would reach using TSTSLS results. In our second application, we evaluate the relationship between students' patience and risk-taking and educational performance across countries. In contrast to our first application, our bounds are informative and, for some of the specifications, exclude the TSTSLS estimates. Taken together, these applications highlight the limitations and, in some cases, misleading nature of the TSTSLS estimates in data combination environments. The partial identification approach we propose in this paper provides a transparent and tractable way to assess the sensitivity of empirical conclusions to the exclusion restrictions underlying the TSTSLS estimates.
\paragraph{Related literature.} To our knowledge, the first paper that considered our problem is pacini2019two. We extend his work in three important dimensions. First, we study the case where some of the common variables are not used as common regressors ($W\ne X_c$). We expect this to be prevalent in practice and we show that it can drastically reduce the identified sets. Second, we show that his bounds are not sharp when $X_o$ is multivariate, and that the difference with the sharp bounds can be substantial. Finally and importantly, we consider estimation and inference. hwang2022bounding also relates closely to our work. While she maintains the restriction that $W= X_c$, she also considers the case where some regressors are only available in the $Y$ dataset, which we do not study here.\footnote{There are still other data combination cases that we do not consider here. kitagawa2023 consider a setup where one observes $(Y,X_1, X_c)$ in one dataset and $(Y, X_2, X_c)$ in another. Yet another possibility, considered by moon2024partial when $X_1, X_2, X_c$ has finite support, is to observe $(Y, X_c)$, $(X_1, X_c)$ and $(X_2, X_c)$ separately.} A third related study is FPPS25. This paper complements ours by studying identification in a more general setup. In particular, they derive sharper bounds than hwang2022bounding in cases where some regressors are observed only in the $Y$ dataset. Their analysis, however, is limited to identification, whereas estimation, inference and empirical applicability are central to our paper.
Our paper is also related to our own previous work d2024partially, in which we consider a similar data combination environment. There are, however, important differences between the two. First, we did not consider previously how auxiliary variables affect identification. Second, we imposed a partially linear model, namely $E[Y|X]=X_o'\beta_o + f(X_c)$. This leads to potentially tighter bounds, but one may be reluctant to improve bounds using such restrictions. Third, for estimation we had to focus on the case for which $X_c$ had finite support, whereas no such assumption is necessary here. Finally, from a technical viewpoint, the restriction on the conditional expectation implies that we relied on entirely different optimal transport results, both for identification and inference.
At a broader level, our paper belongs to a very active literature on data combination problems in econometrics and statistics. See, in particular, RM2007 for a survey of this literature and contributions by fan2014identifying, FSS16, BLL16, BFM24, meango2025combining and, in the context of experimental data under a surrogacy assumption, ACI20, ACIK24, and rambachan2024program. Several of these papers impose restrictions that entail point identification. Following the seminal contribution of cross2002regressions and subsequent article by MP06, our aim is to obtain bounds on parameters of interest under weak restrictions. An important distinction between our work and these last two papers is that we consider different parameters: the best linear parameter in our case versus conditional expectation in theirs. Also, they do not consider the use of auxiliary variables ($W_a$).
From a technical viewpoint, our first identification result can be seen as an extension of the Cambanis-Simons-Stout inequality, see cambanis1976inequalities and, e.g., fan2014identifying (fan2014identifying, FSS16) for an application to data combination problems. Our asymptotic normality result relates to the asymptotic normality of the so-called Wasserstein-2 distance of empirical measures, recently studied in the statistical literature del2019central, berthet2020. Notably, up to a mild strengthening of moment conditions (from order 4 to $4+\varepsilon$ for some $\varepsilon>0$), our result implies asymptotic normality of the Wasserstein-2 distance under weaker conditions than those in berthet2020.
Finally, our paper also speaks to a large and growing empirical literature that deals with data combination problems similar to the one considered here. One important example is voting: given the anonymity of ballots, researchers typically regress average votes on average voter characteristics (e.g., income, hours watching Fox News per week) at the county level martin2017bias. This approach implicitly relies on a TSTSLS strategy, where counties play the role of $W_a$, thereby imposing an exclusion restriction. Another leading example is intergenerational income mobility, which often faces the unavailability of linked income data across generations and similarly relies on exclusion restrictions SantavirtaStuhler22. Data combination issues are also pervasive in consumption research, where income and consumption are often measured in separate datasets CLP22. Similar data combination problems frequently arise in various other subfields, including the economics of education and returns to skill estimation PP16,garcia2016life,hanushek2020culture, health Manski18 and labor ACI20. Finally, gaps in science and innovation by race or gender provide another relevant example, as illustrated in our first application below.
The methods we devise in this paper are broadly applicable in these different contexts, allowing empirical researchers to relax the exclusion restrictions that are typically maintained to achieve point identification. By applying our method to racial disparities in patent approval Dossi23 and the effect of preferences on skill differences hanushek2020culture, our paper also adds to the empirical literature on these questions.
\paragraph{Outline.} Section (ref) introduces the setup and discusses three broad cases for which our analysis is relevant. Section (ref) presents our identification results. Section (ref) develops estimators of the sharp bounds, establishes their asymptotic normality and develops inference on the regression coefficients. Section (ref) examines the finite sample properties of our estimators and confidence intervals through Monte Carlo simulations. We provide in Section (ref) two applications, to racial disparities in patent approval and the effect of preferences on skill differences. Finally, Section (ref) concludes. The appendix includes in particular a discussion of the sharpness of the bounds of pacini2019two and gathers all the proofs of our identification results; the proofs of our inference results appear in the online appendix. Finally, our method can be implemented using our companion R package, RegCombinBLP.\footnote{ This package is available on GitHub at \href{https://github.com/cgaillac/RegCombinBLP}{https://github.com/cgaillac/RegCombinBLP} with a user-friendly guide on how to use it.}
We seek to identify the best linear predictor $EL(Y|X)$ of $Y$ by $X\in\mathbb R^p$, with $X=(X'_o,X'_c)'$. To this end, we assume to have access to two separate datasets that cannot be matched. The first one includes $(Y, W')$, whereas the second one includes $(X_o',W')$; here $W=(W_a',X_c')\in \mathbb R^q$. We call $X_o$ the “outside regressors”, $X_c$ the “common regressors”, $W$ the “common variables” and $W_a$ the “auxiliary variables”. The latter are variables that the researcher does not want to include in the regression of interest, but that may still help for identification since they are included in both datasets. Importantly, they should not be seen as instruments, in the sense that we do not impose below any restrictions on them.
In order for the best linear prediction to be well-defined, we maintain the following assumption hereafter:
Let $b^0=(b^{0,1},...,b^{0,p})\in\mathbb R^p$ be such that $EL(Y|X)=X'b^0$. Usually, researchers are interested in specific components of $b^0$, rather than in the whole vector $b^0$. Therefore, in the following we seek to (partially) identify and estimate $b_d:=d'b^0$, for some $d\in \mathbb R^p$. For instance, if we focus on $b^{0,2}$, the second component of $b^0$, we let $d=(0,1,0,...,0)$.
Our setup includes at least three cases of broad interest.
\paragraph{Proxies for the covariate of interest.} In this case, we are interested in the relationship between a covariate of interest $X_o$ and $Y$. However, we do not observe $X_o$ in the $Y$ dataset, but only proxies $W_a$ of $X_o$. These proxies are also observed in the $X_o$ dataset. This type of situation arises very frequently in empirical microeconomics. A standard strategy in this case is to use two-sample two-stage least squares (TSTSLS). Namely, one first regresses $X_o$ on $W_a$ in the $X_o$ dataset. Then, we regress $Y$ on the predicted $X_o$ in the $Y$ dataset. Importantly though, this strategy implicitly relies on the following exclusion restriction:
This assumption is often restrictive. In our first application below, for instance, $Y$ corresponds to patent approval, $X_o$ is the vector of race dummies and $W_a$ is the vector of applicants' last names. Given that patent reviewers always observe last names but typically do not observe race directly, (ref) seems unlikely to hold. The method we develop in this paper will allow us to (partially) identify $EL(Y|X_o)$ without imposing (ref).
\paragraph{Missing controls.} In this case, we are interested in recovering the effect of $X_c$ on $Y$ using data from a first dataset. However, one or several key control variables ($X_o$) are missing from this dataset. Our setup applies to situations where the control variables are observed in a second dataset, together with $X_c$. In this sense, our framework complements a growing literature that investigates how credible unconfoundedness is, by allowing researchers to rely on unconfoundedness in a broad range of data combination environments AET05,oster2019unobservable, DMP22.
As above, researchers in this context may also have access to auxiliary variables, $W_a$, which are not included as covariates in the main regression, e.g., because they would be “bad controls”. As shown below, these variables may still carry informational content regarding the regression coefficients of interest.
\paragraph{Mediation analysis.} In this case, we are interested in the effect of $X_o$ on a given outcome $Y$. As above, we do not observe $X_o$ and $Y$ in the same dataset. A possible and frequent reason is that $Y$ is a long-run outcome, which is not observed in the data including $X_o$. On the other hand, both datasets may include other outcomes $W_a$, such as short-run outcomes.
To identify in this environment the causal effect of $X_o$ on $Y$ (which is $b^0$ under suitable randomization conditions on $X_o$), a common strategy is to rely on a surrogacy assumption prentice1989surrogate. In our setup, this corresponds to
In other words, one assumes that the effect of $X_o$ on $Y$ is entirely mediated by $W_a$.\footnote{One may also include additional covariates $X_c$ observed in both datasets, in which case (ref) becomes $EL(Y|X, W_a)=EL(Y|X_c, W_a)$.} However, Condition (ref) typically is a strong restriction. For instance, it is reasonable to assume that long-run earnings ($Y$) depend on human capital, even conditional on short-run earnings ($W_a$). Then, if job training ($X_o$) affects human capital, (ref) will fail to hold in general. Our results below imply that one can still obtain simple, sharp bounds on $b^0$, without relying on (ref).
Before presenting our identification results, we introduce additional notation. We denote by $\mathcal{B}_d$ the identified set of $b_d$ and let $\overline{b}_d$ and $\underline{b}_d$ be their sharp upper and lower bounds, namely $$\overline{b}_d =\sup\{d'b: \; b\in \mathcal{B}\}, \quad \underline{b}_d =\inf\{d'b: \; b\in \mathcal{B}\},$$ where $\mathcal{B}$ is the identified set of $b^0$. We focus in the following solely on $\overline{b}_d$, which is without loss of generality since $\underline{b}_d=-\overline{b}_{-d}$.
For any random variables $A$ and $B$, we let $F_A$ denote the cumulative distribution function (cdf) of $A$, $f_A$ its density, and $F_{A|B}$ the cdf of $A$ given $B$. We also let $F_A^{-1}(t):=\inf\{x:F_A(x)\ge t\}$ denote the quantile function of $A$; we denote similarly by $F_{A|B}^{-1}$ the quantile function of $A$ given $B$. We let $\text{Supp}(A)$ (resp. $\text{Supp}(A|B)$) denote the support of the probability distribution of $A$ (resp., of $A$ given $B$). For any vector $v$, we let $v_k$ denote its $k$-th element and $v_{-k}$ the vector obtained by removing $v_k$ from $v$. We also let $e_{k,r}$ denote the $k$-th canonical vector of $\mathbb R^r$. For any set $S$, we let $|S|$ denote its cardinality. Finally, we denote by $\mathcal{U}[0,1]$ the uniform distribution over $[0,1]$ and by $\mathcal{N}(\mu,\Sigma)$ the multivariate normal distribution with mean $\mu$ and covariance matrix $\Sigma$.
We first consider a case without nontrivial common variable ($W=X_c=1$), so that $X=(X'_o,1)'\in \mathbb R^{p}$. Our main result shows that $\mathcal{B}$ is convex and compact, and characterizes $\overline{b}_d$ for any $d\in\mathbb R^p\backslash\{0\}$. Below, we introduce the variable $\eta_d$ as follows. First, let $(d_2,...,d_p)$ be $(p-1)$ vectors in $\mathbb R^{p}$ such that $(d,d_2, ...,d_p)$ forms a basis of $\mathbb R^p$. Let $M$ denote the corresponding matrix and let $T=M^{-1}X$. Then, let $$\eta_d := T_1 - EL[T_1|T_{-1}].$$ In words, $\eta_d$ is the residual of the (population) regression of $T_1$ on $T_{-1}$. Note that $\eta_d$ does not depend on which exact vectors $(d_2,...,d_p)$ are chosen. Also, if $d=e_{k,p}$, $\eta_d$ is simply the residual of the regression of $X_k$ on $X_{-k}$. Finally, if $p=2$ and $d=(d_1,0)'$, $\eta_d=(X_o-E(X_o))/d_1$.
The first part of the theorem states that $\mathcal{B}$ is a convex, compact set included in the ellipsoid $\mathcal{E}$. Also, $(0,...,0,E[Y])'\in\mathcal{B}$: in the absence of common variables, we can always rationalize that $Y$ and $X$ are independent. Since the identified set $\mathcal{B}$ is non-empty, closed, and convex, $\overline{b}_d$ is equal to the so-called support function of $\mathcal{B}$. As a result, the knowledge of $\overline{b}_d$ for all $d\in\mathbb R^p\backslash\{0\}$ characterizes $\mathcal{B}$.
In the case of a single regressor (and the intercept) and $d=(1,0)'$, Equation (ref) reduces to
On the other hand, the true coefficient satisfies $b_d=b^{0,1}=E[(X_o - E(X_o))Y]/V(X_o)$. Thus, (ref) indicates that the sharp upper bound on the unknown term $E[X_oY]$ is $E[F_{X_o}^{-1}(U) F_{Y}^{-1}(U)]$. This is well-known, and corrresponds to the so-called Cambanis-Simons-Stout inequality cambanis1976inequalities. The logic is that (i) $F_{X_o}^{-1}(U)$ and $F_{Y}^{-1}(U)$ are distributed as $X_o$ and $Y$, since $U$ is uniformly distributed, and (ii) these two variables exhibit maximal positive dependence. The exact meaning of (ii) is that the copula of $F_{X_o}^{-1}(U)$ and $F_{Y}^{-1}(U)$ corresponds to the Fr\'echet-Hoeffding upper bound.
With multiple regressors, (ref) cannot be directly deduced from the Cambanis-Simons-Stout inequality. To get some intuition on (ref), suppose that $d=e_{1,p}$. Then, $\eta_d$ is the residual of the linear regression of $X_1$ on $X_{-1}$. If we observed $(Y,X)$, the coefficient of $X_{1}$ in the best linear prediction of $Y$ by $X$ would be $E[\eta_d Y]/E(\eta_d^2)$, by the Frisch-Waugh-Lovell theorem. Now, if we only know the marginal distributions of $\eta_d$ and $Y$, the numerator in (ref) is simply the upper bound of $E[\eta_d Y]$. That the sharp upper bound $\overline{b}_d$ satisfies (ref) is not obvious, however, because we also know the distribution of $X_{-1}$ conditional on $\eta_d$, in addition to the marginal distribution of $\eta_d$. This could, in principle, lead to $\overline{b}_d < E[F_{\eta_d}^{-1}(U) F_Y^{-1}(U)]/E(\eta_d^2)$. Theorem (ref) shows that this is not the case: the conditional distribution of $X_{-1}$ does not carry any additional information about $E[\eta_d Y]$. Although this can be deduced from Lemma 3.3 in delon2023generalized, we propose an alternative proof, which has the advantage of being constructive.
pacini2019two also obtains bounds on $b_d$, see his Theorem 1. However, it turns out that when $X$ is multivariate, his bound is only an outer bound rather than the sharp bound $\overline{b}_d$ on $b_d$. In Appendix (ref), we detail why this is the case, and provide an illustration showing that the sharp bounds given by Theorem (ref) above can in practice be substantially tighter than Pacini's bounds.
Let us now turn to the situation where some covariates are observed in both datasets. We define as above $\eta_d$, with the sole difference that now $X=(X'_o,X'_c)'$. Next, let $\delta_d$ and $\nu_d$ be such that $EL(\eta_d|W)=W'\delta_d$ and $\nu_d:=\eta_d-W'\delta_d$. Define $\delta_Y$ and $\nu_Y$ similarly, with $Y$ in place of $\eta_d$. The following theorem is the counterpart of Theorem (ref) with common variables.
Essentially, the first part of the theorem follows by first applying Theorem (ref) conditional on $W$ and then integrating over $W$. The second part exploits Theorem (ref) but conditioning on $g(W)$ instead of $W$. The sharp bound $\overline{b}_d$ has a simple expression, but it involves the conditional quantile functions $F_{ \nu_d|W}^{-1}$ and $F_{\nu_Y|W}^{-1}$. Thus, estimating this sharp bound involves estimating these two nonparametric functions, which could be cumbersome in practice. On the other hand, when $g(W)$ has a finite support, the outer bound $\overline{b}^g_d$ is elementary to estimate and does not suffer from any curse of dimensionality. Moreover, this bound is actually sharp when $\nu_Y\perp \!\!\! \perp W|g(W)$ and $\nu_d\perp \!\!\! \perp W|g(W)$, as is the case for instance with $g(W)=1$, if $Y$ and $\eta_d$ follow a linear location model in $W$.
\paragraph{How do common regressors affect identification?}
Even without auxiliary variables $W_a$, the identified interval on the coefficients of $X_o$ may exclude 0 in the presence of common regressors, implying that the sign of these coefficients is identified. To see this, suppose that dim$(X_o)=1$, $X_o=f_1(X_c) + \zeta_o$, $Y = g_1(X_c)+\zeta_Y$ and $\zeta_o| X_c\sim\mathcal{N}(0,\sigma^2_o)$, $\zeta_Y| X_c\sim\mathcal{N}(0,\sigma^2_Y)$. Let also $f(X_c):=f_1(X_c) - EL(f_1(X_c)|X_c)$ and $g(X_c):=g_1(X_c) - EL(g_1(X_c)|X_c)$. Using Equation (ref), the fact that by construction $\delta_d=0$, and the normality of $\zeta_o$ and $\zeta_Y$, we obtain that the bounds on $b^{0,1}$ satisfy
where $\eta_{e_1}=f(X_c) + \zeta_o$. In particular, if $|E\left[f(X_c)g(X_c)\right]| > \sigma_o \sigma_Y$, 0 is excluded from the identified set of $b^{0,1}$. This occurs when $X_o$ and $Y$ strongly depend on $X_c$ in a nonlinear way, so that $E\left[f(X_c)g(X_c)\right]$ dominates the contribution from independent terms (namely, $\sigma_o \sigma_Y$). In the extreme case where $X_o$ and $Y$ are deterministic functions of $X_c$, so that $\sigma_o =\sigma_Y=0$, we obtain point identification.
That said, the identified interval on the coefficients of $X_o$ may widen when including common covariates. Even if observing $X_c$ in both dataset does increase the information on the joint distribution of $(Y, X_o)$, the parameter we consider also changes. In particular, the denominator $E[\eta_d^2]$ in (ref) may substantially decrease, if the $R^2$ of the linear regression of $X_o$ on $X_c$ is large. To illustrate this, suppose that $(X_o,X_c')'\sim\mathcal{N}(0,\Sigma_o)$ and $(Y,X_c')\sim\mathcal{N}(0,\Sigma_Y)$, with
Then, some algebra shows that without observing $X_c$, $[\underline{b}_{e_1},\overline{b}_{e_1}]=[-1,1]$. With $X_c$, on the other hand, $$[\underline{b}_{e_1},\overline{b}_{e_1}]=\left[-\sqrt{\frac{1-\rho^2_Y}{1-\rho^2_o}},\; \sqrt{\frac{1-\rho^2_Y}{1-\rho^2_o}}\right].$$ Thus, the interval $[\underline{b}_{e_1},\overline{b}_{e_1}]$ shrinks if $|\rho_Y|>|\rho_o|$ but widens otherwise.
Finally, we may also identify the sign of components of $b_c$, the regression coefficient of $X_c$. In fact, $b_c$ may even be point identified: if $X_c$ and $X_o$ are uncorrelated, $b_c$ is simply the coefficient of the regression of $Y$ on $X_c$.
\paragraph{The role of auxiliary variables.}
By observing auxiliary variables $W_a$ that are not in the regression of interest, we increase the available information without modifying the parameter of interest. Then, the interval $[\underline{b}_{e_1},\overline{b}_{e_1}]$ always shrinks (at least weakly so). This may lead to excluding 0 from $\mathcal{B}$ even without common variables $X_c$, a case that occurs whenever
for some $d\in\mathbb R^p$. Intuitively, (ref) requires enough dependence between $Y$ and $W_a$ and between $X_o$ and $W_a$. For instance, if $(W,X_o)\sim\mathcal{N}(0,\Sigma_o)$ and $(W,Y)\sim\mathcal{N}(0,\Sigma_Y)$, with $\Sigma_o$ and $\Sigma_Y$ as in (ref), we obtain $$[\underline{b}_{e_1},\overline{b}_{e_1}]=\left[\rho_o \rho_Y - \sqrt{(1-\rho_o^2)(1-\rho_Y^2)},\; \rho_o \rho_Y + \sqrt{(1-\rho_o^2)(1-\rho_Y^2)}\right].$$ Then, $0\not\in[\underline{b}_{e_1},\overline{b}_{e_1}]$ if and only if $\rho^2_o \rho^2_Y > (1-\rho_o^2)(1-\rho_Y^2)$. This holds when $W_a$ is sufficiently correlated with $X_o$ and $Y$. For instance, when $\rho_o=\rho_Y$, this occurs if and only if $W_a$ explains more than half of the variance of $X_o$ ($\rho_o^2>1/2$).
A leading case with auxiliary variables is the case of surrogates. Recall that in this case, $X_o$ corresponds to the treatment variable, $Y$ is a long-run outcome while $W_a$ denotes short-run outcomes (surrogate variable). Then, Theorem (ref) yields two sets of bounds, sharp and outer, on the effect of $X_o$ on $Y$, without imposing a surrogacy assumption.
\paragraph{Link with TSTSLS.}
Recall that the TSTSLS estimand identifies $b^0$ if the coefficient of $W_a$ in the “long” regression of $Y$ on $(X,W_a)$ is 0. Now, the discussion above (“How do common regressors affect identification?”) shows that if one views $W_a$ as a common regressor, 0 may not belong to the identified set of the regression coefficient of $W_a$. This implies that the exclusion restriction underlying the TSTSLS estimand can actually be rejected by the data. As a simple example, suppose that $X_o$ and $W_a$ are not correlated. Then, the coefficient of $W_a$ in the “long” regression is equal to the coefficient of $W_a$ in the “short” regression of $Y$ on $W_a$, and this coefficient may not be 0. Beyond this particular case, the TSTSLS estimand for the coefficient $b^{0,k}$ may not belong to the sharp identified set $[\underline{b}_{e_k}, \overline{b}_{e_k}]$, something we illustrate in our second application below.
We have maintained thus far that the two samples at hand are drawn from the same population. While this is a standard assumption in the data combination literature, it is important to consider the extent to which this can be relaxed. To this end, let us introduce the binary variable $D$, with $D=1$ (resp. $D=0$) if we consider the $Y$ dataset (resp. the $X_o$ dataset). Then, our setup implies that we only observe the distributions of $(W,Y)|D=1$ and $(W,X_o)|D=0$, assuming that $D\perp \!\!\! \perp (W,X_o,Y)$. With common variables, this condition can be tested, since it implies $F_{W|D=1}=F_{W|D=0}$. If this implication is rejected, we can weaken the independence assumption by assuming instead that
In words, the first condition imposes that conditional on $W$, the two datasets are drawn from the same population, while the two populations corresponding to $D=0$ and $D=1$ may differ in their marginal distributions of $W$. The second condition in (ref) implies that the joint distribution of $(D,W)$, and thus the “propensity score” $p(W):=P(D=1|W)$, can be retrieved from the knowledge of the distributions of $W|D=0$ and $W|D=1$.
If (ref) holds, the sharp upper bound $\overline{b}_d$ can be obtained by reasoning as in Theorem (ref), using an inverse probability weighting scheme. Specifically, to identify $\delta_Y= E[WW']^{-1} E[WY]$ (and then $\nu_Y$), we cannot directly regress $Y$ on $W$ conditional on $D=1$. Yet, we can recover it by considering instead a weighted regression, as $$\delta_Y = E\left[\frac{DWW'}{p(W)}\right]^{-1}E\left[\frac{DWY}{p(W)}\right].$$ We can identify $\delta_d$ (and then $\nu_d$) similarly, using the weights $(1-D)/(1-p(W))$. Then, Equation (ref) is replaced by: $$\overline{b}_d=\frac{1}{E\left[\frac{(1-D)\eta_d^2}{1-p(W)}\right]} \left\{\delta_d'E(WW')\delta_Y + E\left[F_{\nu_d|W,D=0}^{-1}(U|W)F_{\nu_Y|W,D=1}^{-1}(U|W)\right]\right\}.$$
Another point to note is that if the two populations differ, the parameter of interest may correspond to one of the two populations only. For instance, one may consider, instead of $EL(Y|X)$, $EL(Y|X,D=1)$. In this case, $\delta_Y$ is given by $E[WW'|D=1]^{-1} E[WY|D=1]$ and is thus obtained by an unweighted regression, whereas $\delta_d$ (and then $\nu_d$) is obtained by regressing $\eta_d$ on $W$ with weights $p(W)/(1-p(W))$. The upper bound $\overline{b}_d$ becomes $$\overline{b}_d=\frac{E(D)\big\{ \delta_d'E(WW'|D=1)\delta_Y +E\left[F_{\nu_d|W,D=0}^{-1}(U|W)F_{\nu_Y|W,D=1}^{-1}(U|W)|D=1\right]\big\}}{E\left[(1-D)\eta_d^2p(W)/(1-p(W))\right]}. $$ Finally, another practically relevant situation is one in which one sample is drawn from a subpopulation of the population from which the other sample is drawn. Then, we identify instead (for instance) the distribution of $(Y,W)$ given $D=1$ and the distribution of $(X,W)$. In this case and if we focus as above on $EL(Y|X,D=1)$, we obtain a similar upper bound on $\overline{b}_d$ as above, with just a few differences. First, $\delta_d$ (and then $\nu_d$) is obtained by regressing $\eta_d$ on $W$ with weights $p(W)$. Second, we now have
Note that in this case and the one before, we do not require the joint independence condition in (ref) but only $X_o\perp \!\!\! \perp D|W$.
In practice, one may have access to auxiliary variables that appear in the dataset of $Y$ only, or in the dataset of $X_o$ only. For instance, suppose we identify the distributions of $(W,Y,Z)$ from one dataset and that of $(W,X_o)$ from the other. The following proposition shows that, for identification purposes, knowing the conditional distribution of $Z|W,Y$ provides no additional information. Hereafter, we let $\mathcal{B}_Z$ denote the identified set of $b^0$ when observing some auxiliary non-common variables $Z$.
A similar result clearly holds if we consider instead a variable that appears only in the dataset of $X_o$. The bottom line is that, among variables not included in the regression, only those that are common across the two datasets are relevant for identification.
Consider first the simplest situation where we only observe two independent samples, $\mathcal{S}_1:=(Y_i)_{i=1,...,n}$ and $\mathcal{S}_2:=(X_j)_{j=1,...m}$. Let $\widehat{\eta}_{dj}$ denote $j$'s residual in the sample regression of $T_1$ on $T_{-1}$ (recall the definition of $T$ at the beginning of Section (ref)). To ease notation, we let hereafter $F:=F_Y$ and $G:=F_{\eta_d}$, and let $F_n$ and $\widehat{G}_m$ denote the empirical cdfs of $(Y_i)_{i=1,...,n}$ and $(\widehat{\eta}_{dj})_{j=1,...,m}$, respectively. From Theorem (ref), we have $\overline{b}_d = \int_0^1 F^{-1}(t)G^{-1}(t)dt/E(\eta_d^2)$. Then, we consider the plug-in estimator of $\overline{b}_d$: $$\widehat{\overline{b}}_d = \frac{\int_0^1 F_n^{-1}(t)\widehat{G}^{-1}_m(t)dt}{\widehat{E}(\widehat{\eta}_d^2)},$$ where $\widehat{E}(\widehat{\eta}_d^2)$ denotes the empirical variance of $(\widehat{\eta}_{dj})_{j=1,...,m}$. Remark that when $m=n$, we simply have, denoting by $Y_{(i)}$ the $i$-th order statistic of $(Y_i)_{i=1,...,n}$ (similarly for $\widehat{\eta}_{d(i)})$: $$\int_0^1 F_n^{-1}(t)\widehat{G}^{-1}_m(t)dt = \frac{1}{n}\sum_{i=1}^n Y_{(i)}\widehat{\eta}_{d(i)}.$$ Otherwise, we can still compute the numerator of $\widehat{\overline{b}}_d$ at low cost. To see this, note that for any real-valued variables $U_1$, $U_2$ with finite second moments and cdfs $F_1, F_2$,
where $W_2$ is the Wasserstein-2 distance, $W_2(F_1,F_2):=(\int (F_2^{-1}(t)-F_1^{-1}(t))^2dt)^{1/2}$. For variables with support size of $n$ and $m$ respectively, as is the case here, we can then compute $W_2(F_1,F_2)$ with algorithms of complexity $O(m + n)$, see, e.g., rubner2000earth.
Let us now consider the case where common variables are observed. Specifically, we now assume to observe $\mathcal{S}_1:=\{(Y_1, W^{(1)}_1),....,(Y_n,W^{(1)}_n)\}$ and $\mathcal{S}_2:=\{(X_{o1}, W^{(2)}_1),....,(X_{om},W^{(2)}_m)\}$, where $W^{(1)}_i$ and $W^{(2)}_j$ are both distributed as $W$. We add the exponents $(\ell)$ to indicate that $W^{(\ell)}\in \mathcal{S}_\ell$. Recall from Theorem (ref) that $\overline{b}_d$ involves the nonparametric functions $F^{-1}_{\nu_d|W}$ and $F^{-1}_{\nu_Y|W}$. To avoid their estimation, we consider instead the outer bound $\overline{b}^g_d$ for a function $g$ taking finitely many values $(g_1,...,g_K)$. Then, $$\overline{b}^g_d = \frac{1}{E(\eta_d^2)}\left\{\delta_d' E(WW') \delta_Y + \sum_{k=1}^K p_k F_{|k}^{-1}(U) G_{|k}^{-1}(U)\right\},$$ where $p_k:=P(g(W)=g_k)$, $F_{|k}:=F_{\nu_Y|g(W)}(.|g_k)$ and $G_{|k}:=F_{\nu_d|g(W)}(.|g_k)$. Again, we consider a plug-in estimator of $\overline{b}^g_d$: $$\widehat{\overline{b}}{}^g_d = \frac{1}{\widehat{E}(\widehat{\eta}_d^2)}\left\{\widehat{\delta}_d' \widehat{E}(WW') \widehat{\delta}_Y + \sum_{k=1}^K \widehat{p}_k \int_0^1 \widehat{F}_{|k}^{-1}(u)\widehat{G}_{|k}^{-1}(u)du\right\},$$ where $\widehat{F}_{|k}$ (resp. $\widehat{G}_{|k}$) is the empirical cdf of $\widehat{\nu}_Y$ (resp. $\widehat{\nu}_d$) on the subsample of $\mathcal{S}_1$ satisfying $g(W^{(1)}_i)=g_k$ (resp., the subsample of $\mathcal{S}_2$ satisfying $g(W^{(2)}_j)=g_k$). The estimators $\widehat{E}(WW')$ and $\widehat{p}_k$ are simply obtained by combining the two samples, e.g.,
\paragraph{Choice of $g(.)$.} If $W$ is finitely supported, one can simply let $g(W)=W$. Yet, if $W$ takes many values, it is convenient to group some of these values together, so that none of the $(\widehat{p}_k)_{k=1,...,K}$ is too small and the asymptotic framework below remains a good approximation. When $W$ is not finitely supported, recall from Theorem (ref) that $\overline{b}^g_d$ is sharp if $\nu_Y$ $\perp \!\!\! \perp W|g(W)$ and $\nu_d\perp \!\!\! \perp W|g(W)$. Hence, we can expect tight bounds if $g(W)$ captures most of the dependence between $(\nu_Y,\nu_d)$ and $W$. Since $\nu_Y$ and $\nu_d$ are already residuals, we seek to capture possible heteroskedasticity by regressing $|\nu_Y|$ and $|\nu_d|$ linearly on $W$. This yields two indices, $W'\widehat{\varsigma}_Y$ and $W'\widehat{\varsigma}_d$. The underlying idea is that if $Y$ and $\eta_d$ satisfy a linear location-scale model, namely $Y= W'\delta_Y+(W'\varsigma_Y) \xi_Y$ with $\xi_Y\perp \!\!\! \perp W$ and similarly for $\eta_d$, then $\nu_Y\perp \!\!\! \perp W|g(W)$ and $\nu_d\perp \!\!\! \perp W|g(W)$ hold with $g(W) = (W'\varsigma_Y,W'\varsigma_d)$. However, this construction does not ensure that $g$ is finitely supported. To address this, we perform $K$-means clustering on $(W'\widehat{\varsigma}_Y, W'\widehat{\varsigma}_d)$. This yields a function $g$ taking $K$ values only. The choice of $K$ is discussed in Section (ref) below.
We now turn to the asymptotic properties of $\widehat{\overline{b}}_d$, and the construction of confidence intervals on $b_d$. For conciseness, we focus on the case without common variables; we briefly discuss the effect of these variables at the end of the section.
We first establish the asymptotic normality of $\widehat{\overline{b}}_d$, under the following assumptions.
We consider in Assumption (ref) three possibilities, depending on whether $\eta_d$ and $Y$ are finitely supported or not. The first case corresponds to $\eta_d$ being finitely supported. In such a case, $Y$ can be continuous or discrete, as long as, in the latter case, there is no $(h,y)$ such that $F(y)=G(h) \in (0,1)$. The second case corresponds to $Y$ being finitely supported and $\eta_d$ continuous. The third case corresponds to the two variables being, loosely speaking, continuous (actually, case (iii) is compatible with $Y$ having point masses, if we let $Z=\eta_d$). Then, we impose not only moment conditions but also (ref). This condition holds on $\text{Supp}(Z)\cap[0,\infty)$ for all distributions that have increasing hazard rates, such as log-concave distributions (as their survival function is then log-concave). It also holds for many distributions with decreasing hazard rates, such as Pareto and Weibull distributions. More generally, we expect Condition (ref) to be mild, since for any continuous probability measure $\mu$ with cdf $F$, density $f$ and supremum of support equal to $\overline{x}\le \infty$, we have, for all $A<\overline{x}$ satisfying $F(A)>0$, $$\int_A^{\overline{x}} \frac{f(x)}{F(x)(1-F(x))}dx \ge \int_A^{\overline{x}} (-\ln[1-F(x)])'dx =\infty.$$ On the other hand, for any $C_1,C_2>0$, $$\int_A^{\overline{x}} C_1\wedge \frac{C_2}{|x|\ln(1+|x|)^2}dx < \infty.$$ Thus, one cannot have $f(x)/[F(x)(1-F(x))]\le C_1\wedge C_2/(|x|\ln(1+|x|)^2)$ for all $x$ large enough; and similarly one cannot have $f(x)/[F(x)(1-F(x))]\le C_1\wedge C_2/(|x|\ln(1+|x|)^2)$ for all $x$ small enough.
To define the asymptotic distribution, we introduce additional objects. First, let $h(x):= \int_0^1 F^{-1}[G(x^-)+u(G(x)-G(x^-))]du$ and
These four variables correspond to the influence functions of respectively $\sqrt{m}(\widehat{E}(\widehat{\eta}^2_d)- E(\eta^2_d))$, $\sqrt{m} \int_0^1 F^{-1} (\widehat{G}_m^{-1} - G_m^{-1})dt$, $\sqrt{m}\int_0^1 F^{-1}(G^{-1}_m-G^{-1})dt$, and $\sqrt{n}\int_0^1 G^{-1}(F^{-1}_n-F^{-1})dt$, with $G_m$ the empirical cdf of the $(\eta_{dj})_{j=1,...,m}$ (note that $G_m$ cannot be computed in practice, since the $(\eta_{dj})_{j=1,...,m}$ are unobserved).
\paragraph{Remarks on the result.} First, we comment on the assumptions underlying Theorem (ref). We allow not only for $\lambda\in (0,1)$, but also for $\lambda=0$ or $\lambda=1$, which corresponds to cases where one sample is much larger than the other. In these cases, the asymptotic variance $V_d$ simplifies. Also, when $\min(|\text{Supp}(X)|,|\text{Supp}(Y)|)<\infty$, we obtain weak convergence under minimal conditions; note that $E[\|X\|^4]<\infty$ is close to being necessary for the OLS estimator $\widehat{\gamma}$ of the regression of $T_1$ on $T_{-1}$ to be $\sqrt{m}-$consistent.
When $\min(|\text{Supp}(X)|,|\text{Supp}(Y)|)=\infty$, the conditions we impose are probably not minimal, but note that a moment of order 4 for $Y$ and $\eta_d$ seems necessary in view of (ref) and the discussion of Theorem 1 in del2019central. Moreover, closely related results in the literature on the asymptotic normality of $W_2(F_n, G_m)$ impose strong restrictions.\footnote{By the proof of Point 2 of Theorem (ref), we obtain, under Assumptions (ref)-(ref), the asymptotic normality of $(nm/(n+m))^{1/2} (W_2(F_n, G_m)-W_2(F,G))$.} In particular, instead of Assumption (ref)-(iii), Proposition 2.3 in del2019central imposes strong and high-level conditions (see (2-7)-(2.9) in their paper), while Theorem 14 in berthet2020 also imposes strong regularity conditions. In particular, because their Assumption (FG) must hold for both the left and right tails of the distributions, one can show that their subconditions (FG1) and (FG3) already imply (up to letting $\varepsilon=0$) Assumption (ref)-(iii) for both $Z=Y$ and $Z=\eta_d$.\footnote{On the other hand, both berthet2020 and del2019central also consider more general Wasserstein distances than just $W_2$.}
\paragraph{Sketch of the proof.} In a first step, we account for the fact that $\eta_d$ and $E[\eta_d^2]$ are estimated. This requires in particular to show that $$\sqrt{m} \int_0^1 F_n^{-1}(\widehat{G}^{-1}_m - G^{-1}_m)dt = - E[h(\eta_d)T'_{-1}]\sqrt{m}(\widehat{\gamma}-\gamma_0) + o_P\left(1\right),$$ where $\gamma_0$ is the limit in probability of $\widehat{\gamma}$. This result is not obvious; our proof relies in particular, again, on the Cambanis-Simons-Stout inequality. The second step is to study the asymptotic behavior of $(nm/(n+m))^{1/2} \int_0^1 [F_n^{-1}(t)G_m^{-1}(t)-F^{-1}(t) G^{-1}(t)]dt$. Here, we use the decomposition
where $r_{n,m} := \int_0^1 (F_n^{-1}(t)-F^{-1}(t))(G_m^{-1}(t)-G^{-1}(t))dt$. We prove that the first two terms $T_{1m}$ and $T_{2n}$ are asymptotically linear by adapting results on L-statistics, see in particular Theorem 1 in Chapter 19 of SW86. That the remainder term $r_{n,m}$ is negligible if $Y$ (say) is finitely supported follows from the continuity of $G^{-1}$ at the support points of $Y$. Note that if this continuity condition does not hold, we lose asymptotic normality; see del2024central for the exact distribution in such cases. If Assumption (ref)-(iii) holds, we relate instead the remainder term to bounds on the convergence rate of $W_2(F_n,F)$ and $W_2(G_m,G)$. However, existing results on such rates, and in particular Theorem 1 in fournier2015rate, are not sufficient for our purpose. Here, we improve upon their bound, which holds under weak restrictions, by leveraging in particular Condition (ref). We do this by linking $W_2(F_n,F)$ to the variance of order statistics, and relying on a lemma similar to Corollary 2.12 in boucheron2015; see Lemma (ref) in Online Appendix (ref).
We construct confidence intervals on $b_d$ using the asymptotic normality of $\widehat{\overline{b}}_d$ and a plug-in estimator of $V_d$. Specifically, let $\widehat{h}(x)= \int_0^1 F_n^{-1}[\widehat{G}_m(x^-)+u(\widehat{G}_m(x)-\widehat{G}_m(x^-))]du$ and
Then, define $$\widehat{V}_d := \frac{1}{\left(\frac{1}{m}\sum_{j=1}^m \widehat{\eta}_{dj}^2\right)^2}\times \left[\frac{n}{m(n+m)}\sum_{j=1}^m \left(\widehat{\psi}_{1j} + \widehat{\psi}_{2j} + \widehat{\psi}_{3j}\right)^2 + \frac{m}{n(n+m)} \sum_{i=1}^n \widehat{\psi}_{4i}^2 \right].$$ Note that $\widehat{V}_d$ depends on $d$; in particular, $\widehat{V}_{-d}$ is the estimator of the asymptotic variance of $\overline{b}_{-d}=-\underline{b}_d$. We then consider the following confidence intervals on $b_d$ with nominal level $1-\alpha$: $$\text{CI}_{1-\alpha} := \left[-\widehat{\overline{b}}_{-d} - z_{1-\alpha} \sqrt{\frac{n+m}{nm}\widehat{V}_{-d}}, \; \widehat{\overline{b}}_{d} + z_{1-\alpha} \sqrt{\frac{n+m}{nm} \widehat{V}_d}\,\right],$$ where $z_{1-\alpha}$ is the quantile of order $1-\alpha$ of a standard normal distribution. We can replace the usual quantile $z_{1-\alpha/2}$ by $z_{1-\alpha}$ since by Theorem (ref), the identified interval of $b_d$ is not reduced to a singleton ($\overline{b}_d >0>\underline{b}_{d}$) as long as $V(Y)>0$.
Once again, the proof of Theorem (ref) is not straightforward. In particular, two difficulties are (i) to prove convergence of $(1/m)\sum_{j=1}^m \widehat{h}(\widehat{\eta}_{dj}) T_{-1j}'$; (ii) to handle the terms including $\widehat{\psi}_{3i}$ and $\widehat{\psi}_{4i}$. For (ii), we rely in particular on an extension of Lemma A.1 in del2019central, see Lemma (ref) in Online Appendix (ref).
Given our focus on a finitely supported $g(W)$, the analysis is very similar to the case without common variables, so we mostly highlight the differences here, without providing a formal result for the sake of conciseness. First, the asymptotic variance of $\widehat{\overline{b}}_d$ includes additional terms due in particular to (i) the estimation of $\delta_d'E[WW']\delta_Y$; (ii) the estimation of the residual $\nu_Y$. The exact expression of the asymptotic variance, which includes eleven terms instead of four as above, is given in Online Appendix (ref).
Then, the construction of the confidence interval is similar to that described above, with one important difference, which is to allow for the possibility of point identification. To maintain size control, we rely on stoye2020simple to construct the confidence intervals. This method has the appealing features of not requiring any tuning parameter, being simple to compute, and relying on mild conditions, beyond the joint asymptotic normality of the lower and upper bounds. We implement this inference method in our Monte Carlo simulations (Section (ref)) and in the applications (Section (ref)).
We now study the finite sample performances of our estimators and inference method. We consider a single DGP encompassing three cases of available data: one in which only $Y$ and $X_o$ are available, one in which $X_c$ is also observed jointly and enters the main regression and one in which $W_a$, in addition to $X_c$, is observed. In the latter case, the parameters remain the same as in the second case. The DGP is as follows. We let $W_a \sim \mathcal{U}[0,1]$, $X_c \sim \mathcal{N}(0, 1)$ and
We fix $a_1=1$, $a_2=10$, $d_1=1$, $\sigma_\eta = 1$, $b_1=1$, $b_2=1$, $d_2=0.25$ and $\sigma_\varepsilon = 4$. The true bounds in the first two cases are obtained by simulations, whereas there is a closed-form expression in the last case. We fix $n=m$ and vary it from 400 to 4,800. We construct $g(W)$ as described in Subsection (ref), with $K=\max(2,\lfloor \min(n,m)^{0.2}\rfloor)$, where $\lfloor x\rfloor$ denotes the integer part of $x$; we discuss alternative choices of $K$ below. The results are displayed in Table (ref). We report the average of the estimated bounds (“Bounds”) and the average of the estimated 95% confidence intervals $\text{CI}_{1-\alpha}$ (“95% CI”) for $b^{0,1}$. We also report the mean difference between the length of the confidence sets and that of the identified set, see column “Ex. length” in the table. Finally, the column “Covg” corresponds to the minimum, over $b_1$ in the identified set of $b^{0,1}$, of the estimated probability that $b_1$ belongs to the confidence interval.
A couple of remarks are in order. First, as expected, the 95% confidence intervals shrink with the sample sizes $n$, approximately at the $n^{-1/2}$ rate in the three cases we consider. This is reflected in the evolution of the excess length across sample sizes. Second, the confidence intervals exhibit satisfactory coverage. In particular, coverages for all panels are generally close to the nominal 95% level, even for small sample sizes. Coverage rates are generally conservative, but still very close to the nominal level for the specification reported in Panel 3. This is remarkable: one would in principle need to use the continuous variable $g(W)=1+W_ad_1$ to obtain the sharp bounds, by Theorem (ref), whereas we instead rely on a finitely supported variable $g(W)$ with few points of support (from 3 to 5 when $n$ varies from $400$ to $4,800$). Third, and importantly, the identified set is much tighter in Panel 3 than in Panels 1 and 2. This illustrates the substantial identifying power of the auxiliary variable $W_a$. For this particular DGP, the identifying power - measured by the reduction in the length of the identified set - of $W_a$ is in fact larger than that of the common regressor $X_c$.
Table (ref) reports the computational time needed to compute the estimated bounds and associated confidence intervals. When $W_a$ is observed, this time also includes the $K$-means clustering we perform to compute $g(W)$. The main takeaway is that our procedure is very fast: it takes less than 1 second when observing $(X_c,W_a)$ with $n=m=12,000$, and less than 12 seconds with $n=m$ as large as 120,000.
Finally, we explore the effect of the tuning parameter $K$ on coverage; see Table (ref) below, where we consider two sample sizes ($n=1,200$ and $n=6,000$). As expected, increasing $K$ decreases the length of the CIs, but also reduces coverage. This probably reflects the fact that the estimated bounds become biased for larger $K$. On the other hand, coverage remains above 95% for $K =\max(2,\lfloor\min(n,m)^c\rfloor)$, $c<1/3$, suggesting that our baseline choice of $K$ with $c=0.2$ works well in practice.
We now illustrate our approach with two applications. We first study the influence of race on the probability of patent approval in the United States, revisiting recent work on this question Dossi23. We then investigate the relationship between students' risk and time preferences and educational achievement across countries hanushek2020culture.
In our first application, we investigate the existence and magnitude of racial and ethnicity gaps in science and innovation. This question has attracted much interest in the recent empirical literature Kerr08,ADQW24,Dossi23. A key challenge is that datasets typically do not measure race and ethnicity together with the outcome of interest. Using our notation, race/ethnicity is an outside regressor ($X_o$), with successful patent application being the outcome of interest ($Y$). Instead of $X_o$, we may observe other characteristics, such as the applicant's name. Then, in other datasets, we may observe these characteristics together with race and ethnicity. A commonly used strategy in this context is to impute race and ethnicity using applicant characteristics observed in both datasets. We take a different route and derive bounds that use both datasets without relying on the exclusion restriction implicit in the imputation approach.
Following Dossi23, we rely on two datasets. The first is the publicly available dataset released by the United States Patent and Trademark Office (USPTO) covering the universe of patent applications submitted in the United States. We use the Patent Examination research dataset (PatEx), which contains detailed information on all patent applications, including the full names of the applicants graham2015uspto. We restrict the sample to applications filed between January 2001 and December 2018 and focus on utility patents.\footnote{Utility patents, also referred to as “patents for invention”, constitute 90% of the patent documents issued by the USPTO in recent years.} We further restrict the sample to applicants based in the United States, and as in Dossi23, consider only the first inventor listed on the application.
We combine the PatEx dataset with data from the US Census. Namely, we use the information on the aggregate frequency of last names by race and ethnicity from the 2010 Decennial Census Surname Table comenetz2016frequently. 6.3 million different last names were recorded for 295 million people.\footnote{In the following, we neglect the statistical uncertainty related to this sample.} Among them, we use the publicly released frequency by race and ethnicity of the 162,254 last names that occur more than 100 times, representing 90.1% of the overall population. We consider in our analysis five different categories of race and ethnicity, namely: (i) Black or African American (11.99%), (ii) Asian and native Hawaiian and other Pacific Islander (4.86%), (iii) Hispanic or Latino (16.29%), (iv) American Indian or Alaska Native (1.76%), and (v) others, which includes White (64.40%) and those declaring to belong to two or more races (0.69%), and is used as our reference category.\footnote{Estimation results are robust to splitting the “two or more races” category evenly across the other racial categories.}
Estimation results are reported in Table (ref), where the first column presents the TSTSLS point estimates, the second the point estimates of our bounds, and the last column reports the 95% confidence intervals computed from our asymptotic normality results. We use applicant's last name as an auxiliary variable ($W_a$). Our bounds correspond to $\overline{b}_d^g$ where, to reduce the size of the vector $W$, we define $g(W)$ as $W_a$ unless the names $W_a$ appear $L=5$ times or less in the dataset of inventors, in which case we set $g(W)=0$ (Table (ref) in Appendix (ref) show that our results are robust to choosing $L=3$ or $L=10$ instead). Since inventors are a subset of the whole population, our bounds are plug-in estimates of Equation (ref) above. In contrast to our bounds, the TSTSLS estimates rely on an exclusion restriction. Namely, the applicant's last names is assumed not to be predictive of patent approval once conditioning on applicant's race and ethnicity. While this type of name-based exclusion restriction has frequently been used in applied work, its validity is far from obvious in this particular context. In fact, it does not seem unreasonable to think that, in contrast to the TSTSLS exclusion restriction, any racial discrimination in the patent approval process would operate largely through the applicant's last name. This is consistent with the information available to patent examiners, who always observe applicants' last names but do not observe race, and only infrequently interview them in person (see CKS02, 2002, and Avivi24, 2024 for discussions of the USPTO selection process).\footnote{One may argue that the TSTSLS identifies instead the effect of, e.g., having a Black- or Asian-sounding name. It is unclear whether this interpretation is warranted either. First, names could also predict other relevant characteristics. Second, this interpretation would require another exclusion restriction, namely that the coefficients of race in the “long” regression are zero, which is arguably strong as well.}
Turning to the results, a key takeaway is that the TSTSLS results that are obtained using last names as an exclusion restriction are fragile. Notably, while the TSTSLS estimates point to Black inventors being significantly less likely to be granted a patent, the bounds obtained with our method for this coefficient are wide, with a lower bound as large as $-0.729$ and an upper bound that is positive and large as well ($0.360$). While a similar conclusion holds for Hispanics and American natives, our bounds are somewhat more informative for the coefficient on Asians, with a lower bound of $-0.137$ and an upper bound of $0.187$. At any rate, these results indicate that the conclusions one would reach from the TSTSLS estimates of significant racial differences in the probability of being granted a patent crucially hinge on the underlying exclusion restriction.
A final point is that, although the bounds reported in Table (ref) tend to be wide, using last names as $W_a$ does yield substantial improvements over the simple bounds based solely on $Y$ and $X_o$. In particular, without $W_a$, the sharp lower bound for each of the four coefficients equals -1 and is therefore not informative.\footnote{This occurs here because (i) $P(Y=0)$ is larger than $P($race$)$ for all other races than White and multiracial applicants, and ii) $P(Y=1)$ is larger than $P($White$)$. The corresponding upper sharp bound is equal to $0.432$ for each of the four coefficients. Again, this is due to the particular configurations of $P(Y=1)$ and $P($race$)$. In our setup, 0.432 simply corresponds to $P(Y=0)/P($White$)$.} Hence, even without exclusion restrictions, observing last names in both datasets delivers meaningful informational gains.
Next, we investigate why our bounds are more informative for some races/ethnicities than others. To do so, we report in Figure (ref) the distribution of the racial frequencies conditional on last name, focusing on the names that are the most predictive of race and that together account for 10% of each sub-populations (here, Blacks and Asians - groups for which the bounds are relatively wide and more informative, respectively).This figure shows that last names are highly predictive of being Asian, much more so than for Blacks. This illustrates the connection between the informativeness of our bounds and the extent to which $W_a$ (inventor's last name) is predictive of $X_o$ (race/ethnicity).
We conclude this analysis by exploring further the effect on our bounds of using more auxiliary information, as measured by the auxiliary variables $W_a$. Since the Census Surname table only contains racial characteristics associated with last names at the aggregate level, we cannot use this data for this purpose. Instead, we leverage the fact that voter registration data in North Carolina (Historical Voter Registration Snapshots) records historical individual data about active and inactive voters registered in North Carolina, with information about their full names, city, race and ethnicity. Thus, restricting to the set of inventors residing in North Carolina (62,112 applications associated with 23,689 unique inventors), we are able to merge application data with individual data from the Historical Voter Registration Snapshots.\footnote{To be representative of the population in North Carolina over the period 2001-2018, we actually use snapshots of 2006, 2013, and 2019, keeping only information about full names, city, and recorded race and ethnicity of all the uniquely identified active and inactive voters over this period. This data is openly available at \href{https://www.ncsbe.gov/results-data/voter-history-data.}{www.ncsbe.gov/results-data/voter-history-data}.}
Table (ref) provides a comparison of different point estimates for our bounds, using different sets of $W_a$. A couple of comments are in order. The first one relates to the sample restrictions imposed by the common support requirements when using more comprehensive sets of $W_a$. The underlying reduction in sample size increases from around 3% when using last names only, to as much as 37% when using the complete name and city. Second, we only obtain very informative bounds in the latter case, in which we uniquely identify as much as 98% of the inventors. In other words, (very) high predictive power is needed to obtain tight bounds on the coefficients of interest.
Preference parameters, especially patience and risk taking, play an important role in human capital investment decisions. However, to our knowledge, no single data set jointly measures these preferences and test scores across countries. In the following, we build on the cross-country analysis of hanushek2020culture and combine data from the OECD's Programme for International Student Assessment (PISA) with the Global Preference Survey (GPS) to examine how students' time and risk preferences are associated with educational achievement.
PISA assesses achievement in mathematics, science and reading for random samples of 15-year-old students on a three-year cycle, providing repeated cross-sectional data representative of each country-by-wave cell. In the following, we consider as our main “Y dataset” the standardized math test scores over the seven waves of PISA testing, covering the period 2000-2018. Over this period, a total of 86 countries participated at least once. We combine these test scores with data from the Global Preference Survey falk2018global. The GPS provides scientifically validated data on several preference parameters from representative samples, of around 1,000 respondents in each country surveyed in 2012, measuring patience, risk taking, positive and negative reciprocity, altruism, and trust $(X_{o})$, for 49 different countries. The GPS also records gender for each respondent, which we use as a common regressor $X_c$, together with the country. Restricting the analysis to this subset of 49 countries yields test score data for a total of 1,992,276 students.\footnote{See Appendix (ref) for more details on the GPS dataset, especially on the measurements of preference parameters. As the PISA and GPS datasets are representative of the same common population after reweighting the observations by the corresponding survey weights, we use the survey weights in our analysis.}
In Table (ref), we compare our bounds on the coefficients of patience and risk taking with the TSTSLS estimates considered by hanushek2020culture, where both variables are imputed using country dummies. We consider alternative specifications depending on whether only patience, only risk taking or both are included in the regression. In Panel D, we also include other preference variables (positive and negative reciprocity, altruism and trust) as controls in the regression. Hence, in this last specification, $X_o$ is of dimension 6.
A couple of comments are in order. A first takeaway is that, in contrast to the previous application and despite the absence of any $W_a$ here, our bounds tend to be informative. This holds for both the coefficients of patience and risk-taking, and across all four specifications reported in the table. That the bounds remain informative is particularly noteworthy in Panel D, where we control for four additional preference parameters all included in $X_o$. One might indeed have expected that increasing the dimension of $X_o$ would cause the bounds to widen substantially, yet this is not the case here.
Related to this, for Panels C and D, the TSTSLS estimates of the coefficients associated with patience and risk-taking both lie outside of the estimated sharp bounds, and, for the risk-taking coefficient, outside of the 95% confidence intervals as well. One-sided tests of equality between the lower bound of our identified set on the coefficient of risk-taking and the TSTSLS point estimate in Panels C and D leads us to reject this hypothesis at the 10% level, consistent with a violation of the underlying exclusion restrictions. Recall that the TSTSLS estimator relies on the arguably strong assumption that countries do not affect test scores beyond their effects through risk aversion and patience. In terms of magnitudes, focusing on Panel D where we include additional preference controls, our bounds indicate that a 1 standard deviation (SD) increase in patience is at most associated with 0.977 SD increase in math test scores, against 1.122 SD using TSTSLS. Similarly, it follows from our bounds that a 1 SD increase in risk-taking is, at most, associated with a decline of 1.094 SD in student achievement, against a larger decline of -1.345 using TSTSLS.
Finally and importantly for practice, our inference method can be implemented at a low computational cost. Even though the datasets used in this analysis contain a very large number of observations (1,992,276 and 49,689), it takes 5 minutes only to reproduce the results of Panels A and B, 9 minutes for Panel C, and 21.5 minutes for Panel D (where $X_o$ is of dimension 6), using our R package RegCombinBLP.\footnote{We parallelize the computation over 15 CPUs on an Intel Xeon Gold 6130 CPU 2.10GHz with 382Gb of RAM.}
We study regression coefficients in a context where the outcome of interest and some of the covariates are observed in two different datasets that cannot be matched. This type of data combination environment arises very frequently in various empirical setups. The usual approach, which consists in imputing the outcome $Y$ or the outside regressors $X_o$ using auxiliary variables $W_a$, hinges on exclusion restrictions that may not hold in practice. We take a different route and derive sharp bounds on the regression coefficients using only the observed distributions. As they take a simple form, these bounds can be estimated at a low computational cost; we also derive simple and easy-to-compute confidence intervals.
We illustrate our method with two applications. The first studies racial disparities in patent approval, the second the effects of patience and risk-taking on test scores. The first application highlights that in some cases, results based on an imputation approach crucially rely on the underlying exclusion restriction; without it, uncertainty on the true coefficients of interet remains large. The second application shows that our bounds can be informative on the magnitude of the effects, and can also lead to reject the imputation-based approach.