Marginal treatment effects in the absence of instrumental variables

The marginal treatment effect (MTE), which was developed in a series of seminal papers by Heckman and Vytlacil (1999, 2001, 2005), has become one of the most popular tools in social sciences for describing, interpreting, and analyzing the heterogeneity of the effect of a nonrandom treatment. The MTE is defined as the expected treatment effect conditional on the observed covariates and normalized error term consisting of unobservable determinants of the treatment status. Although the definition of MTE doesn’t necessitate an instrumental variable (IV), nearly all the existing identification strategies for the MTE rely heavily on valid instrumental variation which can induce otherwise similar individuals into different treatment choices. This is mainly because the MTE has been regarded as an extension of and supplement to the standard IV method in the causal inference literature. However, in practical applications, the validity of IVs is usually difficult to justify. In the most common case of just-identification in which only one IV is available, the exogeneity or randomness of the IV is fundamentally untestable. The exclusion restriction, which is another condition underlying the validity of IVs, has also been challenged more often than not (e.g., van den Berg, 2007; Jones, 2015).

Motivated by the potential invalidity of the available instruments, and inspired by the observation that instruments are unnecessary in the definition of MTE, we model, identify, and estimate the MTE without imposing the standard IV assumptions of conditional independence, exclusion, and separability. Namely, we allow all the observed covariates to be statistically correlated to the error term and have a direct impact on the outcome in addition to an indirect impact through the treatment variable. The cornerstone of our model is a normalized treatment equation, which determines treatment participation by the propensity score crossing a reduced-form error term that is statistically independent of (although functionally dependent on) all the covariates. This is comparable to the conventional IV model, in which normalization is performed with respect to only the non-IV covariates. The independence of the reduced-form treatment error from all the covariates can partly justify a conditional mean independence assumption on the potential outcome residuals, which can ensure the additive separability of the MTE into observables and unobservables. This separability then facilitates the semiparametric identification of the MTE, given a linear restriction on the potential outcome regression functions and a nonlinear restriction on the propensity score function. We prove the identification by representing the MTE as a known function of conditional moments of observed variables. Intuitively, the identification power comes from the excluded nonlinear variation in the propensity score, which plays the role of the IV in exogenously perturbing the treatment status. Our identification strategy closely resembles the semiparametric counterpart of identification based on functional form (Dong, 2010; Escanciano et al., 2016; Lewbel, 2019; Pan et al., 2022). We build our result on the notion of identification based on functional form mainly because (i) it can lead to point identification and estimation, (ii) the required assumptions are regular, and (iii) the resulting estimation procedures are the most compatible with the standard procedures in conventional IV-MTE models. After identification, the MTE can be semiparametrically estimated by implementing the kernel-weighted pairwise difference regression for the sample selection model (Ahn and Powell, 1993) for each treatment status.

We contribute to the literature on program evaluation under essential heterogeneity by proposing a method for defining, identifying, and estimating the MTE in the absence of IVs. The main value of this IV-free framework of the MTE is threefold. First, it provides an approach for consistently estimating heterogeneous causal effects when a valid IV is hard to find. Second, it suggests an easy-to-implement means for testing the validity of a candidate IV, because it nests the models that assume exclusion restrictions. For instance, in studies on returns to education, the validity of most instruments for educational attainment is suspect, such as parental education, distance to the school, and local labor market conditions (Kédagni and Mourifié, 2020; Mourifié et al., 2020). Our framework enables a reliable evaluation of returns to education without imposing IV assumptions on the candidate instruments and implies a simple test for exclusion restrictions by $t$-testing the instruments’ coefficients in potential outcome equations. Third, even though the validity of the IV is verified, identification based on functional form can be invoked as a way to increase the efficiency of estimation or to check the robustness of the results to alternative identifying assumptions.

