Semiparametric Single-Index Estimation for Average Treatment Effects

Semiparametric Single-Index Estimation for Average Treatment Effects

\def\spacingset#1{ {#1}} \spacingset{1}

\thispagestyle{empty}

{\it Keywords:} Average treatment effects; Causal inference; Hermite series expansion; Propensity score.

\pagenumbering{arabic} \spacingset{1.45}

Introduction

The literature on program evaluation has attracted significant attention and provided essential tools for empirical studies in social science. A frequently used measure for program evaluation is the average treatment effect (ATE), often estimated with restrictive parametric conditions on either the propensity score function or the outcome equation. However, when these parametric assumptions are violated, the estimator may suffer from non-negligible finite-sample bias.

There are various approaches for estimating the average treatment effects under the assumption of unconfoundedness Abadie2018. Doubly robust estimators, including the inverse probability weighting estimator augmented with additional terms robins1994a, robins1995a and the regression imputation estimator with the propensity score as an additional regressor robins1992estimating,scharfstein1999adjusting, have been widely used due to their robustness against model misspecifications. The estimators are valid as long as the specification of either the propensity score function or the outcome equation is correct. For example, bang2005doubly, robins2000a, and rotnitzky2012improved combine inverse propensity score weighting and matching methods with imputation and projection to estimate treatment effects, while the efficiency of these estimators depends on the choice of tuning parameters and computational implementation. Tan2006,Tan2010 propose a nonparametric doubly robust estimation approach and show that the estimators are efficient if the parametric propensity score model is correctly specified. However, kang2007a and vansteelandt2012on show that the doubly robust estimation approach suffers from non-trivial finite sample bias when one of the two models is misspecified, with the bias being large when both models are misspecified.

To address concerns about bias reduction, a fast-growing literature focuses on improving the robustness of doubly robust estimators. Wang2010 consider a nonparametric method to estimate the propensity score and treatment effects. However, this approach is not feasible for moderately large dimensional covariates. Vermeulen2015,Vermeulen2016 propose a novel data-driven method to reduce estimation bias under potential misspecification of both models. farrell2015a provides an inference approach on average treatment effects that is robust to model selection errors. sloczynski2017a propose a general approach for the doubly robust estimators of average treatment effects under unconfoundedness. Chernozhukov2018 provide a general method to estimate unknown functions in the nonparametric influence function and apply it to debias the estimators of average treatment effects that are feasible in high-dimensional datasets.

In this paper, we propose a semiparametric estimation method for the ATE to address potential model misspecification concerns. Our approach imposes a single-index structure on the propensity score function while allowing for a flexible link function approximated through Hermite polynomials. We simultaneously estimate the index parameter and the link function in the infinitely dimensional function space. Our method is computationally feasible and applicable to models with a large class of error distributions. We establish the model identification under specific regularity conditions, which are incorporated in our estimation procedure, and demonstrate that the proposed average treatment effect estimator achieves the parametric convergence rate and asymptotic normality. Through an extension of the propensity score residual approach in lee2017simple, we propose an estimator that remains valid even when the propensity score approaches zero or one.

This paper contributes to the existing literature on average treatment effect by introducing a novel estimator that relaxes assumptions on both propensity score function and outcome equation specification. Our methodology builds upon the work of Dong2018. In contrast to their least squares approach, we employ maximum likelihood estimation for the propensity score, which yields more efficient estimation of the propensity score function, although estimation results can be obtained using the delta method based on their least squares approach. Our estimator is computationally efficient and readily implementable in empirical research. Moreover, it provides an explicit asymptotic variance form, eliminating the need for bootstrap methods. Additionally, our estimator demonstrates robustness in datasets with limited overlap in covariate distributions and remains stable in the presence of extreme observations. Through the incorporation of a flexible assumption regarding the single index structure of the propensity score function, our estimator exhibits superior performance with numerous covariates compared to nonparametric approaches.

This paper further connects to the literature on binary outcome model estimation. Building on semiparametric single-index models, Ardakani2018 offer an innovative approach to estimating average treatment effects. Wang2010 employ a nonparametric kernel method to approximate the propensity score function; however, their estimation results deteriorate with increasing covariates and become computationally intractable for moderately large dimensional covariates due to the curse of dimensionality. Semiparametric approaches with single-index structures can maintain flexible specifications while circumventing the curse of dimensionality. sun2021estimation introduce a semiparametric single-index model for the propensity score and utilize the kernel method to estimate its functional form semiparametrically. Liu2018 establish a relationship between the treatment indicator and the low-dimensional linear structure of the covariates, estimating the propensity score function with a nonparametric link function. Both studies impose boundedness assumptions on the link function and function support, necessitating additional constraints on asymptotic theory and numerical performance.

Through relaxing the parametric assumption on the propensity score model, we estimate the conditional probability utilizing an orthogonal series-based estimation method. In comparison with the approaches of Liu2018 and sun2021estimation, our estimation methodology eliminates the need for both the boundedness of the regression function's support and the boundedness of the regression function itself Chen2007,dong2016estimation,dong2020a,LiQi2007. The method is computationally tractable and accommodates moderately large dimensional covariates. Our propensity score estimator achieves a super convergence rate $O_{p}(1/N)$ along the direction of the true parameter while maintaining a standard convergence rate $O_{p}(1/\sqrt{N})$ along all other directions, where $N$ denotes the sample size. The methodology can be further extended by assuming the conditional probability depends on the covariates vector through several linear combinations koenker2009parametric,li2016flexible,Ma2012,Ma2013,RACINE2004.

The rest of this paper is organized as follows. Section (ref) introduces treatment effect parameters of interest, and Section (ref) explains our semiparametric estimation method. We provide the large sample properties of the estimators in Section (ref) and simulation results in Section (ref). We present the empirical results regarding the effect of maternal smoking on babies' birth weight and the effect of job training programs on future earnings in Section (ref) and conclude in Section (ref). Appendix A presents the proof of the main theorem. Appendix B provides detailed sieve expansions. Appendix C contains supplementary theoretical proofs. Appendices D and E present additional simulation results and empirical findings, respectively.

Model and Identification

Setup

We consider the setup of the binary treatment and adopt the potential outcome framework proposed by rubin1974a. We use $Y(0)$ and $Y(1)$ to denote the potential outcome without and with the treatment, respectively. Also, $D$ is the treatment indicator taking 1 if an individual receives the treatment and 0 otherwise. For each unit, we observe either $Y(0)$ or $Y(1)$, and the observed outcome is given by:

Suppose that we observe a vector of covariates $X$ with the support of $\mathcal{X} \subset \mathbb{R}^{d}$ and define the propensity score function $\pi(X) := \Pr(D=1|X)$. The covariate vector $X$ does not include a constant term in its specification. The framework accommodates both continuous and discrete or categorical variables as components of $X$, and notably, we do not impose any specific distributional assumptions on the covariate structure.

As a measure of treatment effects, we consider the ATE, defined as:

Alternatively, Angrist1998ECMA discusses the variance weighted ATE, defined as:

where $\omega(X) := \operatorname{var}(D | X) / \mathbb{E}[\operatorname{var}(D | X)]$. Because $\operatorname{var}(D | X) = \pi(X) \big (1 - \pi(X) \big )$, the weight $\omega(X)$ assigns higher weights to a subpopulation if its propensity score is closer to $\frac{1}{2}$.

The propensity score function can be used to obtain consistent estimators for the average treatment effect $\Delta^{ATE}$. rosenbaum1983a show that under the assumption of unconfoundedness given covariates, the confounding bias can be removed using the propensity score function. hahn1998a proposes the regression imputation estimator that achieves the semiparametric efficiency bound for estimating average treatment effect $\Delta^{ATE}$. hirano2003a show that the inverse probability weighting approach can achieve semiparametric efficiency by estimating propensity score in a nonparametric fashion, while this approach may suffer from the curse of dimensionality in empirical analysis.