Early explorations of the identification of endogenous regression models when IVs are unavailable focused on linear systems of simultaneous equations, in which the lack of IVs is associated with hardly justifiable exclusion restrictions and insufficient moment conditions. The typical strategy for addressing this underidentification problem is to impose second- or higher-order moment restrictions to construct instruments by using the available exogenous covariates (e.g., Fisher, 1966; Lewbel, 1997; Klein and Vella, 2010). In particular, imposing restrictions on the correlation between the covariates and the variance matrix of the vector of model errors leads to identification based on heteroscedasticity (Lewbel, 2012). D’Haultfœuille et al. (2021) considered a more general nonseparable, nonparametric triangular model without the exclusion restriction, and established identification for the control function approach based on a local irrelevance condition. However, these results are restricted to the case of continuous endogenous variable or treatment. Lewbel (2018) showed that, in linear regression models, the moment conditions of identification based on heteroscedasticity can be satisfied when the endogenous treatment is binary, at the expense of strong restrictions on the model errors. Alternatively, Conley et al. (2012) and Van Kippersluis and Rietveld (2018) presented a sensitivity analysis approach for performing inference for linear IV models while relaxing the exclusion restriction to a certain extent. However, the linear models implicitly impose an undesirable homogeneity assumption on treatment effects.

The literature on sample selection models without the exclusion restriction is also closely related. A general solution to the problem of lack of excluded variables is partial identification and set estimation (Honoré and Hu, 2020, 2022; Westphal et al., 2022). The limitation of this approach is that the identified or estimated set may be too wide to be informative. Another solution is the approach of identification at infinity that uses the data only for large values of a special regressor (Chamberlain, 1986; Lewbel, 2007) or of the outcome (D’Haultfœuille and Maurel, 2013). Although identification at infinity can lead to point identification, it is typically featured as irregular identification (Khan and Tamer, 2010), and the derived estimate will converge at a rate slower than $n^{-1/2}$, where $n$ is the sample size. Heckman (1979) exploited nonlinearity in the selectivity correction function to achieve point identification and root-$n$ consistent estimation for the linear coefficients of a parametric sample selection model, which is the original version of identification based on functional form. However, Heckman’s model imposes a restrictive bivariate normality assumption on the error terms and thus poses the risk of model misspecification. Escanciano et al. (2016) extended Heckman’s approach to a general semiparametric model and established identification of the linear coefficients by exploiting nonlinearity elsewhere in the model. At the expense of the generality of their model, Escanciano et al. (2016)’s identification result was only up to scale and required a scale normalization assumption. This normalization would be innocuous if we know the sign of the normalized coefficient and are interested in only the sign, but not the magnitude, of the other coefficients. However, the magnitude of the linear coefficients is indispensable to the evaluation of a program or a policy. In addition, the identifying assumption about nonlinearity developed by Escanciano et al. (2016) implicitly rules out the case of the existence of two continuous covariates. We adapt the result of Escanciano et al. (2016) to the MTE model by amending the two defects. First, we take advantage of the specific model structure to relax the scale normalization and identify the magnitude of the linear coefficients. Moreover, the identifying assumption about nonlinearity will be simplified owing to the specific model structure. We give an intuitive interpretation of the new nonlinearity assumption. Second, we combine the nonlinearity assumption with a local irrelevance condition to allow for an arbitrary number of continuous covariates.

The rest of this paper is organized as follows. In Section 2, we introduce our model and the definition of MTE without IVs. In Section 3, we propose a possible set of sufficient conditions for the identification of the MTE, in place of the standard IV assumptions. The key conditions include semiparametric functional form restrictions and the conditional mean independence assumption, of which the latter implies the additive separability of the MTE into observed and unobserved components. Section 4 suggests consistent estimation procedures, and Section 5 provides an empirical application to Head Start. Section 6 concludes. Online appendices comprise (A) a proof of the primary identification result, (B) a discussion on variants of the nonlinearity assumption, and (C) an identification result for limited valued outcomes which entails a slight modification of the assumptions.

In the following, we denote random variables or random vectors by capital letters, such as $U$, and their possible realizations by the corresponding lowercase letters, such as $u$. We denote $F_{U}\left(\cdot\right)$ as the cumulative distribution function (CDF) of $U$ and $F_{U\left|X\right.}\left(\cdot\left|x\right.\right)$ as the conditional CDF of $U$ given $X=x$. Our analysis builds on the potential outcomes framework developed by Roy (1951) in econometrics and Rubin (1973a, b) in statistics. Specifically, we consider a binary treatment, denoted by $D$, and let $Y_{1}$ and $Y_{0}$ denote the potential outcomes if the individuals are treated ($D=1$) or not treated ($D=0$), respectively. The observed outcome is

and the quantity of interest is the counterfactual treatment effect $Y_{1}-Y_{0}$. We suppose that the treatment status is determined by the following threshold crossing rule:

where $X$ is a vector of pretreatment covariates, $1\left\{A\right\}$ is the indicator function of event $A$, $\mu\left(\cdot\right)$ is an unknown function, and $U$ is a structural error term containing unobserved characteristics that may affect participation in the treatment, such as the opportunity costs or intangible benefits of the treatment.

Compared with the conventional MTE model, we relax the independence and separability assumptions in the treatment equation (1) to account for the absence of IVs. First, we don’t assume that $X$ is stochastically independent of $U$; that is, no exogenous covariate is needed. Second, we allow the treatment decision rule to be intrinsically nonseparable in observed and unobserved characteristics, which is equivalent to relaxing the monotonicity assumption in the IV model (Vytlacil, 2002, 2006). Specifically, let $U_{x}$ be the counterfactual error term denoting what $U$ would have been if $X$ had been externally set to $x$, then the nonseparability of (1) means that at least two vectors $x$ and $\tilde{x}$ exist in the support of $X$ such that $U_{x}\neq U_{\tilde{x}}$ with positive probability. In particular, $U$ is allowed to depend functionally on $X$, as in the following example.

Example 1 illustrates that no generality is lost by modeling the treatment variable as equation (1) without imposing the independence and separability assumptions. We define the propensity score function as the conditional probability of receiving the treatment given the covariates,

and define the propensity score variable as $P\equiv\pi\left(X\right)$. Under the regularity condition that $F_{U\left|X\right.}\left(\cdot\left|x\right.\right)$ is absolutely continuous with respect to the Lebesgue measure for all $x$, the structural treatment equation (1) can be innocuously normalized into a reduced form:

where

is a normalized error term, which by definition follows standard uniform distribution conditional on $X$ and thus is stochastically independent of $X$ and $P$. This independence, which seems counterintuitive due to the functional dependence of $V$ on $X$, will be lost if we consider $V_{x}=F_{U\left|X\right.}\left(U_{x}\left|x\right.\right)$, the counterfactual variable of $V$ when $X$ is set to $x$. In general, the conditional distribution of $V_{x}$ given $X=\tilde{x}$ for $\tilde{x}\neq x$ depends on $\tilde{x}$, and the unconditional distribution of $V_{x}$ is not uniform.

Given this independence, we may consider the reduced-form treatment error $V$ as the orthogonalized unobservables with respect to the observables, or the unobservables that are projected onto the subspace orthogonal to the one spanned by the observables. Another interpretation of $V$ is the ranking of the structural error $U$ conditional on $X$. For instance, $V=0.2$ represents a typical individual whose $U$ value ranks above 20% individuals with identical covariates. $V$ enters the normalized crossing rule on the right, making an individual less likely to receive treatment; thus, it refers to resistance or distaste regarding the treatment in the MTE literature. If $V$ is large, then the propensity score $P$ should be large to induce the individual to participate in the treatment. However, an individual with a $V$ value close to zero will participate even though $P$ is small.

In the above instrument-free model, we define the MTE as the expected treatment effect conditional on the observed and unobserved characteristics:

$\Delta^{\text{MTE}}\left(x,v\right)$ captures all the treatment effect heterogeneity that is consequential for selection bias by conditioning on the orthogonal coordinates of the observable and unobservable dimensions. Given $X$ and $V$, the treatment status $D$ is fixed and thus independent of the treatment effect $Y_{1}-Y_{0}$. Similar to the MTE in the IV model, $\Delta^{\text{MTE}}\left(x,v\right)$ can be used as a building block for constructing the commonly used causal parameters, such as the average treatment effect (ATE), the treatment effect on the treated (TT), the treatment effect on the untreated (TUT), and the local average treatment effect (LATE), which can be expressed as weighted averages of $\Delta^{\text{MTE}}\left(x,v\right)$, as follows:

Somewhat surprisingly, compared with the weights on the conventional MTE (e.g., Heckman and Vytlacil, 2005, Table IB), which are generally difficult to estimate in practice, the weights on $\Delta^{\text{MTE}}\left(x,v\right)$ are simpler, more intuitive, and easier to compute.

Heckman and Vytlacil (2001, 2005) considered defining the MTE in a similar nonseparable setting by conditioning on the structural error, and showed how to integrate to generate other causal parameters. However, such an MTE cannot be identified even in the presence of IVs. Our definition, based instead on the reduced-form error, can effectively exploit the statistical independence between the observed and unobserved variables to facilitate identification of the MTE. Zhou and Xie (2019) proposed to redefine the MTE by conditioning on $P$ rather than covariates, which is a more parsimonious specification of all the relevant treatment effect heterogeneity for selection bias. The extension of our identification and estimation procedures to this alternative definition is straightforward.

Our identification strategy relies on a linearity restriction on the potential outcome equations and a nonlinearity restriction on the propensity score function. The intuition is that the propensity score minus the linear outcome index will provide the excluded variation that perturbs treatment status, which plays the role of a continuous IV. Furthermore, to ensure the exogeneity of the excluded variation, it is necessary to impose a certain form of independence between the potential outcome residuals and covariates. We denote $h_{d}\left(x\right)=E\left[\left.Y_{d}\right|X=x\right]$ and $U_{d}=Y_{d}-h_{d}\left(X\right)$, $d=0,1$, as the regression functions and residuals of potential outcomes.

Assumption CMI (Conditional Mean Independence). Assume that $E\left[U_{1}-U_{0}\left|V,X\right.\right]=E\left[U_{1}-U_{0}\left|V\right.\right]$ with probability one.

Assumption CMI is standard in the MTE literature and commonly referred to as separability or additive separability assumption (e.g., Brinch et al., 2017; Mogstad et al., 2018; Zhou and Xie, 2019), because it’s a necessary and sufficient condition for the MTE to be additively separable into observables and unobservables (Zhou and Xie, 2019, p.3074):

Namely, under Assumption CMI, the shape of the MTE curve with respect to $v$ will not vary with covariates. Assumption CMI is partly justified by $E\left[U_{d}\left|X\right.\right]=0$ and $V\perp\!\!\!\!\perp X$, which come directly from the definition. On this basis, it is sufficient for Assumption CMI to hold if the conditional covariance of $U_{d}$ and $V$ is independent of $X$, which is a key assumption in the heteroscedasticity-based identification method as well (Lewbel, 2012). Assumption CMI is implied by and much weaker than the full independence $\left(U_{d},U\right)\perp\!\!\!\!\perp X$ frequently imposed (often implicitly) in applied work, where $U$ is the structural treatment error in (1). In particular, Assumption CMI doesn’t rule out the marginal dependence of $U_{d}$ or $U$ on $X$, as illustrated by Example 3.

For the purpose of identifying the MTE, we first establish the relationship between $\Delta^{\text{MTE}}\left(x,v\right)$ and observed regression functions. Under Assumption CMI, the observed regression functions for treated and untreated individuals are additively separable in a similar way:

where