Identification of the Treatment Effect

This subsection introduces assumptions and presents an identification result. The following conditions are introduced.

(ref)(a), referred as the “unconfoundedness assumption” by rosenbaum1983a, has been widely used in the studies on the treatments effects and program evaluation dehejia1999a,heckman1998a. The validity of this assumption can be assessed by nonparametric approach rosenbaum2002a, ichino2008from. (ref)(b) ensures that the conditional probability of treatment occurrence for both treated and non-treated units is positive, usually imposed in the inverse propensity weighting approach hirano2003a,firpo2007a. rosenbaum1983a refer to the combination of these two assumptions as “strongly ignorable treatment assignment”. Since we can identify the propensity score $\pi(X)$ given the data, we can examine this assumption in practice.

For the estimation of weighted ATE, we consider the following regression model:

For the estimation of ATE, we use the following weighted regression model:

where $\upsilon(X):= \pi(X)\{1-\pi(X)\}$ is the variance of the treatment status $D$ conditional on $X$. In the above equations, $\beta$ and $\gamma$ are scalar parameters, and $U$ denotes the error term.

The variance-based weight function $\upsilon(X)$ optimally assigns weights based on estimation precision, with $\gamma$ attaining its maximum at $\pi(X) = 0.5$ and decreasing towards 0 as $\pi(X)$ approaches either 0 or 1, thereby ensuring higher weights for more precise subsamples and lower weights for less precise ones.

(ref) demonstrates that, under (ref), both the ATE and weighted ATE are identifiable from the data; the formal proof is presented in Appendix C. The ATE represents the unweighted average of individual treatment effects, whereas the weighted ATE employs differential weighting based on two criteria: estimation precision and policy relevance.\footnote{lee2017simple and Li2018JASA demonstrate that the variance-weighted ATE corresponds to the OLS estimand in the regression of the outcome on the treatment and the vector of covariates.}

Estimation

This section proposes the estimator of the treatment effect parameters explained in the previous section. Suppose that we have a sample of size equals to $N$ and observe $\{(Y_{i},X_{i},D_{i})\}_{i=1}^{N}$, which are independent copies of $(Y,X,D)$. From Lemma (ref), we can estimate $\Delta^{ATE}$ and $\Delta_{\omega}^{ATE}$, applying a least square estimation for Equations ((ref)) and ((ref)), once the propensity score is obtained. Thus, we shall focus on the estimation of the propensity score. (ref) presents a semiparametric estimation method for the propensity score function and (ref) explains a practical estimation procedure.

Estimation of the Propensity Score