By multiplying both sides of the above equations with $p$ or $1-p$ in order to eliminate the denominators and then differentiating with respect to $p$, we obtain

where $g_{d}^{\left(1\right)}$ denotes the derivative function of $g_{d}$. Therefore, the MTE can be represented as

and the identification of it can be achieved by identifying $h_{d}$ and $g_{d}$ from the regression equation (5).

Recall that $P=\pi\left(X\right)$, where $\pi$ is a deterministic function, albeit unknown. Given that the argument of $g_{d}$ in (5) exhibits no variation apart from that generated by $X$, it is impossible to distinguish between $g_{d}$ and $h_{d}$ in the absence of additional assumptions. For example, one may choose to let $h_{d}$ absorb $g_{d}\left(\pi\left(\cdot\right)\right)$ and set $g_{d}=0$, or vice versa. Functional form restrictions that differentiate $h_{d}$ and $g_{d}\left(\pi\left(\cdot\right)\right)$ from one another can address this issue and facilitate the identification of both $h_{d}$ and $g_{d}$. To this end, we impose linearity on $h_{d}$ and nonlinearity on $\pi$.

Assumption L (Linearity). Assume that $E\left[\left.Y_{d}\right|X\right]=\alpha_{d}+X^{\prime}\beta_{d}$ with probability one for some fixed $\alpha_{d}$ and $\beta_{d}$, $d=0,1$.

Assumption L imposes a linear restriction on the potential outcome regression functions, which is a common practice in empirical studies. The linear restriction may seem too strong compared with that in the existing identification strategies, which allow highly flexible (especially nonparametric) specifications on the potential outcomes. However, for the estimation of MTE, the linear specification has been nearly universally adopted for tractability and interpretability (e.g., Kirkeboen et al., 2016; Kline and Walters, 2016; Brinch et al., 2017; Heckman et al., 2018; Mogstad et al., 2021; Aryal et al., 2022; Mountjoy, 2022). Noting that Assumption L implicitly requires the potential outcomes to be continuously distributed and supported on the entire real line, we would like to point out that our identification strategy can also be adapted to the case of limited valued outcomes by specifying linear latent index. A detailed discussion is left to Appendix C. With a little abuse of notation, we redenote $U_{d}$ to absorb the intercept $\alpha_{d}$ so that

and thus that the regression equation (5) and the MTE become

To present the nonlinearity assumption, we let $X$ be partitioned as $\left(X^{C},X^{D}\right)$, where $X^{C}$ and $X^{D}$ consist of covariates that are continuously and discretely distributed, respectively. We denote $X_{k}$, $X_{k}^{C}$, and $X_{k}^{D}$ as the $k$-th coordinates of $X$, $X^{C}$, and $X^{D}$, respectively. And we denote $x^{C}$ as a generic element in the support of $X^{C}$; likewise for $x^{D}$, $x_{k}^{C}$, and $x_{k}^{D}$. The nonlinearity assumption requires the propensity score function $\pi$ to be nonlinear in $X^{C}$ given a benchmark value of $X^{D}$. Without loss of generality, we suppose that the vector of zeros is in the support of $X^{D}$ and is the benchmark value. We denote $\pi_{0}\left(x^{C}\right)=\pi\left(x^{C},0\right)$ for notational convenience, where the discrete covariates are equal to zero.

Assumption NL (Non-Linearity). Assume that $\pi_{0}\left(x^{C}\right)$ and $E\left[Y\left|X^{C}=x^{C},X^{D}=0,D=d\right.\right]$, $d=0,1$, are differentiable with respect to $x^{C}$, and that the derivative of $\pi_{0}$ satisfies the following NL2 when $\dim\left(X^{C}\right)\geq 2$ or NL1 when $\dim\left(X^{C}\right)=1.$

– NL2 ($\dim\left(X^{C}\right)\geq 2$): there exist two vectors $x^{C}$, $\tilde{x}^{C}$ in the support of $X^{C}$ and two elements $k$, $j$ in set $\left\{1,2,\cdots,\dim\left(X^{C}\right)\right\}$ such that (i) $\partial_{k}\pi_{0}\left(x^{C}\right)\neq 0$, (ii) $\partial_{j}\pi_{0}\left(x^{C}\right)\neq 0$, (iii) $\partial_{k}\pi_{0}\left(\tilde{x}^{C}\right)\neq 0$, (iv) $\partial_{j}\pi_{0}\left(\tilde{x}^{C}\right)\neq 0$, and (v) $\left.\partial_{k}\pi_{0}\left(x^{C}\right)\right/\partial_{j}\pi_{0}\left(x^{C}\right)\neq\left.\partial_{k}\pi_{0}\left(\tilde{x}^{C}\right)\right/\partial_{j}\pi_{0}\left(\tilde{x}^{C}\right)$, where $\partial_{k}\pi_{0}\left(x^{C}\right)=\left.\partial\pi_{0}\left(x^{C}\right)\right/\partial x_{k}^{C}$ is the partial derivative of $\pi_{0}$ with respect to the $k$-th argument.

– NL1 ($\dim\left(X^{C}\right)=1$): there exists a constant $\tilde{x}^{C}$ in the support of $X^{C}$ such that $\pi_{0}^{\left(1\right)}\left(\tilde{x}^{C}\right)=0$, where $\pi_{0}^{\left(1\right)}\left(x^{C}\right)=\left.d\pi_{0}\left(x^{C}\right)\right/dx^{C}$ is the univariate derivative of $\pi_{0}$.

Assumption NL requires the propensity score function to be nonlinear in continuous covariates in a generalized sense. The combination of Assumptions L and NL will enable identification based on functional form in a semiparametric version, which can realize the identification of linear coefficients by exploiting nonlinearity elsewhere in the model. The nonlinearity assumption has two different forms, depending on the number of continuous covariates. When two or more continuous covariates are available, Assumption NL2 will require some variation in $\pi_{0}$ to distinguish it from a linear-index function. Concretely, NL2 will not hold if $\pi_{0}\left(x^{C}\right)=f\left(x^{C\prime}\gamma\right)$ for a smooth function $f$, because in this case, both sides of the inequality in (v) are equal to $\left.\gamma_{k}\right/\gamma_{j}$. Otherwise, however, it is difficult to construct examples that violate (v). In practice, (v) can be fulfilled even when $\pi_{0}$ is single-index specified, if interaction or/and quadratic terms are added, as shown in Example 4.

Assumption NL2 requires the existence of at least two continuous covariates. However, in empirical studies based on survey data, most of the demographic characteristics are documented as discrete or categorical variables, such as age, gender, race, marital status, educational attainment, and so on. Therefore, we also impose Assumption NL1, as a supplement to NL2, to take into account the situation in which only one continuous covariate is available. NL1 requires $\pi_{0}$ to have a stationary point, which implies that $\pi_{0}$ is necessarily nonlinear. NL1 will hold if the probability of receiving treatment is unaffected by some local change in the continuous covariate. It’s a differential version of the local irrelevance assumption imposed in nonseparable models for attaining point identification (e.g., Torgovitsky, 2015; D’Haultfœuille et al., 2021). NL1 can be extended to the case of no continuous covariate, as discussed in Appendix B. However, using only discrete covariates provides little identifying variation, which may lead to poor performance in the subsequent model estimation (Garlick and Hyman, 2022). Hence, we focus on the case of at least one continuous covariate.

We note that Assumption NL doesn’t exclude the widely specified linear-index treatment equation, as long as the structural error is not independent of covariates.

Unlike in the binary response model, the ubiquitous heteroscedasticity and endogeneity benefit our results while inducing no trouble, because the identification of the MTE is irrelevant to the structural coefficients in $\gamma$. Our MTE is defined by the reduced-form treatment error; thus, all we need from the treatment equation is the propensity score, which has a reduced-form nature.