To estimate the propensity score function, we consider a semiparametric approach, in which the propensity function is assumed to have a single-index structure.\footnote{The single-index structure implies a constant ratio \(\frac{\partial \pi(x) / \partial x_j}{\partial \pi(x) / \partial x_k} = \frac{\theta_{j, 0}}{\theta_{k, 0}}\), when the single index is linear in regressors. While this might appear restrictive for estimating ATE in empirical applications, our Monte Carlo studies in Appendix D demonstrate the approach's robustness to non-single index data generating processes. Furthermore, it is important to note that the single-index structure accommodates non-linearity in regressors, while maintaining linearity in index parameters.}

Specifically, the propensity score is given by

where $\Lambda: \mathbb{R} \to [0,1]$ is a known function. $g_{0}: \mathbb{R} \to \mathbb{R}$ is an unknown link function and $\theta_0$ is a $d \times 1$ vector of unknown parameters. The mapping $\Lambda(\cdot)$ ensures that the propensity score function remains within the unit interval. We employ the logistic function in our analysis, specifically $\Lambda(z)= e^{z}/(1+e^{z})$ for $z \in \mathbb{R}$. Our theoretical results generalize to alternative functions.

The existing literature on single-index models typically imposes restrictive assumptions regarding the boundedness of either the support or the link function. For instance, Ichimura1993 requires compactness on the parameter space and the support of regressors, which consequently restricts the link function to be defined on a compact set. Similarly, Xia2007 and Cai2014 stipulate that the link function be bounded, despite being defined on the entire real line. More recent work by Dong2015 necessitates that the link function be smooth and integrable, thereby implying its boundedness.

Our methodology advances the literature by accommodating an unbounded link function with unbounded support. By conceptualizing the link function $g_{0}$ as a point within an infinitely dimensional function space, we facilitate the simultaneous estimation of the index parameter and the link function. Importantly, the Hilbert space $\mathcal{L}^{2}$ encompasses a wide range of functions---including polynomials, power functions, and bounded functions on $\mathbb{R}$---that are prevalent in both practical applications and econometric theory.

We assume that the true parameter vector $\theta_{0}$ is an element of the compact parameter space $\Theta \subset \mathbb{R}^{d}$ and impose the normalization constraint $ \|\theta_{0}\|=1$ along with $\theta_{0,1}>0$ for identification purposes. The constraint $\theta_{0,1}>0$ serves to pin down the direction of $\theta_{0}$, where the estimation achieves super-convergence rate $O_P(N^{-1})$ in the direction of $\theta_{0}$, while maintaining root-N convergence rate along all directions orthogonal to $\theta_{0}$. This normalization restriction does not affect the optimization procedure under the identification constraint when deriving the index parameter estimator.

The model permits general forms of heteroskedasticity in the function form of the propensity score. The semiparametric approach with a single-index constraint can retain flexible specifications while avoiding the curse of dimensionality. Our method relies on the optimization with the constraint of the identification condition to derive the estimator of parameter $\theta_{0}$. Compared with the semiparametric approach using nonparametric kernel method Liu2018,sun2021estimation, our estimator does not impose restrictions on the link function or its support and establishes the asymptotic theory for inference.

We consider the case where the nonparametric function $g_{0}(\cdot)$ resides in the Hilbert space $\mathcal{L}^2$ with an inner product $\langle g, h \rangle := \int g(w) h(w) \exp (-w^{2} / 2) d w$ for $g, h \in \mathcal{L}^{2}$. The Hilbert space contains all polynomials, all power functions and all bounded, real-valued functions, which are often encountered in both applications and econometric theory Chen2007.\footnote{We explain the properties of Hilbert space used for this paper in Appendix B.} Let $\{h_{j}(\cdot) \}_{j=1}^{\infty}$ be the orthonormal basis for the Hilbert space. Then, we have the orthogonal series expansion for any function $g \in \mathcal{L}^2$:

where the inner product $c_{j}:=\langle g, h_{j}\rangle$. In practice, one has to choose a truncation parameter $k$ to approximate the infinite series, such that:

where $ \widetilde{g}_{k}(\omega) :=\sum_{j=1}^{k} c_{j} h_{j}(\omega) $ and $\epsilon_{k}(\omega) :=\sum_{j=k+1}^{\infty} c_{j} h_{j}(\omega) $.

By virtue of (ref), the $g_0(X'\theta_{0}) $ in (ref) can have the following representation:

We have the maximum likelihood framework and define the log-likelihood function as follows:

For the estimation of the single-index model, we consider the maximum likelihood framework and define the log-likelihood function as follows:

The estimators $\hat{\theta}$ and $\{\hat{c}_{j}\}_{j=1}^{k}$ of $\theta$ and $\{c_{j}\}_{j=1}^{k}$ are defined as the solution to the following maximization problem:

Letting $\widehat{g}_{k}(\cdot):= \sum_{j=1}^{k} \hat{c}_{j} h_{j}(\cdot)$, we have an estimator of the propensity score function:

Our propensity score estimation method is grounded in the work of Dong2018, which explores a semi-parametric single-index model wherein the link function is permitted to be unbounded and possess unbounded support. In this model, the link function is treated as a point in an infinitely-dimensional function space, and both the index parameter and link function can be estimated simultaneously through an optimization process, subject to the identification condition constraint for the index parameter.

While the least squares approach in Dong2018 is prevalent, it may not consistently produce efficient estimates, particularly when the underlying linear regression assumptions are not met. Our maximum likelihood estimation (MLE) approach offers several benefits, including a more versatile estimation framework suitable for various statistical models, encompassing those with non-normal error distributions. This adaptability enables the accommodation of different propensity score model types. Moreover, the MLE method boasts favorable asymptotic properties, resulting in more efficient propensity score function estimation, especially in large samples. Our approach can also be expanded to include prior knowledge through Bayesian estimation techniques, further enhancing efficiency and robustness.

Our estimator has nice finite sample properties. First, compared with a parametric approach, such as probit and logit, our propensity score estimator is more flexible and less prone to misspecification. Second, the estimator can be used for cases with multiple regressors and avoids the curse of dimensionality commonly encountered in nonparametric methods hahn1998a,hirano2003a. Finally, we do not impose any distributional assumptions on the covariates by including continuous, discrete, and categorical variables.

Estimation Procedure

This subsection provides an estimation procedure for the nonlinear constrained optimization problem in Equation ((ref)). Given a truncation parameter $k$, we obtain the estimators $\hat{\theta}$ and $\{\hat{c}_{j}\}_{j=1}^{k}$ as follows:

Algorithm 1. For a fixed truncation parameter $k$,

(i) Use a linear regression estimator as the initial estimator { $\check \theta := (\sum_{i=1}^{N} X_i X_i^{\prime} )^{-1} \sum_{i=1}^{N} X_iD_i.$}

(ii) Obtain the estimator $\{\hat{c}_{j}\}_{j=1}^{k}$ as the solution to the problem: { $\max_{\{c_{j}\}_{j=1}^{k}} \ell_{N}\big (\check{\theta}, \{c_{j}\}_{j=1}^{k} \big)$.}

(iii) Given the estimator $\{\hat{c}_{j}\}_{j=1}^{k}$, obtain the estimator $\hat{\theta}$ as the solution to the constrained optimization problem:

where $\lambda$ is a scalar Lagrange multiplier.

While existing works examine the choice of optimal truncation parameter for non-parametric models Gao2002, there is limited theoretical research on the optimal choice of $k$ for the single-index model. For addressing the choice of the truncation parameter, we propose a leave-one-out cross-validation algorithm. Notably, the truncation error becomes asymptotically negligible under the regularity conditions Dong2018.

Algorithm 2. Let $\mathcal{K}$ be the set of candidate values for the truncation parameter.

(i) For each $k \in \mathcal{K}$, we apply Algorithm 1 to the dataset $\{(Y_{i},X_{i},D_{i})\}$ excluding the $j$-th observation and obtain the leave-one-out estimate, denoted by $\widehat{Y}_{k, -j}$.

(ii) The optimal truncation parameter is determined by minimizing the prediction mean of squares:

To demonstrate the empirical validity of our approach, we present simulation results examining the finite sample performance of the estimation method with various choices of truncation parameter $k$ selected from a set of suitable candidate values in Appendix D.4.

Asymptotic Properties

To establish asymptotic properties for the average treatment effect estimator, we assume the following conditions and provide justifications for the following assumptions.

(ref)(b) covers the conditions that are commonly used in the existing literature on sieve estimation Chen2007 and binary choice model coppejans2001estimation,cosslett1983distribution,Ichimura1993,Klein1993 as a special case and includes many commonly used distributions such as normal distribution and those possessing compact support. Recall we impose $\|\theta_0\|=1$ for identification purpose. The estimator is derived from an optimization with the constraints for index parameter, and we do not need to impose additional parameter normalization conditions.\footnote{ma2016semiparametric develops a semiparametric single index estimator for the propensity score with a one-step maximum likelihood-type estimator that imposes parameter normalizations. Our approach relax such restrictions on parameters including function form $g$ and parameter $\theta$ as long as the identification condition $\|\theta_0\|=1$ with $\theta_{0,1}>0$ holds.} (ref)(c) imposes the smoothness restrictions on the link function and the divergence restrictions on the truncation parameter to guarantee the truncation error is negligible.

For theoretical developments, we need to assume additional regularity conditions. Since these conditions are standard in the literature, we state the conditions as (ref).

(ref)(a) covers commonly used conditions in the existing literature and establishes the uniform convergence of $(\theta_{0}, g_{0})$ in a small neighborhood. The uniform convergence of $(\theta_{0}, g_{0})$ in a small neighborhood is required by yu2002penalized, which can be derived using the results of Lemma A2 in Newey2003. For simplicity, we directly impose conditions (ii) and (iii) in (ref)(a). (ref)(b) matches Assumption 2 of newey1997convergence and Assumption 3 of Su2012, and the restrictions on the derivatives of orthogonal Hermite functions can be derived using the results in belloni2015some. (ref)(c) is the under-smoothing condition belloni2015some,chang2015high that ensures the negligibility of the truncation error terms and the asymptotic normality of the estimator.

(ref) imposes restrictions on the probability distribution functions of the regressors $h_{m}(\omega)$ and $h_{n}(\omega)$ to make the functions reduce to $k^{2d}/N$ as we use Hermite polynomials to decompose the link function $g_0(\omega)$ in the function space $\mathcal{L}^{2}$.

The subsequent lemma establishes the consistency and asymptotic properties of the estimators $\hat{\theta}$. The proof of (ref) is presented in Appendix C.

In (ref) below, we establish the convergence rate and the limit distribution of the estimator $\hat{\theta}$. To present the result, let $P_{\theta_{0}} := I_{d}-\theta_{0}\theta_{0}^{\prime}$ be a $d \times d$ matrix. Under the constraint of $\|\theta_{0}\|=1$, $P_{\theta_{0}}$ represents the projection matrix that maps any $d \times 1$ vector into the orthogonal complement space of the true parameter $\theta_0$. Let $V:=[v_{1}, \dots, v_{d-1}]$ be a $d \times (d-1)$ matrix consisting of eigenvectors $v_{j}$ associated with the eigenvalue of 1.

(ref)(a) demonstrates the asymptotic normality of $V^{\prime}(\hat{\theta}-\theta_{0})$ through the transformation of $\hat{\theta}-\theta_{0}$ into $\mathbb{R}^{d-1}$, which is necessary since both $\hat{\theta}$ and $\theta_{0}$ are constrained to the unit sphere in $\mathbb{R}^{d}$. (ref)(b), which follows from the constraint $\|\theta_{0}\|=1$, characterizes the convergence behavior along the radial direction. Specifically, it reveals that the convergence rate along the direction of $\theta_{0}$ is faster ($O_p(1/N)$) compared to the directions orthogonal to $\theta_{0}$ ($O_p(1/\sqrt{N})$).

We are now ready to establish the large sample properties of the estimators for ATE ($\widehat{\Delta}^{ATE}$) and variance weighted ATE ($\widehat{\Delta}_w^{ATE}$) defined in (ref) where the propensity score is estimated based on the series estimation approach in (ref). The proof of (ref) below is given in Appendix A.

The variance weighted ATE estimator $\widehat{\Delta}_w^{ATE}$ has the form:

The ATE estimator $\widehat{\Delta}^{ATE}$ is:

Our proposed estimator has several advantages over other estimators including the inverse propensity score weighting hirano2003a,imbens2010a, matching abadie2006a,abadie2011a,abadie2016a, regression adjustment lane1982analysis, and doubly robust estimators Tan2010,tsiatis2006semiparametric,Wooldridge2007,Liu2018. First, our estimator is computationally simple and easy to implement. There is no need to specify the propensity score function form. Second, our estimator has a simple explicit form of asymptotic variance that works well even in small samples, and we do not need to estimate the asymptotic variance through bootstrap procedures. Third, compared with the inverse propensity score approach, our estimation method works well in the datasets that are limited overlap in the covariate distributions and is not influenced by extreme observations in empirical applications.\footnote{Wooldridge2007 proposes the inverse probability weighted M-estimation under a general missing data scheme and shows the misspecication of the error distribution is negligible under some regularity conditions, including the existence of pseudo-parameters. Our approach relaxes the requirement and our simulation and empirical results, reported in Sections (ref) and (ref), show some advantages of our approach in finite samples.} Fourth, by adopting a flexible assumption on the single index structure of propensity score function, our estimator performs well when there are many covariates.

For statistical inferences, we need to estimate the asymptotic variances $\sigma_{\omega}^2$ and $\sigma^2$. Following the sample analog principles, we can estimate the asymptotic variances by using the sample moments based on the observables: $\hat{U}_{i}$, $D_{i} - \widehat{\pi}(X_{i})$, and $\widehat{\upsilon}(X_{i}):= \widehat{\pi}(X_{i})\{1-\widehat{\pi}(X_{i})\}$. We denote by $\hat{\sigma}^2_{\omega}$ and $\hat{\sigma}^2$ the estimator of $\sigma_{\omega}^2$ and $\sigma^2$, respectively. The proposition below shows that those estimators are consistent. The proof of (ref) is given in Appendix C.

Monte Carlo Studies

In this section, we conducted comprehensive Monte Carlo simulations to evaluate our semiparametric estimator proposed in (ref) for the propensity score function and semiparametric estimator proposed in (ref) for the average treatment effect. We repeat each experiment 10,000 times with sample size $N$ = 400, 800, or 1600 and set the truncation parameter $k$ equal to $\lfloor N^{1 / 5}\rfloor$.\footnote{We provide an algorithm for calculating truncation parameter and our simulation results are valid and robust to a range of truncation parameters, as demonstrated in Appendix D.4.}

Propensity Score Function

We conducted Monte Carlo simulations to evaluate our semiparametric estimator proposed in (ref) for the propensity score function, with detailed simulation setup and results in Appendix D.1. Our first simulation examined four functional forms with two-dimensional covariates $X \sim N(0, I_{2})$, considering both cases where $\|\theta_{0}\|=1$ holds and where it is violated. As shown in (ref), our estimator demonstrates small biases that decrease with sample size, achieving the theoretical $O(1/N)$ convergence rate when $\|\theta_{0}\|=1$ and $O(1/\sqrt{N})$ otherwise. We then extended to six-dimensional covariates $X \sim N(0, I_{6})$. The results in (ref) confirm our estimation approach remains valid with higher dimensions, maintaining small biases and decreasing RMSEs as sample size increases. Finally, we examined the estimator's performance for average treatment effects under various challenging scenarios, including normalization condition violations, probit specifications, heavy-tailed errors, and propensity score misspecification. Results in (ref) demonstrate the robustness of our approach, with small biases and RMSEs maintained across different sample sizes.

We further compare our semiparametric (SP) estimator with existing approaches including covariate balancing (CB) and high-dimensional selection (HD) methods for propensity score estimation, with detailed simulation setup and results in Appendix D.2. Our simulation study considers four challenging DGPs featuring unbounded link functions and unbounded support for propensity score functions, including both single-index and non-single-index structures. As shown in (ref), while all three estimators demonstrate near-zero bias, our SP estimator outperforms CB Anna2021JAE,imai2014covariate,Zubizarreta and HD belloni2017a,Chernozhukov2018,SunTan2021 approaches in terms of standard deviation and RMSE, particularly for non-single-index DGPs.

Average Treatment Effect

We conduct an extensive simulation study for the finite sample properties of our semiparametric estimator proposed in (ref) for the average treatment effect.

{\bf Setting 5A}: $g_0(X^{\prime}\theta_0) = \sin(X^{\prime}\theta_0)$, where $\theta_0 = (\theta_{0,1}, \theta_{0,2})^{\prime} = (0.8,-0.6)^{\prime}$ and the vector of covariates $X=(X_{1}, X_{2})^{\prime}$ are generated from independent standard normal distributions $N(0,1)$. The propensity score function is defined as $D = \Lambda \big( g_0(X'\theta_{0}) \big)$ and the outcome model is generated as $Y=\beta_d D + X_{1}+ X_{2} + \epsilon$, where $\epsilon \sim N(0,1)$ and $\beta_d = 1$.

{\bf Setting 5B}: $g_0(X^{\prime}\theta_0) = 0.5\{(X^{\prime}\theta_0)^3 - (X^{\prime}\theta_0)\}$, where $\theta_0 = (\theta_{0,1}, \theta_{0,2})^{\prime} = (0.8,-0.6)^{\prime}$ and the vector of covariates $X=(X_{1}, X_{2})^{\prime}$ are generated from independent standard normal distributions $N(0,1)$. The propensity score function is defined as $D = \Lambda \big( g_0(X'\theta_{0}) \big)$ and the outcome model is generated as $Y=\beta_d D + X_{1}+ X_{2} + \epsilon$, where $\epsilon \sim N(0,1)$ and $\beta_d = 1$.

{\bf Setting 5C}: $g_0(X^{\prime}\theta_0) = \sin(X^{\prime}\theta_0)$, \\where $\theta_0 = (\theta_{0,1}, \theta_{0,2}, \theta_{0,3}, \theta_{0,4}, \theta_{0,5}, \theta_{0,6})^{\prime} = (\sqrt{0.2},\sqrt{0.3},\sqrt{0.25},-\sqrt{0.1},\sqrt{0.08},-\sqrt{0.07})^{\prime}$ and the vector of covariates $X=(X_{1}, X_{2},X_{3}, X_{4},X_{5}, X_{6})^{\prime}$ are generated from independent standard normal distributions $N(0,1)$. The propensity score function is defined as $D = \Lambda \big( g_0(X'\theta_{0}) \big)$ and the outcome model is generated as $Y=\beta_d D + X_{1}+ X_{2} + X_{3} + X_{4} + X_{5} + X_{6} + \epsilon$, where $\epsilon \sim N(0,1)$ and $\beta_d = 1$.

{\bf Setting 5D}: $g_0(X^{\prime}\theta_0) = 0.5\{(X^{\prime}\theta_0)^3 - (X^{\prime}\theta_0)\}$, \\where $\theta_0 = (\theta_{0,1}, \theta_{0,2}, \theta_{0,3}, \theta_{0,4}, \theta_{0,5}, \theta_{0,6})^{\prime} = (\sqrt{0.2},\sqrt{0.3},\sqrt{0.25},-\sqrt{0.1},\sqrt{0.08},-\sqrt{0.07})^{\prime}$ and the vector of covariates $X=(X_{1}, X_{2},X_{3}, X_{4},X_{5}, X_{6})^{\prime}$ are generated from independent standard normal distributions $N(0,1)$. The propensity score function is defined as $D = \Lambda \big( g_0(X'\theta_{0}) \big)$ and the outcome model is generated as $Y=\beta_d D + X_{1}+ X_{2} + X_{3} + X_{4} + X_{5} + X_{6} + \epsilon$, where $\epsilon \sim N(0,1)$ and $\beta_d = 1$.

We summarize the simulation results in (ref). In all four simulation settings, our estimators perform well according to bias and standard deviation criteria. For all the sample sizes under consideration, the estimators have relative small biases and standard deviations. The finite--sample simulation results support that the asymptotic normal approximation is accurate and the rate of convergence is consistent with the results in (ref).

{

}

We compare our semiparametric (SP) estimator with the locally efficient estimator Liu2018 and nonparametric cross-validation estimator sun2021estimation for treatment effects estimation, with detailed simulation setup and results in Appendix D.3. While these alternative approaches assume bounded link functions and support, we examine performance under four challenging DGPs with unbounded properties. As shown in (ref), our SP estimator demonstrates smaller bias when boundedness assumptions are violated, and achieves lower standard deviation and RMSE across all sample sizes compared to both alternatives. The results confirm our method's superior performance for data with unbounded link functions and support.

We investigate the sensitivity of our estimator to different truncation parameters $k$, with detailed simulation setup and results in Appendix D.4. Results in (ref) demonstrate the performance of propensity score estimation across truncation parameters ranging from $k = 2$ to $k = 6$. Similarly, (ref) show the average treatment effect estimation results across the same range of truncation parameters. The results confirm that our estimators maintain validity across different choices of truncation parameters.

Empirical Study

In this section, we apply the proposed semiparametric method to analyze two real data examples. We study the effect of maternal smoking on babies' birth weight in the main context and consider the effects of the job training program on future earnings in Appendix E.\footnote{As shown in cattaneo2015a, the assumption that maternal smoking is exogenous to the babies' birth outcome may not hold. We apply our approach for this empirical study as it is widely used in treatment effects studies as a benchmark. We also provide the effects of the job training program on future earnings in Appendix E as an additional empirical study.} In both cases, we demonstrate the validity of our estimation and inference method for the datasets with potential misspecification of propensity score and datasets with limited overlap.

We apply our semiparametric method to analyze the ATE of maternal smoking on babies' birth weight. Low birth weight may increase infant mortality rates and economic costs, including late entry into kindergarten, repeated grades, and longer-term labor market outcomes arcidiacono2011a,Bao2021JFQA,permutt1989simultaneous,rosenzweig1991inequality.

Although a large body of research confirms the negative effect of maternal smoking on babies' birth weight, there is no consensus on its exact magnitude abrevaya2001the,abrevaya2006estimating,chernozhukov2011inference,evans1999can. We use the dataset of low birth weight initially by Almond2005. This is a rich database of singletons in Pennsylvania with 4,642 detailed observations of mothers and their infants' birth information.

The outcome variable $Y$ is infant birth weight measured in grams. The binary treatment variable $D$ is the mother's smoking status ($D = 1$ indicates mother is a smoker, $D = 0$ indicates mother is a non-smoker). The covariates $X$ include mother's age, mother's marital status, an indicator variable for alcohol consumption during pregnancy, an indicator for whether there was a previous birth where the newborn died, mother's education, father's education, number of prenatal care visits, mother's race, an indicator of firstborn baby, and months since last birth cattaneo2010a.

The critical assumption in our application is that the mother's smoking status is independent of the infant's birth weight conditional on all observed demographic variables, implying that maternal smoking may impact the babies' birth weights only through its effect on observed covariates. We show the summary statistics for all variables in Table (ref), which provides the summary statistics for baby birth weight data, where $N$ denotes sample size. For each variable, we report the sample average (Mean) and sample standard deviation (Std). *, **, and *** indicate the significance level at 10%, 5%, and 1%, respectively. As it reveals, the comparison between these two groups shows that the control group is quite different from the treated group.

(ref) shows the estimation of the average treatment effect based on our approach. The naive difference in the weight of babies belonging to non-smoking and smoking mothers is 275.3 grams. Given the potential confounding effects of covariates on the potential outcome, this result may not be a valid estimate of the average treatment effect. We next compare the results of average treatment effect estimation using the augmented inverse propensity weighting (AIPW) estimator Tan2010,tsiatis2006semiparametric, inverse propensity weighting (IPW) estimator hirano2003a,imbens2010a, inverse propensity weighting with regression adjustment (IPW-RA) estimator Wooldridge2007, bias-corrected matching (MBC) estimator abadie2006a,abadie2011a, propensity score matching (PSM) estimator rosenbaum1983a,abadie2016a, regression adjustment (RA) estimator lane1982analysis, doubly robust with efficient propensity score estimation (Efficient), locally efficient propensity score estimation (Local), or parametric logistic estimation (Logistic) estimators Liu2018,Ma2013. We summarize the above estimation results for the average treatment effect of maternal smoking on babies' birth weight in (ref) with the mean, standard deviation, $Z$-statistics, $p$-value, and 95% confidence interval.

We note that the doubly robust estimator using the propensity score estimated by logistic regression is different from other estimators, suggesting that the parametric logistic form may not be suitable for the propensity score function in this data example. In comparison, our estimation approach is valid for potential misspecification of the propensity score function and provides a valid estimation and inference for the average treatment effect.

To further illustrate the potential misspecification of propensity score, we plot the estimated link function $\hat{g}$ of propensity score using the wild bootstrap simulation with 500 bootstrap repetitions. As we use the logistic function $\Lambda$ in (ref), if the true DGP for propensity score in this data example is logistic, then the estimated link function $\hat{g}$ should be the identity function. We calculate the optimal truncation parameter $k$ based on the minimization of the prediction mean of squares. In this data example, the optimal truncation parameter $k = 2$. We estimate the approximated second-order polynomial fitted function using the ordinary least squares estimation approach as follows (with standard errors in brackets):

where $\omega_{i}$ is the $x_{i}'\hat{\theta}$ for $1\leq i \leq N$.

{

}

We show the estimation in (ref) with the estimated link function in solid line and the corresponding 95% confidence interval in dashed line. We also plot the approximated second-order polynomial fitted function using the ordinary least squares estimation approach as a reference. Our results show that the estimated link function $\hat{g}$ does not seem to be an identity function. Such potential misspecification of the propensity score may explain the differences in treatment effect estimation between our estimator and the doubly robust estimator using the propensity score estimated by logistic regression. Both the estimated link function and the approximated second-order polynomial fitted function are close to each other, demonstrating that the estimation method works well in this empirical application.

Conclusion

We propose a semiparametric estimation method for determining the average treatment effect. Our approach assumes a single-index structure for the propensity score function, allowing for a flexible link function. We simultaneously estimate the index parameter and link function in an infinitely dimensional function space, and calculate the propensity score through optimization, taking into consideration the constraints of identification conditions. We also establish the large sample properties of our estimator. Our method is empirically applicable and suitable for models encompassing a wide range of error distributions and offers computational flexibility for problems of moderately large dimensions. We conduct an extensive simulation study to assess the finite sample performance of our proposed estimator in finite samples. We present empirical results concerning the impact of maternal smoking on infant birth weight, as well as the effect of job training programs on future earnings. Our approach is less susceptible to misspecification of the propensity score function and remains valid when the propensity score approaches zero or one.

{

\pagenumbering{arabic}

Appendix A provides the proof for (ref). Appendix B gives some preliminary results about Hilbert space, Appendix C gives the proofs of the main results, Appendix D provides some additional numerical results, and Appendix E presents additional empirical application results.

\setcounter{equation}{0} {A.\arabic{equation}} {A.\arabic{subsection}} {A.\arabic{figure}} {A.\arabic{table}} {A.\arabic{lemma}} {A.\arabic{remark}} {A.\arabic{corollary}} {A.\arabic{assumption}} {A.\arabic{proposition}}

\setcounter{equation}{0} {B.\arabic{equation}} {B.\arabic{subsection}} {B.\arabic{figure}} {B.\arabic{table}} {B.\arabic{lemma}} {B.\arabic{remark}} {B.\arabic{corollary}} {B.\arabic{assumption}} {B.\arabic{proposition}}

\setcounter{equation}{0} {C.\arabic{equation}} {C.\arabic{subsection}} {C.\arabic{figure}} {C.\arabic{table}} {C.\arabic{lemma}} {C.\arabic{remark}} {C.\arabic{corollary}} {C.\arabic{assumption}} {C.\arabic{proposition}}

\setcounter{equation}{0} {D.\arabic{equation}} {D.\arabic{subsection}} {D.\arabic{figure}} {D.\arabic{table}} {D.\arabic{lemma}} {D.\arabic{remark}} {D.\arabic{corollary}} {D.\arabic{assumption}} {D.\arabic{proposition}}

\setcounter{equation}{0} {E.\arabic{equation}} {E.\arabic{subsection}} {E.\arabic{figure}} {E.\arabic{table}} {E.\arabic{lemma}} {E.\arabic{remark}} {E.\arabic{corollary}} {E.\arabic{assumption}} {E.\arabic{proposition}}

{\bf Appendix A: Proof of (ref)}

\pagenumbering{arabic}

{\bf Appendix B: Preliminary Properties for Hilbert Space}

In this section, we present the properties of the Hilbert Space $\mathcal{L}^{2}$ that includes all polynomials, all power functions, all bounded functions, and even some exponential functions. More importantly, all functions in the space are defined on the entire real line.

An inner product for $g_{1}(w), g_{2}(w) \in \mathcal{L}^{2}$ is given by: $$ \langle g_{1}(w), g_{2}(w)\rangle=\int g_{1}(w) g_{2}(w) \exp(-w^{2} / 2) d w. $$

Hermite polynomials form a complete orthogonal sequence in $\mathcal{L}^{2}$ with each elements defined by:

The orthogonality of this basis system satisfies $\int h_{m}(w) h_{n}(w) \exp (-w^{2} / 2) d w=\sqrt{2 \pi} \delta_{m n}$, where $\delta_{m n}$ is the Kronecker delta. For $\forall g(w) \in \mathcal{L}^{2}$, we have orthogonal series expansion: $$ g(w)=\sum_{j=1}^{\infty} c_{j} h_{j}(w), \quad c_{j}=\langle g(w), h_{j}(w)\rangle. $$

For any truncation parameter $k \geq 1$, we split the orthogonal series expansion into two parts as following:

where $g_{k}(w) :=H(w)^{\prime} C_{k}$, $H(w)=\big(h_{1}(w), \ldots, h_{k}(w)\big)^{\prime}$, $C_{k}=(c_{1}, \ldots, c_{k})^{\prime}$, and the approximation residual is given by $\epsilon_{k}(w)=\sum_{j=k+1}^{\infty} c_{j} h_{j}(w)$.

The quantity of $\|H(w)\|$ is crucial for the asymptotic properties of our estimators, using the results in the Theorem 1.1 of levin1992christoffel and dong2020a, we have:

uniformly for $k \geq 1$ and $w \in\{u:|u| \leq \sqrt{2 k}(1+L k^{-2 / 3})\},$ where $L>0$ is given.\footnote{In our context, $f \asymp g \Longleftrightarrow \exists C, D>0: C|g| \leq|f| \leq D|g|$.} In addition, we can show that the order of $\|H(w)\|$ is $O(\sqrt{k})$ when the condition for $1-|x| / \sqrt{2 k}>k^{-2 / 3}$ satisfies as long as the truncation parameter $k$ diverge.

\setcounter{equation}{0}

{\bf Appendix C: Proofs of Main Results}

\setcounter{equation}{0} \setcounter{table}{0}

{\bf Appendix D: Additional Simulation Results}

{\bf Appendix D.1: Estimating the Propensity Score Function}

We first conduct a Monte Carlo study to investigate the finite sample properties of our semiparametric estimator proposed in (ref) for the propensity score. We consider the following four functional forms (DGPs) with a vector of regressors $X=(X_{1}, X_{2})^{\prime} \sim N(0, I_{2})$:

The propensity score function is defined as $D = \Lambda \big( g_0(X'\theta_{0}) \big)$. The DGPs in Settings 1A and 1B satisfy the $\|\theta_{0}\|=1$ for identification purpose, while the DGPs in Settings 2A and 2B violate $\|\theta_{0}\|=1$.

Let $\hat{\theta}= (\hat{\theta}_{1}, \hat{\theta}_{2})'$ be the estimator of $\theta_{0}$. We evaluate the estimators in terms of bias, standard deviation (std) and root mean-squared error (RMSE). We summarize the simulation results in (ref). For the DGPs in Settings 1A and 1B that satisfy $\|\theta_{0}\|=1$, the simulation results show that our estimator performs well. Moreover, by focusing on the bias and RMSE separately, it is clear that both terms vanish quickly as the sample size increases. The asymptotic normal approximation is accurate as the standard deviations become smaller as the sample size grows. Our estimation approach is valid for the polynomial single index structure or more complicated periodic single index form. More importantly, the convergence speed is $O(1/N)$ along the direction of true parameter $\theta_{0}$, consistent with the results in (ref).

{

}

For the DGPs in Settings 2A and 2B that violate $\|\theta_{0}\|=1$, the estimated value $\hat{\theta}_{1}$ and $\hat{\theta}_{2}$ converge to $1/\sqrt{2}$ and -$1/\sqrt{2}$, respectively. For all sample sizes we studied, the estimators have relative small biases. The standard deviations of $\hat{\theta}_{1}$ and $\hat{\theta}_{2}$ become smaller as the sample size increases. The asymptotic normal approximation is accurate and the convergence speed is $O(1/\sqrt{N})$ along all other direction orthogonal to the true value $\theta_{0}$, in accordance with the results in (ref).

We further conduct the simulation study to investigate the performance of our estimators for six dimensions of covariates. We consider the following four DGPs for the propensity score with the vector of regressors $X=(X_{1}, X_{2},X_{3},X_{4},X_{5}, X_{6})^{\prime} \sim N(0,I_{6})$:

The propensity score function is defined as $D = \Lambda \big( g_0(X'\theta_{0}) \big)$. The DGPs in Settings 3A and 3B satisfy $\|\theta_{0}\|=1$, while the DGPs in Settings 4A and 4B violate $\|\theta_{0}\|=1$.

We show the simulation results for Settings 3A and 3B in (ref) and the results for Settings 4A and 4B in (ref), respectively. Our estimation approach is valid for the six dimensions of covariates we considered. All estimators have small biases for the three sample sizes we considered. The biases vanish quickly when the sample size increases and RMSEs generally decrease as the sample size increases. The asymptotic normal approximation is accurate as the standard deviations become smaller as the sample size grows. The convergence speed is consistent with the results in (ref).

{

}

We further conduct a Monte Carlo study to investigate the finite sample properties of our semiparametric estimator proposed in (ref), and examine the average treatment effect under potential misspecification of propensity score function and the relaxation of $\|\theta_{0}\|=1$ in (ref). We consider the following eight DGPs for both the propensity score and potential outcome.

{\bf Setting 7A}: $g_0(X^{\prime}\theta_0) = \sin(X^{\prime}\theta_0)$, where $\theta_0 = (\theta_{0,1}, \theta_{0,2})^{\prime} = (0.5,-0.5)^{\prime}$, the vector of covariates $X=(X_{1}, X_{2})^{\prime}$ are generated from standard normal distributions $N(0,1)$. The propensity score function is defined as $D = \Lambda \big( g_0(X'\theta_{0}) \big)$ and the outcome model is generated as $Y=\beta_d D + X_{1}+ X_{2} + \epsilon$, where $\epsilon \sim N(0,1)$ and $\beta_d = 1$.

{\bf Setting 7B}: $g_0(X^{\prime}\theta_0) = 0.5\{(X^{\prime}\theta_0)^3 - (X^{\prime}\theta_0)\}$, where $\theta_0 = (\theta_{0,1}, \theta_{0,2})^{\prime} = (0.5,-0.5)^{\prime}$, the vector of covariates $X=(X_{1}, X_{2})^{\prime}$ are generated from standard normal distributions $N(0,1)$. The propensity score function is defined as $D = \Lambda \big( g_0(X'\theta_{0}) \big)$ and the outcome model is generated as $Y=\beta_d D + X_{1}+ X_{2} + \epsilon$, where $\epsilon \sim N(0,1)$ and $\beta_d = 1$.

{\bf Setting 8A}: $D=1\{X^{\prime}\theta_0+\varepsilon>0\}$, where $\varepsilon \sim N(0,1)$ is independent of $X$, $\theta_0 = (\theta_{0,1}, \theta_{0,2})^{\prime} = (0.8,-0.6)^{\prime}$, the vector of covariates $X=(X_{1}, X_{2})^{\prime}$ are generated from standard normal distributions $N(0,1)$. The outcome model is generated as $Y=\beta_d D + X_{1}+ X_{2} + \epsilon$, where $\epsilon \sim N(0,1)$ and $\beta_d = 1$.

{\bf Setting 8B}: $D=1\{(X^{\prime}\theta_0)^3 + (X^{\prime}\theta_0)^2+ X^{\prime}\theta_0 + \varepsilon>0\}$, where $\varepsilon \sim N(0,1)$ is independent of $X$, $\theta_0 = (\theta_{0,1}, \theta_{0,2})^{\prime} = (0.8,-0.6)^{\prime}$, the vector of covariates $X=(X_{1}, X_{2})^{\prime}$ are generated from standard normal distributions $N(0,1)$. The outcome model is generated as $Y=\beta_d D + X_{1}+ X_{2} + \epsilon$, where $\epsilon \sim N(0,1)$ and $\beta_d = 1$.

{\bf Setting 9A}: $g_0(X^{\prime}\theta_0) = \sin(X^{\prime}\theta_0)$, where $\theta_0 = (\theta_{0,1}, \theta_{0,2})^{\prime} = (0.8,-0.6)^{\prime}$ and the vector of covariates $X=(X_{1}, X_{2})^{\prime}$ are generated from standard normal distributions $N(0,1)$. The propensity score function is defined as $D = \Lambda \big( g_0(X'\theta_{0}) \big)$ and the outcome model is generated as $Y=\beta_d D + X_{1}+ X_{2} + \epsilon$, where $\epsilon$ is the chi-square with 1 degree of freedom minus its median and $\beta_d = 1$.

{\bf Setting 9B}: $g_0(X^{\prime}\theta_0) = 0.5\{(X^{\prime}\theta_0)^3 - (X^{\prime}\theta_0)\}$, where $\theta_0 = (\theta_{0,1}, \theta_{0,2})^{\prime} = (0.8,-0.6)^{\prime}$ and the vector of covariates $X=(X_{1}, X_{2})^{\prime}$ are generated from the standard normal distribution $N(0,1)$. The propensity score function is defined as $D = \Lambda \big( g_0(X'\theta_{0}) \big)$ and the outcome model is generated as $Y=\beta_d D + X_{1}+ X_{2} + \epsilon$, where $\epsilon$ is Cauchy distribution and $\beta_d = 1$.

{\bf Setting 10A}: $g_0(X^{\prime}\theta_0) = 2X_{1}X_{2}$, where the vector of covariates $X=(X_{1}, X_{2})^{\prime}$ are generated from standard normal distributions $N(0,1)$. The propensity score function is defined as $D = \Lambda \big( g_0(X'\theta_{0}) \big)$ and the outcome model is generated as $Y=\beta_d D + X_{1}+ X_{2} + \epsilon$, where $\epsilon \sim N(0,1)$ and $\beta_d = 1$.

{\bf Setting 10B}: $D=1\{2X_{1}X_{2}+\varepsilon>0\}$, where $\varepsilon \sim N(0,1)$ is independent of $X$ and the vector of covariates $X=(X_{1}, X_{2})^{\prime}$ are generated from standard normal distributions $N(0,1)$. The outcome model is generated as $Y=\beta_d D + X_{1}+ X_{2} + \epsilon$, where $\epsilon \sim N(0,1)$ and $\beta_d = 1$.

The DGPs in Settings 7A and 7B violate $\|\theta_{0}\|=1$ in (ref), the DGPs in Settings 8A and 8B have probit propensity score function and are supposed to be handled by our semiparametric single-index models, DGPs in Settings 9A and 9B have heavy-tail error distribution of the outcome equations, and DGPs in Settings 10A and 10B impose additional misspecification of propensity score function with respect to the semiparametric single-index model.

{

}

We summarize the simulation results in (ref). Our estimator performs well under the relaxation of $\|\theta_{0}\|=1$ and misspecification of propensity score and outcome function. All estimators have relatively small biases, standard deviations, and RMSEs for the three sample sizes we studied, and all three terms generally decrease as the sample size increases. Our simulation results show that our estimators are valid for certain degree of misspecification in propensity score and outcome equation.

{\bf Appendix D.2: Comparison with Other Estimators in Propensity Score Approximation}

We employ a single-index structure as an approximation for the unknown propensity score function. In this section, we investigate its accuracy and compare the finite sample performance of our proposed method with recent approaches to propensity score estimation based on covariates balancing Anna2021JAE,imai2014covariate,Zubizarreta and high-dimensional selection belloni2017a,Chernozhukov2018,SunTan2021. To challenge these assumptions, we employ four data generating processes (DGPs) that feature unbounded link functions and unbounded support for propensity score functions.

{\bf Setting 9A}: $g_0(X^{\prime}\theta_0) = \sin(X^{\prime}\theta_0)$, where $\theta_0 = (\theta_{0,1}, \theta_{0,2})^{\prime} = (0.8,-0.6)^{\prime}$ and the vector of covariates $X=(X_{1}, X_{2})^{\prime}$ are generated from standard normal distributions $N(0,1)$. The propensity score function is defined as $D = \Lambda \big( g_0(X'\theta_{0}) \big)$.

{\bf Setting 9B}: $g_0(X^{\prime}\theta_0) = 0.5\{(X^{\prime}\theta_0)^3 - (X^{\prime}\theta_0)\}$, where $\theta_0 = (\theta_{0,1}, \theta_{0,2})^{\prime} = (0.8,-0.6)^{\prime}$ and the vector of covariates $X=(X_{1}, X_{2})^{\prime}$ are generated from the standard normal distribution $N(0,1)$. The propensity score function is defined as $D = \Lambda \big( g_0(X'\theta_{0}) \big)$.

{\bf Setting 10A}: $g_0(X^{\prime}\theta_0) = 2X_{1}X_{2}$, where the vector of covariates $X=(X_{1}, X_{2})^{\prime}$ are generated from standard normal distributions $N(0,1)$. The propensity score function is defined as $D = \Lambda \big( g_0(X'\theta_{0}) \big)$.

{\bf Setting 10B}: $D=1\{2X_{1}X_{2}+\varepsilon>0\}$, where $\varepsilon \sim N(0,1)$ is independent of $X$ and the vector of covariates $X=(X_{1}, X_{2})^{\prime}$ are generated from standard normal distributions $N(0,1)$.

(ref) displays the finite sample performance of the three approaches, where we denote our estimator as SP, the covariates balancing estimator as CB, and the high-dimensional selection estimator as HD. All three estimators exhibit near-zero bias in all cases. Our estimator also works well when the DGPs are non-single index and outperforms the covariates balancing Anna2021JAE,imai2014covariate,Zubizarreta and high-dimensional selection belloni2017a,Chernozhukov2018,SunTan2021 in terms of standard deviation and root mean squared error (RMSE) for all considered sample sizes. Overall, our method is superior to the existing methods when dealing with DGPs that are non-single index.

{\bf Appendix D.3: Comparison with Other Estimators in Treatment Effects Performance}

In this section, we compare the finite sample performance of our approach with two alternative approaches in Liu2018 and sun2021estimation. Both of these papers assume boundedness of the link function and the function support, which imposes additional constraints on their asymptotic theory and numerical performance. We use the following four DGPs with unbounded link function and unbounded support for propensity score functions.

{\bf Setting 11A}: $g_0(X^{\prime}\theta_0) = 10\exp(X^{\prime}\theta_0)$, where $\theta_0 = (\theta_{0,1}, \theta_{0,2})^{\prime} = (0.8,-0.6)^{\prime}$ and the vector of covariates $X=(X_{1}, X_{2})^{\prime}$ are generated from independent standard normal distributions $N(0,1)$. The propensity score function is defined as $D = \Lambda \big( g_0(X'\theta_{0}) \big)$ and the outcome model is generated as $Y=\beta_d D + X_{1}+ X_{2} + \epsilon$, where $\epsilon \sim N(0,1)$ and $\beta_d = 1$.

{\bf Setting 11B}: $g_0(X^{\prime}\theta_0) = 10 \{(X^{\prime}\theta_0)^5 - (X^{\prime}\theta_0)^3\} + 10\exp(X^{\prime}\theta_0)$, where $\theta_0 = (\theta_{0,1}, \theta_{0,2})^{\prime} = (0.8,-0.6)^{\prime}$ and the vector of covariates $X=(X_{1}, X_{2})^{\prime}$ are generated from independent standard normal distributions $N(0,1)$. The propensity score function is defined as $D = \Lambda \big( g_0(X'\theta_{0}) \big)$ and the outcome model is generated as $Y=\beta_d D + X_{1}+ X_{2} + \epsilon$, where $\epsilon \sim N(0,1)$ and $\beta_d = 1$.

{\bf Setting 12A}: $g_0(X^{\prime}\theta_0) = 10\exp(X^{\prime}\theta_0)$, where $\theta_0 = (\theta_{0,1}, \theta_{0,2})^{\prime} = (0.8,-0.6)^{\prime}$ and the vector of covariates $X=(X_{1}, X_{2})^{\prime}$ are generated from independent Cauchy distribution. The propensity score function is defined as $D = \Lambda \big( g_0(X'\theta_{0}) \big)$ and the outcome model is generated as $Y=\beta_d D + X_{1}+ X_{2} + \epsilon$, where $\epsilon \sim N(0,1)$ and $\beta_d = 1$.

{\bf Setting 12B}: $g_0(X^{\prime}\theta_0) = 10 \{(X^{\prime}\theta_0)^5 - (X^{\prime}\theta_0)^3\} + 10\exp(X^{\prime}\theta_0)$, where $\theta_0 = (\theta_{0,1}, \theta_{0,2})^{\prime} = (0.8,-0.6)^{\prime}$ and the vector of covariates $X=(X_{1}, X_{2})^{\prime}$ are generated from independent Cauchy distribution. The propensity score function is defined as $D = \Lambda \big( g_0(X'\theta_{0}) \big)$ and the outcome model is generated as $Y=\beta_d D + X_{1}+ X_{2} + \epsilon$, where $\epsilon \sim N(0,1)$ and $\beta_d = 1$.

Similarly, (ref) presents the finite sample performance of the three approaches under the same conditions, where we denote our estimator as SP, the locally efficient estimator as Efficient Liu2018, and the nonparametric estimator using cross-validation least squares to select bandwidth as NP sun2021estimation. Our estimator consistently exhibits a smaller bias compared to the alternative methods when the boundedness assumptions are not met. Furthermore, our estimator demonstrates a lower standard deviation and RMSE than those of Liu2018 and sun2021estimation across all sample sizes. In conclusion, our proposed method proves to be more advantageous than the existing methods for analyzing data with unbounded link functions and unbounded support.

{\bf Appendix D.4: Additional Simulation with Various Choices of Truncation Parameters}

We show simulation results for the finite sample performance of the estimation approach with various choices of truncation parameter $k$. (ref) show the simulation results of the propensity score estimation with the choices of truncation parameter $k$ ranging from $k = 2$ to $k = 6$. (ref) show the simulation results of the average therapeutic effect with the choices of truncation parameter $k$ ranging from $k = 2$ to $k = 6$. Overall, the estimators are valid against different choices of truncation parameters.

{

}

\setcounter{equation}{0} \setcounter{table}{0} \setcounter{figure}{0}

{\bf Appendix E: Additional Empirical Results}

We apply our semiparametric methods to study the effect of job training on future earnings. Our dataset is based on the classic National Supported Work (NSW) demonstration dataset, as analyzed by lalonde1986evaluating and reconstructed by dehejia1999a. The NSW experiment took place from March 1975 until June 1977 and randomly assigned participants to the treatment group who received a guaranteed job for 9 to 18 months and frequent counselor meetings or control groups who were left in the labor market by themselves.

The outcome variable $Y$ is the participant's earnings in 1978. The binary treatment variable $D$ is the job training status ($D = 1$ indicates participant is assigned training, $D = 0$ indicates participant is not assigned training). The covariates $X$ include the education, ethnicity, age, and employment variables before treatment, including earnings in 1974 and 1975. Following the imbens2015a, we augment this dataset with observations from the Current Population Survey (CPS). We focus on estimating the average treatment effect of the program for participants using the comparison group from the CPS. Due to data availability for the randomly assigned control group, we can assess whether the non-experimental estimators are accurate. We provide the summary statistics for these two datasets in (ref) and see much larger differences between the two groups, suggesting that it is important to carefully adjust for these differences in estimating the average treatment effects.

{

}

As the dataset has 7,657 observations in the control group with propensity scores less than 1.00e-05, the propensity score-related estimators may not be readily applied in this case. We compare the results of average treatment effect estimation based on bias-corrected matching (MBC) estimator abadie2006a,abadie2011a and regression adjustment (RA) estimator lane1982analysis. We summarize the estimation results for the average treatment effect of job training on earnings in 1978 in (ref) along with the mean and standard deviation. We first show the average difference between treatment group and experiment control group is 1794.34 and statistically significant. The RA estimator suggests that the job training increases the further earning by 3680.60 but is statistically insignificant. The MBC estimator shows that the average effect of 6293.11 is statistically significant but quite different from the experimental effect of 1794.34. Our ATE estimator is 4778.24 while variance-weighted ATE estimator is 2914.72. Both estimators are statistical significant, and the estimate of variance-weighted ATE is closest to the difference between treatment and experiment control group among all estimators we considered, suggesting that the our estimators may provide relatively accurate estimation and inference for the dataset with limited overlap.

{

}

We plot the estimated link function $\hat{g}$ of propensity score using the wild bootstrap simulation with 500 bootstrap repetitions in (ref). The estimated link function $\hat{g}$ is shown in a solid line, the corresponding 95% confidence interval is drawn in dashed lines. We calculate the optimal truncation parameter $k$ based on the minimization of the prediction mean of squares. In this data example, the optimal truncation parameter $k = 2$. We estimate the approximated second-order polynomial fitted function using the ordinary least squares estimation approach as a reference as following (with standard errors in bracket):

where $\omega_{i}$ is the $x_{i}'\hat{\theta}$ for $1\leq i \leq N$.

For reference, we also plot the approximated second-order polynomial fitted function using the ordinary least squares estimation approach. As shown, the estimated link functions $\hat{g}$ in the data sample of job training program are quite far away from the identity function, meaning that the logistic parametric assumption for propensity score may not be valid. One reason is that in the job training program data sample, there are 185 observations in treated group, 15,992 observations in control group, and 7,657 observations in the control group with propensity scores less than 1.00e-05. Therefore, most of the observations have the zero value for propensity score, making the regular parametric assumptions not feasible in this particular data sample. The estimated link function is very close to the approximated second-order polynomial fitted function, showing our model and method work well in this data sample.

}