Assumption NL ensures identification of coefficients of continuous covariates. In order to identify coefficients of discrete covariates, we impose a mild support condition. For any $x_{k}^{D}\neq 0$, we denote $x^{Dk}$ as the $\dim\left(X^{D}\right)\times 1$ vector with the $k$-th coordinate being equal to $x_{k}^{D}$ and all the other coordinates being equal to zero.

Assumption S (Support). For each $k\in\left\{1,2,\cdots,\dim\left(X^{D}\right)\right\}$, assume for some $x_{k}^{D}\neq 0$ in the support of $X_{k}^{D}$ that there exists $x^{C}\left(k\right)$ in the support of $X^{C}$ such that $\pi\left(x^{C}\left(k\right),x^{Dk}\right)$ is in the support of $\pi_{0}\left(X^{C}\right)$.

A sufficient condition for Assumption S to hold is that $\pi_{0}\left(X^{C}\right)$ has a full support on the unit interval, or, more generally, that the support of $\pi_{0}\left(X^{C}\right)$ is overlapped with those of $\pi\left(X^{C},x^{D}\right)$ for all $x^{D}\neq 0$. Otherwise, we have to find an $x_{k}^{D}\neq 0$ for each $k$ such that the support of $\pi\left(X^{C},x^{Dk}\right)$ is overlapped with that of $\pi_{0}\left(X^{C}\right)$. The following theorem establishes our main identification result.

The proof of Theorem 1 is grounded on the observed regression functions (7), which summarize the information from the data. As $g_{d}$ is unknown, we need to eliminate it through some subtraction to realize the identification of $\beta_{d}$. When only one continuous covariate exists, the subtraction can be carried out locally around the point satisfying Assumption NL1. Otherwise, our strategy is to perturb two continuous covariates, such as $x_{k}^{C}$ and $x_{j}^{C}$, in such a way that $\pi\left(x\right)$ remains unchanged. Specifically, for each group $d$, we increase $x_{k}^{C}$ by a small $\epsilon$ and simultaneously change $x_{j}^{C}$ by $\epsilon$ multiplied by $\left.-\partial_{k}\pi_{0}\left(x^{C}\right)\right/\partial_{j}\pi_{0}\left(x^{C}\right)$, resulting in a perturbed value of the regression function. Subtracting the perturbed regression function from the original (7) will cancel out $g_{d}$ due to the equality of $\pi\left(x\right)$, giving rise to an equation for $\beta_{d,k}^{C}$ and $\beta_{d,j}^{C}$. Note that the multiplier $\left.-\partial_{k}\pi_{0}\left(x^{C}\right)\right/\partial_{j}\pi_{0}\left(x^{C}\right)$ is the partial derivative of $x_{j}^{C}$ with respect to $x_{k}^{C}$ if we view $x_{j}^{C}$ as an implicit function of the other continuous covariates by equating $\pi_{0}\left(x^{C}\right)$ to a constant. Accordingly, Assumption NL2.(v) implies two linearly independent equations, ensuring an exact solution (i.e., identification) of $\beta_{d,k}^{C}$ and $\beta_{d,j}^{C}$. The detailed proof of Theorem 1 is presented in Appendix A. It is worth mentioning that our strategy naturally features overidentification in the sense that generally more than one pair of points in the support of $X^{C}$ satisfy Assumption NL2 or more than one point satisfies NL1, because $X^{C}$ is continuously distributed. Moreover, there may be more than one value of $X^{D}$ that satisfies Assumptions NL and S. Consequently, the identification can also be represented as the average of solutions over all points of $X^{C}$ satisfying Assumption NL and over all values of $X^{D}$ satisfying Assumptions NL and S.

According to (8), Theorem 1 implies the identification of MTE in the absence of IVs. Specifically, This result allows practitioners to include all the relevant observed characteristics into both treatment and outcome equations, without imposing any exclusion or full independence assumptions. Under Theorem 1, the conventional causal parameters can also be identified without instruments, provided that the support of $P$ contains 0 and/or 1 (which implies identifiability of $g_{d}\left(0\right)$ and/or $g_{d}\left(1\right)$):

Our identification strategy for the MTE implies a separate estimation procedure that works with the partially linear regression (7) by each treatment status. In particular, we recommend the kernel-weighted pairwise difference estimation method proposed by Ahn and Powell (1993) because of its computational simplicity, well-established asymptotic properties, and, most importantly, its capacity to effectively leverage the overidentifying information that underlies the data. Suppose that $\left\{\left(Y_{i},D_{i},X_{i}\right):i=1,2,\cdots,n\right\}$ is a random sample of observations on $\left(Y,D,X\right)$. In the first step, we estimate the nonparametrically specified propensity score using the kernel method, that is,

and

where $h_{1l}$, $l=1,2,\cdots,\dim\left(X^{C}\right)$, are bandwidths and $k_{1}$ is a univariate kernel function. If the dimension of $X^{D}$ is large, then a smoothed kernel for discrete covariates (Racine and Li, 2004) can be applied as a substitute for the indicator function, to alleviate the potential problem of inadequate observations in each data cell divided by the support of $X^{D}$. If the number of continuous covariates is not small either, then the well-known curse of dimensionality will appear, and a linear-index specification may thus be practically more relevant when modeling the propensity score. The index should include a series of interaction terms and quadratic or even higher-order terms of continuous covariates to supply sufficient nonlinear variation. The linear-index propensity score can be estimated by parametric probit/logit or semiparametric methods (e.g., Powell et al., 1989; Ichimura, 1993; Klein and Spady, 1993; Lewbel, 2000), depending on the distributional assumption on the error term.

In the second step, we estimate $\beta_{d}$ for each $d$ through a weighted pairwise difference least squares regression:

where the weights are given by

with $h_{2}$ and $k_{2}$ being the bandwidth and kernel function, respectively, which can be different from those in the first step. Given $\hat{\beta}_{d}$, the nonlinear function $g_{d}\left(p\right)$ and its derivative function $g_{d}^{\left(1\right)}\left(p\right)$ at any $p$ in the support of $P$ can be estimated by the local linear method, namely,

where

Finally, we plug $\hat{\beta}_{d}$, $\hat{g}_{d}$, $\hat{g}_{d}^{\left(1\right)}$, and $\hat{\pi}$ into identification equations (8) and (9) to estimate the MTE and other causal parameters, as follows:

The kernel-weighted pairwise difference estimator has the advantage of having a closed-form expression, so we need not solve any formidable optimization problems. However, it faces the challenging problem of bandwidth selection, similar to most semiparametric estimation methods. Alternatively, we can consider imposing a parametric specification on the unobservable heterogeneity of the MTE such that $E\left[\left.U_{d}\right|V=v\right]=E\left[\left.U_{d}\right|V=v;\theta_{d}\right]$ for finite dimensional $\theta_{d}$, e.g., the polynomial specification that $E\left[\left.U_{d}\right|V=v\right]=\sum_{j=0}^{J}\theta_{dj}v^{j}$ or the normal polynomial specification that $E\left[\left.U_{d}\right|V=v\right]=\sum_{j=0}^{J}\theta_{dj}\Phi^{-j}\left(v\right)$. In the latter, setting $J=1$ will match Heckman’s normal sample selection model. Under the parametric restriction, the second step becomes a global regression $E\left[Y\left|X,D=d\right.\right]=X^{\prime}\beta_{d}+g_{d}\left(P;\theta_{d}\right)$ with parameterized selection bias correction term $g_{d}\left(p;\theta_{d}\right)$ for each $d$; therefore, the tuning of $h_{2}$ and $h_{3}$ is circumvented. Another advantage of a parametrically specified second step is the lower computational burden compared with that of the pairwise difference estimator defined by double summation which requires a squared amount of calculation (Pan and Xie, 2023).