# Estimation of a score-explained non-randomized treatment effect in fixed and
high dimensions

Debarghya Mukherjee, Moulinath Banerjee and Ya’acov Ritov

Department of Statistics, University of Michigan

###### Abstract

Non-randomized treatment effect models are widely used for the assessment of
treatment effects in various fields and in particular social science
disciplines like political science, psychometry, psychology. More
specifically, these are situations where treatment is assigned to an
individual based on some of their characteristics (e.g. scholarship is
allocated based on merit or antihypertensive treatments are allocated based on
blood pressure level) instead of being allocated randomly, as is the case, for
example, in randomized clinical trials. Popular methods that have been largely
employed till date for estimation of such treatment effects suffer from slow
rates of convergence (i.e. slower than $\sqrt{n}$). In this paper, we present
a new model coined SCENTS: Score Explained Non-Randomized Treatment Systems,
and a corresponding method that allows estimation of the treatment effect at
$\sqrt{n}$ rate in the presence of fairly general forms of confoundedness,
when the ‘score’ variable on whose basis treatment is assigned can be
explained via certain feature measurements of the individuals under study. We
show that our estimator is asymptotically normal in general and semi-
parametrically efficient under normal errors. We further extend our analysis
to high dimensional covariates and propose a $\sqrt{n}$ consistent and
asymptotically normal estimator based on a de-biasing procedure. Our analysis
for the high dimensional incarnation can be readily extended to analyze
partial linear models in the presence of noisy variables corresponding to the
non-linear part of the model, where the noise can be correlated with the
variables corresponding to the linear part. We analyze two real datasets via
our method and compare our results with those obtained by using previous
approaches. We conclude this paper with a discussion on some possible
extensions of our approach.

##  1 Introduction and background

Estimation of treatment effect under non-random treatment allocation has been
extensively studied in the statistics, biomedicine and econometrics
literatures. As an introduction to the idea, imagine that a scholarship
granting agency tests a group of high school students and assigns scholarship
to those whose scores are above some pre-determined cutoff (w.l.o.g. 0, after
centering). Of interest is to determine whether the scholarship has any
tangible outcome on future academic performance. Letting $Y_{i}$ be the score
of student $i$ in a subsequent semester, we can write down a model of the
type:

$Y_{i}=\alpha_{0}\mathds{1}_{Q_{i}\geq 0}+X_{i}^{\top}\beta_{0}+\nu_{i}\,$
(1.1)

where $X_{i}$ is a covariate vector including demographic information on the
students, $Q_{i}$ is the centered score on the test, and $\nu_{i}$ is a
residual term. The parameter $\alpha_{0}$ represents the effect of the
treatment: scholarship. If the hypothesis $\alpha_{0}=0$ is rejected by a
statistical test, one concludes that the scholarship has a significant impact
on subsequent academic score. A simple multiple regression cannot be performed
to estimate the treatment $\alpha_{0}$ owing to the possible dependence
between $Q$ and $\nu$. Similar examples are also prevalent in biomedicine,
especially to assess the efficacy (or possible side effects) of existing drugs
(for example, to quantify the side effects of the statin prescribed if some
threshold related to cholesterol level is exceeded or anti-depressant
prescribed when the Hamilton’s score is higher than some threshold).

It is instructive to take a brief moment to compare and contrast this setup
with traditional randomized clinical trials. In both, we estimate the effect
of a treatment which applies to a subset of the participants. The main
difference is that the allocation of treatment in the latter case is
independent of the covariates and error terms in the model, while in the
former case, the applied treatment is non-trivially correlated with the
covariate and error terms and is therefore endogenous. In our example, the
score $Q_{i}$ on the assessment test, based on which the treatment
(scholarship when $Q_{i}>0$) is applied, is not only correlated with the
indicator of treatment and the background covariates $X_{i}$, but _also_ the
unmeasured sources of variation (say, native abilities of the individual not
captured by the observed covariates) contained in the residual $\nu_{i}$. More
meritorious students are more likely to get the scholarship and are also more
likely to have higher values of $Y$, above what is predicted by the
scholarship effect and their demographics! This typically translates to what
is called the ‘endogeneity assumption’: $E(\nu_{i}|Q_{i})\neq 0$, and hence
the model in equation (1.1) cannot be estimated using simple linear regression
of $Y$ on $(1(Q>0),X)$ unlike the randomized trial framework. This model was
initially studied under the name of regression discontinuity design in ([44],
[12]), who analyzed data from a national merit competition (see [19] and [43])
and have since found varied applications, for example [41], [15] in education,
[26] in health related research, [6] in social research, [16] in epidemiology,
to name a few.

Most of the analysis in such jump-effect models, till date, is based off of a
local analysis, since observations from a small neighborhood of the cut-off
$Q$ are considered as almost free from endogeneity and only those are used to
get an approximately unbiased estimate of the treatment effect. This is the
main reason why these methods only provide an estimation of local average of
the treatment effect around the cut-off. However, in many practical scenarios,
the interest is _not solely_ on the jump at the discontinuity but also on the
_global effect_ of the treatment, which also requires assessing treatment
effectivity for individuals that lie substantially away from the
discontinuity. If we consider for example a prestigious program, most of the
investment is in students who are clearly above the threshold, and the social
implications of the program are far from being restricted to those near it.
Thus, measuring the effect of the jump should ideally keep the large
implication in mind, whilst in the standard regression discontinuous design
analysis, it is only the local phenomenon around the jump that is the object
of focus. The aim of this manuscript is to propose and analyze a model which
captures this overall phenomenon.

To highlight our contribution and contrast it with the extant literature on
the estimation of treatment effect in presence of endogeneity, let’s first
consider one of the standard approaches, as exemplified by the model below:

$Y_{i}=(\alpha_{0}+b_{0}Q_{i})\mathds{1}_{Q_{i}<\tau_{0}}+(\alpha_{1}+b_{1}Q_{i})\mathds{1}_{Q_{i}\geq\tau_{0}}+X_{i}^{\top}\beta_{0}+\nu_{i}\,,$
(1.2)

where $Q_{i}$ is the score variable which determines treatment via the known
threshold $\tau_{0}$, and $X_{i}$’s are background covariates. The parameter
of interest is $\alpha_{1}-\alpha_{0}$ which encodes the effect of the
treatment. Generally, only weak assumptions are made on the conditional
expectation of $E(Y\mid Q,X)$ to encode possible endogeneity. As the main
takeaway from the model is that both the intercept and the slope of $Q$ change
at the cut-off, the traditional approach for the estimation of the treatment
effect is as follows: Select a possibly (data-driven) bandwidth $h_{n}$ and
look at all the observations $Y_{i}^{\prime}s$ for which $|Q_{i}-\tau_{0}|\leq
h_{n}$. Now run a weighted local polynomial regression on these observations
to estimate $\alpha_{1}-\alpha_{0}$. As this approach effectively uses
$O_{p}(nh_{n})$ observations for estimation purposes, the rate of convergence
of the estimator is _slower that $\sqrt{n}$_ (typically $\sqrt{nh_{n}}$, where
the bandwidth $h_{n}\to 0$ is chosen by standard bias-variance trade-off),
see, for example [10]. However, the precision with which the treatment-effect
is estimated can be improved in situations where it is possible to exploit the
relationship between the score $Q$ and the measured covariates $X$ on the
individuals of interest. In such cases, we argue below that it becomes
possible to use _all our samples_ to estimate the treatment effect at
$\sqrt{n}$ rate.

Our main contribution:  We assume that the linear equation (1.1) can be
augmented by a second equation that explains the score $Q$ through the
background covariates $X$ to obtain:

$\displaystyle Y_{i}=\alpha_{0}\mathds{1}_{Q_{i}\geq
0}+X_{i}^{\top}\beta_{0}+\nu_{i}$ $\displaystyle Q_{i}=g(X_{i})+\eta_{i}\,.$
(1.3)

The effect of unmeasured covariates that can affect both $Q$ and $Y$ is
encoded by the mean 0 error vector $(\nu_{i},\eta_{i})$ on which we place no
parametric assumptions, they may be dependent and $b(\eta)\equiv
E(\nu|\eta)\neq 0.$ In fact, one can (and indeed this may be necessary under
certain circumstances) generalize the second equation and write:

$\displaystyle Y_{i}=\alpha_{0}\mathds{1}_{Q_{i}\geq
0}+X_{i}^{\top}\beta_{0}+\nu_{i}$ $\displaystyle Q_{i}=g(Z_{i})+\eta_{i}\,.$
(1.4)

where $X$ and $Z$ may be identical or completely disjoint or may share several
covariates. The intuition behind this generalization is that it is quite
possible that some extra covariates are available while measuring $Y$ but not
while measuring $Q$ or vice-versa. The updated equation (1) takes care of
that.

At first glance, the augmentation of the second equation in (1) may suggest
that $Z$ is being used as an instrumental variable. However, this is not
necessarily the case in our model since the _exclusion restriction_ [27] that
is critical in the instrumental variables approach, namely that the instrument
should influence $Y$ only through $Q$ and not directly, is no longer satisfied
(especially when $X=Z$, which is typically the case for many real
applications). Consequently, the two-stage least squares procedure [1]
typically employed in the instrumental variable literature is _invalid_ in our
problem. We emphasize that, in contrast to instrumental variable regression,
we _do not need_ to deploy a completely new set of covariates to explain $Q$,
rather using at least part of the already available background information $X$
that influences $Y$ to also explain the score $Q$. This is typically the case
in many real life examples, e.g. in the scholarship example, the background
information, that explains the future academic performance of a student (Y)
also influences their scholarship test score $Q$.

In this paper, we analyze a simpler version of the above model (equation (1))
for technical simplicity:

$\displaystyle Y_{i}=\alpha_{0}\mathds{1}_{Q_{i}\geq
0}+X_{i}^{\top}\beta_{0}+\nu_{i}$ $\displaystyle
Q_{i}=Z_{i}^{\top}\gamma_{0}+\eta_{i}\,.$ (1.5)

where we assume the error $(\nu,\eta)\perp\perp(X,Z)$, $(X,Z)$ can be
arbitrarily related, and $\nu$ is _correlated_ with $\eta$. However our
methods can be extended to the more general model (1) as discussed in
Subsection 6.1. 111Note that if $Z\perp\perp X$, then one may use standard
techniques from instrumental variable regression. However, this fails to hold
in most our the intended applications.

We next discuss how this augmentation helps us obtain a $\sqrt{n}$-consistent
estimator. _If we observe $\eta$_, writing
$\nu_{i}=b(\eta_{i})+{\epsilon}_{i}$, our model would reduce to a simple
partial linear model (see e.g. [48], [7], [30] and references therein) with
$(\alpha_{0},\beta_{0})$ being the parametric component and the unknown mean
function $b$ being the nonparametric component:

$Y_{i}=\alpha_{0}\mathds{1}_{Q_{i}\geq
0}+X_{i}^{\top}\beta_{0}+b(\eta_{i})+{\epsilon}_{i}\,.$ (1.6)

Hence, we could acquire a $\sqrt{n}$ consistent estimator of $\alpha_{0}$
following the standard analysis of the partial linear model _provided
$\mathbb{E}({\sf var}(\mathds{1}(Q>0)\mid\eta))>0)$_ (for more details about
semiparametric efficient estimation in partial linear models, see [35] or
[39]). Indeed, without this assumption it is easy to see that there is an
identifiability problem and the effect of $\alpha_{0}$ cannot be separated
from the effect of $b$ in the above model.

As _we don’t observe $\eta$_, we cannot use the standard partial linear model
analysis proposed in [35] or [39] to estimate $\alpha_{0}$. We can
nevertheless approximate it by (1), using the residuals obtained via
regressing $Q$ on $Z$, i.e. plug-in the estimates of $\eta_{i}$’s in equation
(1.6) and treat the resulting equation as an _approximate partial linear
model_. Indeed, this idea lies at the heart of our method. However, a naive
replacement of $\eta$ by $\hat{\eta}$ (obtained by regressing $Q$ on $Z$) is
not sufficient, as the approximation error of $\eta$ is
$\hat{\eta}-\eta=O_{p}(n^{-1/2})$, whilst we need the approximation error to
be of the order $o_{p}(n^{-1/2})$ for $\sqrt{n}$-consistent estimation of
$\alpha_{0}$. Hence, some more fine-tuning is needed to remove the bias. We
elaborate our method of estimation in Section 2.

The assumption $\gamma_{0}\neq 0$ is critical to our analysis, because if not,
$Q_{i}=\eta_{i}$ and equation (1.6) becomes:

$Y_{i}=\alpha_{0}\mathds{1}_{Q_{i}\geq
0}+X_{i}^{\top}\beta_{0}+b(Q_{i})+{\epsilon}_{i}\,.$

It is now no longer possible to estimate $\alpha_{0}$ at $\sqrt{n}$ rate as it
is hard to separate the first and third term of the RHS. Thus, the augmenting
equation $Q=Z^{\top}\gamma_{0}+\eta$ with $\gamma_{0}\neq 0$ is what prevents
the variable corresponding to the parametric component of interest from
becoming a measurable function of the variable corresponding to the non-
parametric component of the model, and enables converting the estimation
problem to an _approximate_ partial linear model. As noted earlier, for
$\alpha_{0}$ to be estimated as $\sqrt{n}$ rate, it is not necessary for $Q$
to have a linear relation with $Z$, in fact any non-linear parametric relation
like $Q=g(Z)+\eta$ for some known link function $g$ will work as long as
$\mathbb{E}\left({\sf var}(\mathds{1}(Q>0)\mid\eta)\right)>0$. Also the
assumption $(\nu,\eta)\perp\perp(X,Z)$ is not necessary, all we need is that
$\mathbb{E}[\eta\mid Z]=0$ and $\mathbb{E}[{\epsilon}\mid X,Z,\eta]=0$ for
$\sqrt{n}$-consistent estimation of $\alpha_{0}$.

There is also a degree of similarity between our model and _triangular
simultaneous equation models_ , which are well-studied in the economics
literature. Interested readers may take a look at [31], or [32] and references
therein for more details. However, there are two key differences: (1) The
triangular simultaneous equation models generally assume a smooth link
function between $Y_{i}$ and $(X_{i},Z_{i},Q_{i})$, whilst our model presents
an inherent discontinuity. Hence our work _cannot be derived_ from analyses of
the triangular simultaneous equation model. (2) The direct influence of X on Y
would ruin the identification in these non-separable models.

Estimation of the treatment effect via partial linear model in the context of
regression discontinuity design has also been mentioned in [33] through the
lens of the following model:

$Y_{i}=\alpha_{0}\mathds{1}_{Q_{i}>0}+b(Q_{i})+{\epsilon}_{i},\ \ \
\mathbb{E}[{\epsilon}\mid Q,\mathds{1}_{Q>0}]=0\,.$

The author argued that one _cannot estimate_ $\alpha_{0}$ at $\sqrt{n}$ rate
as $S_{i}=\mathds{1}_{Q_{i}>0}$ can be completely explained by $Q_{i}$, i.e.
$\mathbb{E}[{\sf var}(S\mid Q)]=\mathbb{E}\left[(S-\mathbb{E}[S\mid
Q])^{2}\right]=0$ (see the discussion after equation (1.6)). The main
difference between our model and that of [33] is the usage of the background
information $X$ (or $(X,Z)$) which influences both $Q$ and $Y$. This, on one
hand, prevents us from using $X$ as simple instruments, and on the other hand
ensures that $S$ is not completely explained by $\eta$ (the random variable
corresponding non-parametric component $b$) at it also depends on $X$, which
enables us to estimate the parametric component at $\sqrt{n}$ rate.

Recently [2] presented an idea where one may exploit background covariates to
estimate the treatment effect at $\sqrt{n}$ rate under the assumption
$\mathbb{E}[Y\mid Q,X]=\mathbb{E}[Y\mid X]$, i.e. the relation between $Q$ and
$Y$ can be explained fully by $X$, which permits them to perform simple OLS
method to obtain a $\sqrt{n}$-consistent estimator of the treatment effect.
However, their assumption is not satisfied in our case: the main component of
our model is the unobserved $\eta$ that quantifies the _innate_ ability of an
individual (i.e. intelligence), which, along with $X$’s, influences both $Q$
and $Y$ (in-fact it might even be argued that $\eta$ is as important as $X$,
especially in the scholarship example like one presented in [2], as innate
abilities plays a primary role in the performance of an individual in the
scholarship test and beyond).

Extension to high dimensional setup: We have also extended our analysis to the
high dimensional model, where we assume both the dimension of $X$ and $Z$ are
much larger than the available sample size. Analysis of treatment effect in
presence of high dimensional covariates is a relatively new topic. For
example, [5] proposed a debiased approach for inference on the treatment
effect in presence of high dimensional controls, and, recently, [34] studied
heterogeneous treatment effect in the presence of high dimensional covariates
under the availability of instrumental variables; [17] proposed a modeling
strategy for personalized medicine in presence of high dimensional covariates
and [3] presented an approach to estimate the parameters of a standard
regression discontinuity design in high dimensions.

Following the previous discussion, our model in the high dimensional setup can
be viewed as a variant of the high dimensional partial linear model where we
do not observe the random variable corresponding to the non-linear part, but
instead a noisy version of it, with the noise being correlated with both the
random variables corresponding to the linear and non-linear parts. Several
efforts have been made to estimate both the linear and the non-linear part of
a high dimensional partial linear model, of which we mention a few here. [24]
used a projected spline estimator along with the $\ell_{1}$ penalty to
estimate the parametric and non-parametric parts of a high dimensional partial
linear model; [28] extended the analysis to the distributed setup; [30] used
the $\ell_{1}$ penalty for the linear part and a smoothing penalty (i.e.
penalty involving the double derivative of the non-parametric function) for
the non-linear part; [18] proposed an adaptive approach based on [20]’s
pairwise difference based method which does not require knowledge of any
function class a-priori; [50] proposed a method to de-bias the estimator of
the parameter corresponding to the linear part which facilitates inference for
any particular co-ordinate or some linear combination of the parameters
corresponding to the linear part. To the best of our knowledge, the closest
paper that studies a similar situation (i.e. partial linear model with
erroneous observations corresponding to the non-linear part) is [49]. Our
analysis is different from that of [49] for the following reasons: (i) our
basic estimation procedure (elaborated in the subsequent sections) is
different as our model is quite different from that considered in [49] and
(ii) we provide a de-biased estimate of the treatment effect $\alpha_{0}$ and
establish asymptotic normality for purposes of inference, whereas [49] only
calculate the rate of the estimation error in their model.

To summarize, our work makes the following contributions:

  1. 1.

It is able to take care of endogeneity among the errors in a general manner
and provide a $\sqrt{n}$-consistent estimator of the treatment effect. Indeed,
the main feature of our approach lies in modeling $Q$ itself in terms of
covariates up to error terms, which enables _the use of the entirety of data
available_ and _not just the observations in a small vicinity of the boundary
defined by the $Q$ threshold_, on which existing approaches are typically
based.

  2. 2.

Our estimate achieves semiparametric efficiency under an appropriate submodel.

  3. 3.

Our method does not depend on tuning parameters, in the sense that the use of
tuning parameters to estimate the treatment effect is _secondary_. As will be
seen in Section 2, we do require tuning parameter specifications for non-
parametric estimation of $b(\eta)$, but as long as those parameters satisfy
some minimal conditions, our estimate of $\alpha_{0}$ – in terms of both rate
of convergence and asymptotic distribution– does not depend on it.

  4. 4.

Our analysis of SCENT in presence of high dimensional covariates appears to be
the first systematic attempt to deal with de-biasing in the high dimensional
partial linear model where we observe a noisy version of the random variable
corresponding to the non-linear part, and the noise is correlated with both
the random variables corresponding to the linear and the non-linear parts.

Organization of the paper:  In Section 2 we describe the estimation procedure.
Section 3 provides the theoretical results along with a brief outline of the
proof of asymptotic normality of our estimator. In Section 5, we provide
analyses of two real data examples using our method, as well as comparisons to
previous methods. In Section 6 we present some possible future research
directions based on this work. In particular, we discuss the scenario where
the treatment effect, which is assumed to be constant for the model studied in
this paper, can depend on the innate ability $\eta$: in other words, what
happens when $\alpha_{0}$ is replaced by $\alpha(\eta)$, a generalization that
may be warranted in certain applications. We point out that
$\sqrt{n}$-consistent estimates of the integrated treatment effect can be
obtained in this situation as well. Rigorous proofs of the main results are
established in Section B \- Section D of the Supplementary document. In
Section E of the Supplementary document, we present proofs of auxiliary
results that are required to prove the main results. In Section F of the
Supplementary document, we provide a few details on the spline estimation
techniques used in our analysis and in Section G of the Supplementary
document, an algorithm integrating the main steps of the proposed methodology
for the ease of implementation.

Notation: For any matrix $A$, we denote by $A_{i,*}$, the $i^{th}$ row of $A$
and by $A_{*,j}$ the $j^{th}$ column of $A$. Both $a_{n}\lesssim_{P}b_{n}$ and
$a_{n}=O_{p}(b_{n})$ bear the same meaning, i.e. $a_{n}/b_{n}$ is a tight
(random) sequence. Also, for two non-negative sequences $\\{a_{n}\\}$ and
$\\{b_{n}\\}$, we denote by $a_{n}\gg b_{n}$ (respectively $a_{n}\ll b_{n}$),
the conditions that $\liminf_{n}a_{n}/b_{n}\to\infty$ (respectively
$\limsup_{n}a_{n}/b_{n}\to 0$). For any random variable (or object) $X$, we
denote by $\mathcal{F}(X)$, the sigma-field generated by $X$. For two random
variables $X$ and $Y$, $X\perp\perp Y$ indicates that $X$ and $Y$ are
independent, $X\perp Y$ denotes that they are uncorrelated.

##  2 Estimation procedure for $\alpha_{0}$

As mentioned in Section 1, our model for estimating the treatment effect can
be written as:

$\displaystyle Y_{i}$
$\displaystyle=\alpha_{0}\mathds{1}_{Z_{i}^{\top}\gamma_{0}+\eta_{i}\geq
0}+X_{i}^{\top}\beta_{0}+\nu_{i}\,.$

where $\nu_{i}$ and $\eta_{i}$ are correlated. Defining
$b(\eta_{i})=\mathbb{E}(\nu_{i}\mid\eta_{i})$, write
$\nu_{i}=b(\eta_{i})+\epsilon_{i}$ where $\epsilon_{i}\perp\eta_{i}$. Using
this we can rewrite our model as:

$\displaystyle Y_{i}$
$\displaystyle=\alpha_{0}\mathds{1}_{Z_{i}^{\top}\gamma_{0}+\eta_{i}\geq
0}+X_{i}^{\top}\beta_{0}+b(\eta_{i})+\epsilon_{i}\,.$ (2.1)

We first divide the whole data in three (almost) equal parts, say
$\mathcal{D}_{1},\mathcal{D}_{2},\mathcal{D}_{3}$. Henceforth for simplicity
we assume each $\mathcal{D}_{i}$ has $n/3$ observations. Denote the dimension
of $X$ and $Z$ by $p_{1}$ and $p_{2}$ respectively. The first two data sets
are used to obtain (consistent) estimates of several nuisance parameters which
are then plugged into the equations corresponding to the third data set from
which an estimator of the treatment effect is constructed. The data splitting
technique makes the theoretical analysis more tractable as one can use
independence among the three subsamples to our benefit. Furthermore, by
rotating the samples (to be elaborated below), we obtain three asymptotically
independent and identically distributed estimates which are then averaged to
produce a final estimate that takes advantage of the full sample size. While
we believe that the estimator obtained without data-splitting achieves the
same asymptotic variance, a point that is corroborated by simulation studies
(not reported in the manuscript), an analysis of this estimator would be
incredibly tedious with minimal further insight.

We, first, estimate $\gamma_{0}$ from $\mathcal{D}_{1}$ via standard least
squares regression of $Z$ on $Q$:

$\hat{\gamma}_{n}=(Z^{\top}Z)^{-1}Z^{\top}Q\,.$

Using $\hat{\gamma}_{n}$, we expand our first model equation as :

$\displaystyle Y_{i}$
$\displaystyle=\alpha_{0}S_{i}+X_{i}^{\top}\beta_{0}+b(\hat{\eta}_{i})+(\eta_{i}-\hat{\eta}_{i})b^{\prime}(\hat{\eta}_{i})+R_{1,i}+\epsilon_{i}$
$\displaystyle=\alpha_{0}S_{i}+X_{i}^{\top}\beta_{0}+b(\hat{\eta}_{i})+b^{\prime}(\hat{\eta}_{i})Z_{i}^{T}(\hat{\gamma}_{n}-\gamma_{0})+R_{1,i}+\epsilon_{i}\,.$
(2.2)

where $S_{i}=\mathds{1}_{Q_{i}>0}$,
$\hat{\eta}_{i}=Q_{i}-Z_{i}^{\top}\hat{\gamma}_{n}$ and
$R_{1,i}=(\hat{\eta}_{i}-\eta_{i})^{2}b^{\prime\prime}(\tilde{\eta}_{i})/2$ is
the residual obtained from a two-step Taylor expansion with $\tilde{\eta}_{i}$
lying between $\eta_{i}$ and $\hat{\eta}_{i}$. It should be pointed out that
if $\eta_{i}$ were known, our model (equation (2.1)) would reduce to a simple
partial linear model and estimation of $\alpha_{0}$ would become straight-
forward. As we _don’t_ observe $\eta_{i}$, but rather use $\hat{\eta}_{i}$ as
its proxy, the corresponding approximation error needs careful handling, as we
need to show that the estimator does not inherit any resulting bias. Indeed,
this is one of the core technical challenges of this paper, and explained in
the subsequent development.

Note that the function $b^{\prime}$ in equation (H.2) is unknown. Since
$\hat{\gamma}_{n}=\gamma_{0}+O_{p}(n^{-1/2})$ (which, in turn, implies
$\hat{\eta}_{i}=\eta_{i}+O_{p}(n^{-1/2})$), as long as we have a consistent
estimate of $b^{\prime}$, say $\hat{b}^{{}^{\prime}}$, the approximation error
$(b^{\prime}(\hat{\eta}_{i})-\hat{b}^{\prime}(\hat{\eta}_{i}))Z_{i}(\hat{\gamma}_{n}-\gamma_{0})$
is $o_{p}(n^{-1/2})$ and therefore asymptotically negligible.

We now elaborate how we use B-spline basis to estimate $b^{\prime}$ from
$\mathcal{D}_{2}$. Equation (2.1) can be rewritten as:

$\displaystyle Y_{i}$
$\displaystyle=\alpha_{0}S_{i}+X_{i}^{\top}\beta_{0}+b(\hat{\eta}_{i})+\epsilon_{i}+R_{i}\,.$
(2.3)

where $R_{i}=b(\eta_{i})-b(\hat{\eta}_{i})$. Ignoring the remainder term
(which is shown to be asymptotically negligible under some mild smoothness
condition on $b$ to be specified later), we estimate $b^{\prime}$ via a
B-spline basis from equation (2.3). The theory of spline approximation is
mostly explored for estimating compactly supported non-parametric regression
functions. Our errors $\eta$ are, of course, assumed to be unbounded as
otherwise the problem would become artificial. However, there are certain
technical issues with estimating $b^{\prime}$ on the entire support (see
Remark 2.4), and to circumvent that we restrict ourselves to a compact (but
arbitrary) support $[-\tau,\tau]$, i.e. we consider those observations for
which $|\hat{\eta}_{i}|\leq\tau$. We then use a cubic B-spline basis
appropriately scaled to the interval of interest with equispaced knots to
estimate $b^{\prime}$ 222While we work with cubic splines, one may certainly
use higher degree polynomials. However, in practice, it has been observed that
cubic spline works really well in most of the scenarios. Notationally, we use
$K-1$ knots to divide $[-\tau,\tau]$ into $K$ intervals of length $2\tau/K$
where $K=K_{n}$ increases with $n$ at an appropriate rate (see remark 2.2
below), giving us in total $K+3$ spline basis functions. For any $x$, we use
the notation $\tilde{N}_{K}(x)\in\mathbb{R}^{(K+3)}$ to denote the vector of
scaled B-spline basis functions $\\{\tilde{N}_{K,j}\\}_{j=1}^{K+3}$ evaluated
at $x$333A brief discussion on B-spline basis and scaled B-spline basis is
presented in Section F of the Supplementary document for the ease of the
reader.. Using these basis functions we further expand equation (2.3) as
follows:

$Y_{i}=\alpha_{0}S_{i}+X_{i}^{\top}\beta_{0}+\tilde{N}_{K}(\hat{\eta}_{i})^{\top}\omega_{b,\infty,n}+\epsilon_{i}+R_{i}+T_{i}\,.$
(2.4)

for all those observations with $|\hat{\eta}_{i}|\leq\tau$, where $T_{i}$ is
the spline approximation error, and
$T_{i}=b(\hat{\eta}_{i})-N_{K}(\hat{\eta}_{i})^{\top}\omega_{b,\infty,n}$ and
$\omega_{b,\infty,n}$ is the (population) parameter defined as:

$\omega_{b,\infty,n}=\mathop{\rm
argmin}_{\omega\in\mathbb{R}^{(K+3)}}\sup_{|x|\leq\tau}\left|b(x)-\tilde{N}_{K}^{\top}(x)\omega\right|\,.$

Suppose we have $n_{2}\leq n/3$ observations in $\mathcal{D}_{2}$ with
$|\hat{\eta}_{i}|\leq\tau$. Denote by $\mathbf{Y}$ the vector of all the
corresponding $n_{2}$ responses, by $\mathbf{X}\in\mathbb{R}^{(n_{2},p_{1})}$
the covariate matrix, and by
$\tilde{\mathbf{N}}_{K}\in\mathbb{R}^{(n_{2},K+3)}$ the approximation matrix
with rows $\tilde{\mathbf{N}}_{K,i*}=\tilde{N}_{K}(\hat{\eta}_{i})$.
Regressing $Y$ on $(\mathbf{S},\mathbf{X},\tilde{\mathbf{N}}_{K})$ we estimate
$\omega_{b,\infty,n}$ (details can be found in the proof of Proposition 2.1)
and set
$\hat{b}^{\prime}(x)=\nabla\tilde{N}_{K}(x)^{\top}\hat{\omega}_{b,\infty,n}$
where $\nabla\tilde{N}_{K}(x)$ is the vector of derivates of each of the co-
ordinates of $\tilde{N}_{K}(x)$. The following theorem establishes consistency
of our estimator of $b^{\prime}$ (proof can be found in Section D of the
supplementary document):

######  Proposition 2.1.

Under Assumptions 3.1 \- 3.4 (elaborated in Section 3) we have:

$\sup_{|x|\leq\tau}\left|b^{\prime}(x)-\hat{b}^{\prime}(x)\right|=\sup_{|x|\leq\tau}\left|b^{\prime}(x)-\nabla\tilde{N}_{K}(x)^{\top}\hat{\omega}_{b,\infty,n}\right|=o_{p}(1)\,.$

######  Remark 2.2.

Henceforth, we choose $K\equiv K_{n}$ such that $n^{1/8}\ll K\ll n^{1/3}$ to
control the approximation errors of certain non-parametric functions
(including $b(\eta)$) that appear in our analysis via the B-spline basis.
However the bounds can be improved in presence of additional derivations of
the non-parametric functions involved.

The final (key) step involves $\alpha_{0}$ from $\mathcal{D}_{3}$. Suppose
there are $n_{3}\leq n/3$ observations in $\mathcal{D}_{3}$ with
$|\hat{\eta}_{i}|\leq\tau$. Replacing $b^{\prime}$ by $\hat{b}^{\prime}$
obtained from $\mathcal{D}_{2}$ in equation (H.2) we obtain:

$\displaystyle Y_{i}$
$\displaystyle=\alpha_{0}S_{i}+X_{i}^{\top}\beta_{0}+b(\hat{\eta}_{i})+\hat{b}^{\prime}(\hat{\eta}_{i})Z^{\top}_{i}(\hat{\gamma}_{n}-\gamma_{0})+R_{1,i}+R_{2,i}+\epsilon_{i}$
$\displaystyle\triangleq\alpha_{0}S_{i}+X_{i}^{\top}\beta_{0}+b(\hat{\eta}_{i})+\tilde{Z}^{\top}_{i}(\hat{\gamma}_{n}-\gamma_{0})+R_{1,i}+R_{2,i}+\epsilon_{i}\,,$
(2.5) $\displaystyle Q_{i}-Z_{i}^{\top}\hat{\gamma}_{n}$
$\displaystyle=-Z_{i}^{\top}(\hat{\gamma}_{n}-\gamma_{0})+\eta_{i}\,.$

where we define $\tilde{Z}_{i}=\hat{b}^{\prime}(\hat{\eta}_{i})Z_{i}$, the
residual term $R_{2,i}$ as:

$R_{2,i}=\left(b^{\prime}(\hat{\eta}_{i})-\hat{b}^{\prime}(\hat{\eta}_{i})\right)Z_{i}^{\top}\left(\hat{\gamma}_{n}-\gamma_{0}\right)=\left(b^{\prime}(\hat{\eta}_{i})-\hat{b}^{\prime}(\hat{\eta}_{i})\right)\left(\hat{\eta}_{i}-\eta_{i}\right)\,.$

and $R_{1,i}$ is same as in equation (H.2). An inspection of the first part of
Equation (2) shows that up to the remainder terms
$\left\\{\left(R_{1,i},R_{2,i}\right)\right\\}_{i=1}^{n_{3}}$, our model is a
partial linear model with parameters $(\alpha_{0},\beta_{0},b)$. These
remainder terms are asymptotically negligible as shown in the proof of our
main theorem (Theorem 3.6). We estimate $\alpha_{0}$ from equation (2) using
standard techniques for the partial linear model which, again, involve
approximating the function $b$ with the same B-spline basis as before:

$\displaystyle Y_{i}$
$\displaystyle=\alpha_{0}S_{i}+X_{i}^{\top}\beta_{0}+\tilde{Z}_{i}^{\top}(\hat{\gamma}_{n}-\gamma_{0})+b(\hat{\eta}_{i})+R_{1,i}+R_{2,i}+{\epsilon}_{i}$
$\displaystyle=\alpha_{0}S_{i}+X_{i}^{\top}\beta_{0}+\tilde{Z}_{i}^{\top}(\hat{\gamma}_{n}-\gamma_{0})+\tilde{N}_{K}(\hat{\eta}_{i})^{\top}\omega_{b,\infty,n}+R_{1,i}+R_{2,i}+R_{3,i}+{\epsilon}_{i}\,.$
(2.6)

where $R_{3,i}$ is the spline approximation error, i.e.
$R_{3,i}=b(\hat{\eta}_{i})-\tilde{N}_{K}(\hat{\eta}_{i})^{\top}\omega_{b,\infty,n}$.
Combining equation (2.6) along with the second equation of (2) we formulate
the following linear model equation:

$\displaystyle\begin{pmatrix}\mathbf{Y}\\\
\mathbf{Q}-\mathbf{Z}\hat{\gamma}_{n}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}\mathbf{S}&\mathbf{X}&\tilde{\mathbf{Z}}&\tilde{\mathbf{N}}_{K}\\\
0&0&-\mathbf{Z}&0\end{pmatrix}\begin{pmatrix}\alpha_{0}\\\ \beta_{0}\\\
\hat{\gamma}_{n}-\gamma_{0}\\\
\omega_{b,\infty,n}\end{pmatrix}+\begin{pmatrix}\mathbf{R}\\\
0\end{pmatrix}+\begin{pmatrix}\epsilon\\\ \eta\end{pmatrix}$
$\displaystyle=\begin{pmatrix}\mathbf{W}&\tilde{\mathbf{N}}_{K,a}\end{pmatrix}\begin{pmatrix}\theta_{0}\\\
\omega_{b,\infty,n}\end{pmatrix}+\begin{pmatrix}\mathbf{R}\\\
0\end{pmatrix}+\begin{pmatrix}\epsilon\\\ \eta\end{pmatrix}$ (2.7)

where $\theta_{0}=(\alpha_{0},\beta_{0},\hat{\gamma}_{n}-\gamma_{0})$ ,

$\mathbf{W}=\begin{bmatrix}\mathbf{W}_{1}\\\ \mathbf{W}_{2}\end{bmatrix}\ \
\text{ with }\ \
\mathbf{W}_{1}=\begin{pmatrix}\mathbf{S}&\mathbf{X}&\tilde{\mathbf{Z}}\end{pmatrix},\
\ \ \mathbf{W}_{2}=\begin{pmatrix}0&0&-\mathbf{Z}\end{pmatrix}\,,$

where
$\mathbf{W}_{1}\in\mathbb{R}^{(n_{3},1+p_{1}+p_{2})},\mathbf{W}_{2}\in\mathbb{R}^{(n/3,1+p_{1}+p_{2})}$
(as we are using all the observations in $\mathcal{D}_{3}$ in the second
regression equation of (2)), $\mathbf{R}$ is the vector of the sum of three
residuals $R_{1,i},R_{2,i},R_{3,i}$ mentioned in equation (2.6), and the
matrix $\tilde{\mathbf{N}}_{K,a}$ (read $\tilde{\mathbf{N}}_{K}$ appended) has
the form:

$\tilde{\mathbf{N}}_{K,a}=\begin{bmatrix}\tilde{\mathbf{N}}_{K}\\\
0\end{bmatrix}\,.$

From equation (2), we estimate $\theta_{0}$ via ordinary least squares
methods:

$\hat{\theta}=\left(\mathbf{W}^{\top}{\sf
proj}^{\perp}_{\tilde{\mathbf{N}}_{K,a}}\mathbf{W}\right)^{-1}\mathbf{W}^{\top}{\sf
proj}^{\perp}_{\tilde{\mathbf{N}}_{K,a}}\begin{pmatrix}\mathbf{Y}\\\
\mathbf{Q}-\mathbf{Z}\hat{\gamma}_{n}\end{pmatrix}$

and set $\hat{\alpha}$ as the first co-ordinate of $\hat{\theta}$, where for
any matrix $A$, we define ${\sf proj}_{A}$ as the projection matrix on the
columns space of $A$. Note that because of data-splitting, $\hat{\gamma}_{n}$,
$\hat{b}^{\prime}$ and $\mathcal{D}_{3}$ are mutually independent, which
provides a significant technical advantage in dealing with the asymptotics of
our estimator.

Finally, we apply the same methodology on permutations of the three sets of
data: i.e. we estimate $\hat{\gamma}_{n}$ from $\mathcal{D}_{2}$,
$\hat{b}^{\prime}$ from $\mathcal{D}_{3}$, $\hat{\alpha}\in\mathcal{D}_{1}$
and $\hat{\gamma}_{n}$ from $\mathcal{D}_{3}$, $\hat{b}^{\prime}$ from
$\mathcal{D}_{1}$, $\hat{\alpha}\in\mathcal{D}_{2}$. Denote by
$\hat{\alpha}_{i}$, the estimator $\hat{\alpha}$ estimated from
$\mathcal{D}_{i}$ for $1\leq i\leq 3$. Our final estimate is then
$\overline{\hat{\alpha}}=(1/3)\sum_{i}\alpha_{i}$.

######  Remark 2.3.

In our estimation procedure, we effectively estimate the conditional mean
function $b(\cdot)$ twice: once while estimating $b^{\prime}$ from
$\mathcal{D}_{2}$ and again while estimating $\alpha_{0}$ from
$\mathcal{D}_{3}$. Note that, the second re-estimation of $b$ is quite
critical (i.e. we cannot _use_ the estimator of $b$ obtained from
$\mathcal{D}_{2}$) due to presence of higher order bias (slower than
$n^{-1/2}$).

######  Remark 2.4.

As described in this section, we only use observations for which
$|\hat{\eta}_{i}|\leq\tau$, losing some efficiency in the process. One way to
circumvent this issue is to use a sequence $\\{\tau_{n}\\}$ slowly increasing
to $\infty$ and considering all those observations for which
$|\hat{\eta}|\leq\tau_{n}$. Although, this will acquire efficiency in the
limit, we need a stronger set of assumptions to make it work: for starters, we
need conditions on the decay of the density of $\eta$, we need it not to
vanish anywhere in $[-\tau_{n},\tau_{n}]$, and knowledge of the rate at which
$\min_{|x|\leq\tau_{n}}f_{\eta}(x)$ approaches $0$. We also need stronger
conditions on some conditional expectation functions (i.e. conditional
expectation of $\begin{pmatrix}S&\mathbf{X}&\mathbf{Z}\end{pmatrix}$ given
$\eta+a^{\top}Z$ for some vector $a$), e.g., bounded derivatives in both co-
ordinates over the entire space (see Lemma B.1 of the Supplementary document.)
With more technical nuances we believe our method can be extended to the
entire real line by using a growing interval, but from a purely statistical
angle, it will not bring anything insightful to the methodology that we
propose here.

##  3 Analysis for fixed dimensional covariates

In this section we presents our main theorems with broad outline of the
proofs. Details are provided in the Supplementary document. To establish the
theory, we need the following assumptions:

######  Assumption 3.1.

The errors $(\eta,\nu)$ are independent of the distribution of $(X,Z)$ and
have zero expectation.

######  Assumption 3.2.

The distribution of $(\eta,\nu)$ satisfies the following conditions:

  1. i)

The density of $\eta$, denoted by $f_{\eta}$, is continuously differentiable
and both $f_{\eta}$ and its derivatives are uniformly bounded.

  2. ii)

The conditional mean function $b(\eta)=\mathbb{E}[\nu\mid\eta]$ is $3$ times
differentiable with $b^{\prime\prime}$ and $b^{\prime\prime\prime}$ are
uniformly bounded over real line.

  3. iii)

The variance function $\sigma^{2}(\eta)={\sf var}({\epsilon}\mid\eta)$ is
uniformly bounded from above.

  4. iv)

There exists some $\xi>0$ such that: $\min_{|x|\leq\tau+\xi}f_{\eta}(x)>0$.

######  Assumption 3.3.

Define the matrices $\Omega$ and $\Omega^{*}$ as:

$\displaystyle\Omega$ $\displaystyle=\mathbb{E}\left[{\sf
var}\left(\begin{bmatrix}S&X&Zb^{\prime}(\eta)\end{bmatrix}\left|\right.\eta\right)\right]+{\sf
var}\left(\begin{bmatrix}0&0&Z\end{bmatrix}\right)$ $\displaystyle\Omega^{*}$
$\displaystyle=\mathbb{E}\left[\sigma^{2}(\eta){\sf
var}\left(\begin{bmatrix}S&X&Zb^{\prime}(\eta)\end{bmatrix}\mid\eta\right)\right]+{\sf
var}(\eta){\sf var}\left(\begin{bmatrix}0&0&Z\end{bmatrix}\right)\,.$

Then, the minimum eigenvalues of $\Omega$ and $\Omega^{*}$ are strictly
positive.

######  Assumption 3.4.

The distribution of $(X,Z)$ satisfies the following conditions:

  1. i)

$(X,Z)$ has bounded continuous density function and have zero expectation.

  2. ii)

The first four moments of $X,Z$ are finite.

######  Remark 3.5.

Assumption 3.2 provides a low-level assumption on the smoothness of the
density of $\eta$, the conditional mean function $b(\eta)$ and the conditional
variance profile $\sigma^{2}(\eta)$, which is required for the standard
asymptotic analysis of the partial linear model. Assumption 3.3, again is a
standard assumption in partial linear model literature. It is essential for
the asymptotic normality of the treatment effect as the asymptotic variance of
our estimator is a function of these variances. If this assumption is
violated, then the asymptotic variance will be infinite and that estimation at
$\sqrt{n}$ is not possible. (Note that if $\gamma_{0}=0$, then Assumption 3.3
is violated.) As our method does not use all observations, but a fraction
depending on the interval $[-\tau,\tau]$, our limiting variance comprises of
the following truncated versions of $\Omega$ and $\Omega^{*}$:

$\displaystyle\Omega_{\tau}$ $\displaystyle=\mathbb{E}\left[{\sf
var}\left(\begin{bmatrix}S&X&Zb^{\prime}(\eta)\end{bmatrix}\mid\eta\right)\mathds{1}_{|\eta|\leq\tau}\right]+{\sf
var}\left(\begin{bmatrix}0&0&Z\end{bmatrix}\right)\,,$
$\displaystyle\Omega^{*}_{\tau}$
$\displaystyle=\mathbb{E}\left[\sigma^{2}(\eta){\sf
var}\left(\begin{bmatrix}S&X&Zb^{\prime}(\eta)\end{bmatrix}\mid\eta\right)\mathds{1}_{|\eta|\leq\tau}\right]+{\sf
var}(\eta){\sf var}\left(\begin{bmatrix}0&0&Z\end{bmatrix}\right)\,.$

It is immediate that if $\tau\to\infty$, then $\Omega_{\tau}\to\Omega$ and
$\Omega^{*}_{\tau}\to\Omega^{*}$. Hence, in light of Assumption 3.3, by
continuity, the minimum eigenvalues of $\Omega_{\tau}$ and $\Omega^{*}_{\tau}$
are also positive for all large $\tau$.

###  3.1 Asymptotic normality

We now state the main result of the paper:

######  Theorem 3.6.

Consider the estimates obtained at the end of the previous section. Under
assumptions 3.1-3.4:

$\sqrt{n}\left(\hat{\alpha}-\alpha_{0}\right)\overset{\mathscr{L}}{\implies}\mathcal{N}(0,3e_{1}^{\top}\Omega_{\tau}^{-1}\Omega^{*}_{\tau}\Omega_{\tau}^{-1}e_{1})\,,$

whilst

$\sqrt{n}\left(\bar{\hat{\alpha}}-\alpha_{0}\right)\overset{\mathscr{L}}{\implies}\mathcal{N}(0,e_{1}^{\top}\Omega_{\tau}^{-1}\Omega^{*}_{\tau}\Omega_{\tau}^{-1}e_{1})$

Sketch of proof: We present a high-level outline of the key steps of the
proof, deferring all technical details to Subsection B of the Supplementary
document. From (2), on a set of probability approaching 1, our estimator can
be written as:

$\displaystyle\hat{\alpha}$
$\displaystyle=e_{1}^{\top}\left(\mathbf{W}^{\top}{\sf
proj}^{\perp}_{\tilde{\mathbf{N}}_{K,a}}\mathbf{W}\right)^{-1}\mathbf{W}^{\top}{\sf
proj}^{\perp}_{\tilde{\mathbf{N}}_{K,a}}\begin{pmatrix}\mathbf{Y}\\\
\mathbf{Q}-\mathbf{Z}\hat{\gamma}_{n}\end{pmatrix}$
$\displaystyle=\alpha_{0}+e_{1}^{\top}\left(\mathbf{W}^{\top}{\sf
proj}^{\perp}_{\tilde{\mathbf{N}}_{K,a}}\mathbf{W}\right)^{-1}\mathbf{W}^{\top}{\sf
proj}^{\perp}_{\tilde{\mathbf{N}}_{K,a}}\left[\begin{pmatrix}\mathbf{R}\\\
0\end{pmatrix}+\begin{pmatrix}{\epsilon}\\\ \eta\end{pmatrix}\right]\,.$

This implies:

$\sqrt{n}\left(\hat{\alpha}-\alpha_{0}\right)=e_{1}^{\top}\left(\frac{\mathbf{W}^{\top}{\sf
proj}^{\perp}_{\tilde{\mathbf{N}}_{K,a}}\mathbf{W}}{n}\right)^{-1}\frac{\mathbf{W}^{\top}{\sf
proj}^{\perp}_{\tilde{\mathbf{N}}_{K,a}}}{\sqrt{n}}\left[\begin{pmatrix}\mathbf{R}\\\
0\end{pmatrix}+\begin{pmatrix}{\epsilon}\\\ \eta\end{pmatrix}\right]$ (3.1)

which is our main estimating equation. We next outline the key steps of our
proof.

Step 1:  First show that:

$\frac{\mathbf{W}^{\top}{\sf
proj}^{\perp}_{\tilde{\mathbf{N}}_{K,a}}\mathbf{W}}{n}\overset{P}{\longrightarrow}\frac{1}{3}\Omega_{\tau}\,.$

where $\Omega_{\tau}$ is as defined in Remark 3.5.

Step 2: Next, establish the following asymptotic linear expansion:

$\frac{\mathbf{W}^{\top}{\sf
proj}^{\perp}_{\tilde{\mathbf{N}}_{K,a}}}{\sqrt{n}}\begin{pmatrix}{\bm{\epsilon}}\\\
\eta\end{pmatrix}=\frac{1}{\sqrt{n}}\sum_{i=1}^{n/3}\varphi\left(X_{i},Z_{i},\eta_{i},\nu_{i}\right)\
+o_{p}(1)\,.$

for some influence function $\varphi$.

Step 3:  Apply central limit theorem to obtain:

$\frac{1}{\sqrt{n}}\sum_{i=1}^{n/3}\varphi\left(X_{i},Z_{i},\eta_{i},\nu_{i}\right)\overset{\mathscr{L}}{\implies}\mathcal{N}(0,\frac{1}{3}\Omega^{*}_{\tau})\,.$

where $\Omega^{*}_{\tau}$ is as defined in Remark 3.5.

Step 4: Finally ensure that ‘residual term’ is asymptotically negligible:

$\frac{\mathbf{W}^{\top}{\sf
proj}^{\perp}_{\mathbf{N}_{K,a}}}{\sqrt{n}}\begin{pmatrix}\mathbf{R}\\\
0\end{pmatrix}\overset{P}{\longrightarrow}0\,.$

Now, combining above steps we can conclude:

$\displaystyle\sqrt{n}(\hat{\alpha}-\alpha_{0})$
$\displaystyle=e_{1}^{\top}\left(\frac{1}{3}\Omega_{\tau}\right)^{-1}\frac{1}{\sqrt{n}}\sum_{i=1}^{n/3}\varphi\left(X_{i},Z_{i},\eta_{i},\nu_{i}\right)+o_{p}(1)$
$\displaystyle\overset{\mathscr{L}}{\implies}\mathcal{N}\left(0,3e_{1}^{\top}\Omega_{\tau}^{-1}\Omega^{*}_{\tau}\Omega^{-1}_{\tau}e_{1}\right)$

where the leading term only depends on the observations in $\mathcal{D}_{3}$
and is consequently independent of $\mathcal{D}_{1},\mathcal{D}_{2}$. Finally,
rotating the dataset and taking average of the $\hat{\alpha}$’s we further
conclude:

$\sqrt{n}(\bar{\hat{\alpha}}-\alpha_{0})\overset{\mathscr{L}}{\implies}\mathcal{N}\left(0,e_{1}^{\top}\Omega_{\tau}^{-1}\Omega^{*}_{\tau}\Omega^{-1}_{\tau}e_{1}\right)\,.$

###  3.2 Semi-parametric efficiency

We further show that our estimator is semi-parametrically efficient under
certain restrictions. As our estimator is based on least square approach, it
can not be shown to be efficient unless the error ${\epsilon}$ is normal. We
prove the following theorem (proof can be found in Section C.2 of the
supplementary document):

######  Theorem 3.7.

Suppose the model is the following:

$Y_{i}=\alpha_{0}\mathds{1}_{Z_{i}^{\top}\gamma_{0}+\eta>0}+X_{i}^{\top}\beta_{0}+b(\eta_{i})+{\epsilon}_{i}\,.$

where ${\epsilon}_{i}\sim\mathcal{N}(0,\tau^{2})\perp\perp\eta_{i}$. Then our
estimator of $\alpha_{0}$ is semi-parametrically efficient, i.e. its
asymptotic variance $\sigma_{0}^{2}=\tau^{2}e_{1}^{\top}\Omega^{-1}e_{1}$
attains the semi-parametric information bound for this model.

######  Remark 3.8.

The assumption of the normality of ${\epsilon}$ is necessary to establish
semiparametric efficiency for least squares type methods, but the assumption
of homoskedasticity is essential only if we use ordinary least squares method.
One may easily take care of heteroskedasticity by using weighted least squares
instead. The first step towards that direction is to approximate the variance
profile $\sigma(\eta)$ using $\hat{\sigma}(\hat{\eta})$ for some non-
parametric estimate $\hat{\sigma}(\cdot)$ of $\sigma(\cdot)$. Then defining
$\mathbf{D}\in\mathbb{R}^{(n_{3}+n/3)\times(n_{3}+n/3)}$ to be the diagonal
matrix with first $n/3$ diagonal entries being
$\hat{\sigma}(\hat{\eta}_{i})$’s (i.e. for all those $\hat{\eta}_{i}$’s such
that $|\hat{\eta}_{i}|\leq\tau$) and last $n/3$ diagonal entries being
$\hat{\sigma}_{\eta}$’s (an estimate of variance of $\eta$) we estimate the
treatment effect as:

$\hat{\alpha}=e_{1}^{\top}\left(\mathbf{W}^{\top}\mathbf{D}^{-1/2}{\sf
proj}^{\perp}_{\mathbf{D}^{-1/2}\tilde{\mathbf{N}}_{K,a}}\mathbf{D}^{-1/2}\mathbf{W}\right)^{-1}\mathbf{W}^{\top}\mathbf{D}^{-1/2}{\sf
proj}^{\perp}_{\mathbf{D}^{-1/2}\tilde{\mathbf{N}}_{K,a}}\mathbf{D}^{-1/2}\begin{pmatrix}\mathbf{Y}\\\
\mathbf{Q}-\mathbf{Z}\hat{\gamma}_{n}\end{pmatrix}$

A more tedious analysis establishes that this estimator is asymptotic normal
and semi-parametrically efficient under the error structure:
$\nu_{i}=b(\eta_{i})+\sigma(\eta_{i}){\epsilon}_{i}$, where
${\epsilon}_{i}\sim\mathcal{N}(0,1)$. As this does not add anything of
significance to the core idea of the paper, we confine ourselves to use OLS
instead of WLS for ease of presentation.

##  4 Analysis for high dimensional covariates

As highlighted in our analysis for fixed dimensional covariates in the
previous section, our model is an example of a partial linear model, where the
observations corresponding to the non-linear part are noisy and the noise is
correlated with the covariates corresponding to the linear part. Recall that
our is model is defined as:

$\displaystyle Y_{i}$ $\displaystyle=\alpha_{0}\mathds{1}_{Q_{i}\geq
0}+X_{i}^{\top}\beta_{0}+b(\eta_{i})+{\epsilon}_{i}$ $\displaystyle Q_{i}$
$\displaystyle=Z_{i}^{\top}\gamma_{0}+\eta_{i}\,.$ (4.1)

with $\mathbb{E}[\eta\mid Z]=0$ and $\mathbb{E}[{\epsilon}\mid\eta,X,Z]=0$. In
this section, we assume that the dimension of ($X,Z$) is larger than the
sample size. More specifically, denoting $p_{1}:=\dim(X)$ and
$p_{2}:=\dim(Z)$, we assume that $p_{1}\wedge p_{2}\gg n$. As before, we only
observe $(Y_{i},X_{i},Z_{i},Q_{i})$ but not $\eta_{i}$. As in the case for
almost all high dimensional statistical analysis, we assume that both
$\beta_{0}$ and $\gamma_{0}$ are sparse vectors (with
$\|\beta_{0}\|_{0}=s_{\beta}$ and $\|\gamma_{0}\|_{0}=s_{\gamma}$ where
$s_{\beta}\vee s_{\gamma}\ll n$) to enable consistent estimation of the
treatment effect $\alpha_{0}$. We will quantify the precise assumptions needed
for our theory later. The core estimation procedure in this high dimensional
regime is similar to that for its fixed dimension counterpart, but with
certain changes to take care of the effect of high dimensional covariates. As
before, we divide the whole data $\mathcal{D}$ into three parts
$\mathcal{D}_{1},\mathcal{D}_{2},\mathcal{D}_{3}$ (with $n_{i}$ data in
$\mathcal{D}_{i}$, where $n_{1}\sim n_{2}\sim n_{3}\sim n/3$). From
$\mathcal{D}_{1}$, we estimate $\gamma_{0}$ by performing LASSO regression of
$Q$ on $Z$:

$\hat{\gamma}=\mathop{\rm
argmin}_{\gamma\in{\mathbb{R}}^{p_{2}}}\left[\frac{1}{2n_{3}}\left\|\mathbf{Q}-\mathbf{Z}\gamma\right\|^{2}+\lambda\left\|\gamma\right\|_{1}\right]\,,$

where $\lambda\asymp\sqrt{\log{p_{2}}/n}$ and set
$\hat{\eta}_{i}=Q_{i}-Z_{i}^{\top}\hat{\gamma}$ for all observations in
$\mathcal{D}_{2}\cup\mathcal{D}_{3}$. This estimator of $\gamma_{0}$ is
consistent and rate optimal under the restricted eigenvalue assumption
(henceforth RE) on $\mathbf{Z}$. Just like in our analysis for fixed
dimension, we only consider all the $\hat{\eta}$’s (both in $\mathcal{D}_{2}$
and $\mathcal{D}_{3}$) such that $|\hat{\eta}|\leq\tau$ and ignore the
remaining data from $\mathcal{D}_{2}$ and $\mathcal{D}_{3}$. Therefore, the
rest of our analysis on $\mathcal{D}_{2}$ and $\mathcal{D}_{3}$ are solely
based on those observation for which $|\hat{\eta}|\leq\tau$. Using this
approximation of the unknown $\eta_{i}$’s, we expand the first equation of
(4.1) as:

$Y_{i}=\alpha_{0}\mathds{1}_{Q_{i}\geq
0}+X_{i}^{\top}\beta_{0}+b(\hat{\eta}_{i})+\underbrace{(\eta_{i}-\hat{\eta}_{i})}_{Z_{i}^{\top}(\hat{\gamma}_{n}-\gamma_{0})}b^{\prime}(\hat{\eta}_{i})+(\eta_{i}-\hat{\eta}_{i})^{2}b^{\prime\prime}(\tilde{\eta}_{i})+{\epsilon}_{i}\,.$
(4.2)

The logic behind this two step Taylor expansion is similar to that for our
fixed dimensional analysis as articulated in Section 2. We next estimate
$b^{\prime}$ using $\mathcal{D}_{2}$ using a different expansion of (4.1).
Equation (4.2) will be used later for estimating $\alpha_{0}$ based on the
data in $\mathcal{D}_{3}$. To estimate $b^{\prime}$ (using the observations in
$\mathcal{D}_{2}$) we expand the first equation of (4.1) as:

$Y_{i}=\alpha_{0}\mathds{1}_{Q_{i}\geq
0}+X_{i}^{\top}\beta_{0}+b(\hat{\eta}_{i})+R_{i}+{\epsilon}_{i}\,.$ (4.3)

The above equation is obtained simply by replacing $b(\eta_{i})$ in the first
equation of (4.1) with $b(\hat{\eta}_{i})$ giving the residual term
$R_{i}=b(\hat{\eta}_{i})-b(\eta_{i})$. We next invoke techniques from the
analysis of high dimensional partial linear models to estimate $b^{\prime}$
using the above equation. Replacing $b^{\prime}(\hat{\eta}_{i})$ by
$\hat{b}^{\prime}(\hat{\eta}_{i})$ in equation (4.2) we obtain the following
representation for the observations in $\mathcal{D}_{3}$:

$\displaystyle Y_{i}$ $\displaystyle=\alpha_{0}\mathds{1}_{Q_{i}\geq
0}+X_{i}^{\top}\beta_{0}+b(\hat{\eta}_{i})+Z_{i}^{\top}(\hat{\gamma}_{n}-\gamma_{0})\hat{b}^{\prime}(\hat{\eta}_{i})$
$\displaystyle\qquad\qquad+Z_{i}^{\top}(\hat{\gamma}_{n}-\gamma_{0})\left(b^{\prime}(\hat{\eta}_{i})-\hat{b}^{\prime}(\hat{\eta}_{i})\right)+(\eta_{i}-\hat{\eta}_{i})^{2}b^{\prime\prime}(\tilde{\eta}_{i})+{\epsilon}_{i}\,.$

Further, an approximation of $b(\hat{\eta}_{i})$ by B-spline basis yields:

$Y_{i}=\alpha_{0}\mathds{1}_{Q_{i}\geq
0}+X_{i}^{\top}\beta_{0}+\mathbf{N}_{k}(\hat{\eta}_{i})^{\top}\omega_{b}+Z_{i}^{\top}(\hat{\gamma}_{n}-\gamma_{0})\hat{b}^{\prime}(\hat{\eta}_{i})+\tilde{R}_{i}+{\epsilon}_{i}$
(4.4)

where $\mathbf{N}_{k}(\eta)$ consists of the B-spline basis functions
evaluated at $\hat{\eta}_{i}$. The residual term $\tilde{R}_{i}$ can be
decomposed as $R_{1,i}+R_{2,i}+R_{3,i}$ where the individual residuals are:

$\displaystyle R_{1,i}$
$\displaystyle=b(\hat{\eta}_{i})-\mathbf{N}_{k}(\hat{\eta}_{i})^{\top}\omega_{b}$
$\displaystyle R_{2,i}$
$\displaystyle=Z_{i}^{\top}(\hat{\gamma}_{n}-\gamma_{0})\left(b^{\prime}(\hat{\eta}_{i})-\hat{b}^{\prime}(\hat{\eta}_{i})\right)$
$\displaystyle R_{3,i}$
$\displaystyle=(\eta_{i}-\hat{\eta}_{i})^{2}b^{\prime\prime}(\tilde{\eta}_{i})\,.$

The first residual is the B-spline approximation error of $b$ while the second
is the product of two different error terms: (i) the error in estimation of
$\gamma_{0}$ from $\mathcal{D}_{3}$ and (ii) the error in estimation of
$b^{\prime}$ from $\mathcal{D}_{2}$ and the last is the Taylor approximation
error. We use equation (4.4) as our main estimating equation for $\alpha_{0}$.
Before delving into the estimation procedure we introduce some notation:

  1. 1.

$\breve{Z}_{i}$ will be used to denote
$Z_{i}\hat{b}^{\prime}(\hat{\eta}_{i})$.

  2. 2.

The vector $W$ will be used to denote the random vector
$(S,X,Zb^{\prime}(\eta))$.

  3. 3.

$\breve{W}$ will be used to denote the random vector $(S,X,\breve{Z})$ and the
$j^{th}$ element of $W$ (resp. $\check{W}$) will be denoted as $W_{j}$ (resp.
$\check{W}_{j}$).

  4. 4.

$m_{j}(\eta)$ will be used to denote $\mathbb{E}[W_{j}\mid\eta]$, for all
$0\leq j\leq 1+p_{1}+p_{2}$ where $W_{0}=Y$.

  5. 5.

$\check{m}_{j}(\hat{\eta})$ will be used to denote
$\mathbb{E}[W_{j}\mid\hat{\eta}]$, for all $0\leq j\leq 1+p_{1}+p_{2}$ where
$W_{0}=Y$.

  6. 6.

$\widehat{\tilde{W}}$ will be used to denote $W-\mathbb{E}[W\mid\hat{\eta}]$.

  7. 7.

$\tilde{W}$ will be used to denote the random vector
$W-\mathbb{E}[W\mid\eta]$.

  8. 8.

$\breve{\mathbf{W}},\hat{\tilde{\mathbf{W}}},\tilde{\mathbf{W}}$ will be used
to denote the matrix version of $\breve{W},\hat{\tilde{W}},\tilde{W}$
respectively by concatenating all of the $n_{3}$ observations of
$\mathcal{D}_{3}$ row-wise. For any matrix $\mathbf{A}$, we use
$\mathbf{A}_{*,i}$ to denote the $i^{th}$ row of $\mathbf{A}$, and
$\mathbf{A}_{*,j}$ to denote the $j^{th}$ column of $\mathbf{A}$.

  9. 9.

$\theta_{0}\equiv\theta_{0,n}=(\alpha_{0},\beta_{0},\hat{\gamma}_{n}-\gamma_{0})$.

  10. 10.

$\mathbf{N}_{k}$ will be used to denote the $n\times K$ matrix, whose $i^{th}$
row consists of $K$ B-spline bases evaluated at $\hat{\eta}_{i}$.

  11. 11.

For any vector $v$ (or matrix $\mathbf{A}$), the notation $v^{\perp}$ (resp.
$\mathbf{A}^{\perp}$) is used to denote $P^{\perp}_{\mathbf{N}_{k}}v$ (resp.
$P^{\perp}_{\mathbf{N}_{k}}A$).

Using the above notation, equation (4.4) can be rewritten (in matrix form
concatenating all the observations along the rows) as:

$\mathbf{Y}=\breve{\mathbf{W}}\theta_{0}+\mathbf{N}_{k}\omega_{b}+\tilde{\mathbf{R}}+{\bm{\epsilon}}\,.$

Projecting out the effect of $\mathbf{N}_{k}$ from both sides yields:

$P^{\perp}_{\mathbf{N}_{k}}\mathbf{Y}=P^{\perp}_{\mathbf{N}_{k}}\breve{\mathbf{W}}\theta_{0}+P^{\perp}_{\mathbf{N}_{k}}\tilde{\mathbf{R}}+P^{\perp}_{\mathbf{N}_{k}}{\bm{\epsilon}}\,.$

Henceforth we use $\perp$ notation as superscript to denote that the
vector/matrix is pre-multiplied by $P^{\perp}_{\mathbf{N}_{k}}$. Our estimate
of $\alpha_{0}$ is defined as follows:

$\hat{\alpha}=\frac{\left(\mathbf{Y}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,Y}\right)^{\top}\left(\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,S}\right)}{\left(\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,S}\right)^{\top}\left(\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,S}\right)}$

where $\hat{\theta}_{-1,Y}$ denotes the LASSO estimator obtained by regressing
$\mathbf{Y}$ on $\breve{\mathbf{W}}_{-1}$, i.e.

$\hat{\theta}_{-1,Y}=\mathop{\rm
argmin}_{\theta}\left[\frac{1}{2n_{3}}\left\|\mathbf{Y}-\breve{\mathbf{W}}_{-1}\theta\right\|^{2}+\lambda_{0}\left\|\theta\right\|_{1}\right]$
(4.5)

and $\hat{\theta}_{-1,S}$ denotes the LASSO estimator obtained by regressing
$\mathbf{S}$ on $\breve{\mathbf{W}}_{-1}$, i.e.

$\hat{\theta}_{-1,S}=\mathop{\rm
argmin}_{\theta}\left[\frac{1}{2n_{3}}\left\|\mathbf{S}-\breve{\mathbf{W}}_{-1}\theta\right\|^{2}+\lambda_{1}\left\|\theta\right\|_{1}\right]$
(4.6)

for some appropriate choice of $\lambda_{0},\lambda_{1}$ to be specified
later. Here $n_{3}$ is the number of observations in $\mathcal{D}_{3}$ for
which $|\hat{\eta}_{i}|\leq\tau$. It is immediate that $n_{3}\asymp n/3\asymp
n$, i.e. the number of observations in $\mathcal{D}_{3}$ with
$|\hat{\eta}_{i}|\leq\tau$ is of the order $n$, the total number of
observations. Henceforth, we will ignore this difference and state all our
results in terms of $n$. We show that
$\sqrt{n_{3}}\left(\hat{\alpha}-\alpha_{0}\right)$ is asymptotically normal
under certain assumptions which are stated below:

######  Assumption 4.1 (Smoothness of $b$).

The function $b$ in equation (4.1) is assumed to be $\upsilon\geq 3$ times
differentiable with bounded derivates.

######  Assumption 4.2.

The density $f_{\eta}$ of $\eta$ is bounded, continuously differentiable and
lower bounded by some $f_{-}>0$ on $[-\tau-\xi,\tau+\xi]$ for some small
$\xi>0$.

######  Assumption 4.3 (Smoothness of conditional expectation).

Define the function $g_{j}(a,t)$ as:

$g_{j}(a,t)=\mathbb{E}\left[W_{j}\mid\eta+a^{\top}Z=t\right]$

for $0\leq j\leq 1+p_{1}+p_{2}$. Assume for any fixed $\|a\|\leq r$, the
collection $\\{g_{j}(a,\cdot)\\}_{1\leq j\leq p}$ belongs to a function class
$\Sigma(\alpha,\mathbf{l})$ for some $\alpha\geq 3$ and
$\mathbf{l}\in{\mathbb{R}}^{\lfloor\alpha\rfloor}$, where
$\Sigma(\alpha,\mathbf{l})$ is defined as the collection of all the functions
$f$, which are $\lfloor\alpha\rfloor$ times differentiable and the
$\lfloor\alpha\rfloor^{th}$ derivative satisfies:

$\left|f^{\lfloor\alpha\rfloor}(x)-f^{\lfloor\alpha\rfloor}(y)\right|\leq\mathbf{l}_{\lfloor\alpha\rfloor}\left|x-y\right|^{\alpha-\lfloor\alpha\rfloor}$

and $\|f^{(i)}\|_{\infty}\leq\mathbf{l}_{i}$ for all $1\leq
i\leq\lfloor\alpha\rfloor-1$. Moreover, this $\mathbf{l}$ is uniform over all
$\|a\|\leq r$ where $r$ is independent of the underlying dimension.

######  Remark 4.4.

The conditional mean functions $m_{j}$ and $\check{m}_{j}$ are special cases
of $g_{j}$ as $\check{m}_{j}(t)=g_{j}(\hat{\gamma}_{n}-\gamma_{0},t)$ and
$m_{j}(t)=g_{j}(0,t)$.

######  Assumption 4.5 (Sub-gaussianity).

Assume that for $0\leq j\leq p_{1}+p_{2}$,
$\tilde{W}_{j},\hat{\tilde{W}}_{j},m_{j}(\eta),\check{m}_{j}(\hat{\eta}),\eta$
and ${\epsilon}$ are subgaussian random-variables with subgaussianity constant
uniformly bounded by $\sigma_{W}$. Furthermore assume that ${\sf
var}(W_{j}\mid\hat{\eta})$ is uniformly (over $j$) bounded on $[-\tau,\tau]$.

For our next assumption, we define a covariance matrix $\Sigma$ and two
vectors $\theta^{*}_{Y}\in{\mathbb{R}}^{p_{1}+p_{2}}$ and
$\theta^{*}_{S}\in{\mathbb{R}}^{p_{1}+p_{2}}$ which play a crucial role in our
analysis:

$\displaystyle\Sigma_{\tau}$
$\displaystyle=\mathbb{E}_{\eta}\left[\left(W-\mathbb{E}[W\mid\eta]\right)\left(W-\mathbb{E}[W\mid\eta]\right)^{\top}\mathds{1}_{|\eta|\leq\tau}\right]$
$\displaystyle=\mathbb{E}[\tilde{W}\tilde{W}^{\top}\mathds{1}_{|\eta|\leq\tau}]$
$\displaystyle=\mathbb{E}_{\eta}\left[{\sf
var}(W\mid\eta)\mathds{1}_{|\eta|\leq\tau}\right]$
$\displaystyle=\mathbb{E}_{\eta}\left[{\sf
var}\left(\begin{pmatrix}S&X&b^{\prime}(\eta)Z\end{pmatrix}\mid\eta\right)\mathds{1}_{|\eta|\leq\tau}\right]$
$\displaystyle\theta^{*}_{S}$ $\displaystyle=\mathop{\rm
argmin}_{\delta}\mathbb{E}\left[\left(\tilde{S}-\tilde{W}_{-1}^{\top}\delta\right)^{2}\mathds{1}_{|\eta|\leq\tau}\right]=\left(\Sigma_{\tau_{-1,-1}}\right)^{-1}\mathbb{E}[\tilde{W}_{-1}\tilde{S}\mathds{1}_{|\eta|\leq\tau}]$
(4.7) $\displaystyle\theta^{*}_{Y}$ $\displaystyle=\mathop{\rm
argmin}_{\delta}\mathbb{E}\left[\left(\tilde{Y}-\tilde{W}_{-1}^{\top}\delta\right)^{2}\mathds{1}_{|\eta|\leq\tau}\right]=\left(\Sigma_{\tau_{-1,-1}}\right)^{-1}\mathbb{E}[\tilde{W}_{-1}\tilde{Y}\mathds{1}_{|\eta|\leq\tau}]\,.$
(4.8)

The matrix $\Sigma_{\tau}$ is the expectation of conditional variance of $W$
given $\eta$ on the set $|\eta|\leq\tau$, which arises in the estimation of
the linear part of a partial linear model. This matrix can be thought as high
dimensional analogue of $\Omega_{\tau}$ defined in our analysis for fixed
dimensional covariates. $\theta_{Y}^{*}$ (resp. $\hat{\theta}^{*}_{S}$) is the
best linear estimator in the population for regressing $\tilde{Y}$ (resp.
$\tilde{S}$) on $\tilde{W}$. In the definition of our estimator
$\hat{\alpha}$, we have regressed $\mathbf{Y}^{\perp}$ on
$\breve{\mathbf{W}}_{-1}^{\perp}$ and $\mathbf{S}^{\perp}$ on
$\breve{\mathbf{W}}_{-1}^{\perp}$. Intuitively speaking, projecting out the
column space of $\mathbf{N}_{k}(\hat{\eta})$ is asymptotically equivalent to
centering around the conditional expectation with respect to $\eta$.
Therefore, it is expected that $\hat{\theta}_{-1,Y}$ (respectively
$\hat{\theta}_{-1,S}$) should be asymptotically consistent for regressing
$\theta_{Y}^{*}$ (respectively $\hat{\theta}^{*}_{S})$ under certain sparsity
assumption and Restricted eigenvalue (RE) condition on the covariate matrix
$\breve{\mathbf{W}}_{-1}^{\perp}$. We now state our assumptions on
$\Sigma,\theta^{*}_{Y},\theta^{*}_{S}$:

######  Assumption 4.6 (Asymptotic variance).

Assume that there exists $C_{\min}>0$ and $C_{\max}<\infty$ such that:

$C_{\min}\leq\lambda_{\min}\left(\Sigma_{\tau}\right)\leq\lambda_{\max}\left(\Sigma_{\tau}\right)\leq
C_{\max}<\infty\,.$

Furthermore, define

$\sigma^{2}_{n,1}=\mathbb{E}\left[{\epsilon}^{2}\left(\tilde{S}-\tilde{W}_{-1}^{\top}\theta^{*}_{S}\right)^{2}\mathds{1}_{|\eta|\leq\tau}\right]$

and assume that:

$C_{\min}\leq\liminf_{n}\sigma^{2}_{n,1}\leq\limsup_{n}\sigma^{2}_{n,1}\leq
C_{\max}<\infty\,.$

Furthermore assume that:

$\limsup_{n\to\infty}\mathbb{E}\left[\left(\tilde{S}-\tilde{W}_{-1}^{\top}\theta^{*}_{S}\right)^{2+\xi}\right]<\infty\,.$

for some small $\xi>0$.

######  Assumption 4.7 (Sparsity).

Assume there exists $s_{0},s_{1}>0$ such that $\|\theta^{*}_{Y}\|_{0}\leq
s_{0}$, $\|\theta^{*}_{S}\|_{0}\leq s_{1}$ and $s_{i}\log{(p_{1}\wedge
p_{2})}/n\to 0$ .

######  Remark 4.8.

Although it is apparent from equation (4.5) and (4.6) that RE condition is
needed on $\breve{\mathbf{W}}_{-1}^{\perp}$ for the consistency and rate
optimality of LASSO estimates $\hat{\theta}_{-1,S}$ and $\theta^{*}_{-1,1}$,
from the sub-gaussianity of $\tilde{W}$, we have that
$\tilde{\mathbf{W}}_{-1}$ satisfies RE condition with high probability (see
[38]). We prove in Proposition A.5 that this implies
$\breve{\mathbf{W}}_{-1}^{\perp}$ also satisfies RE condition with high
probability and consequently $\hat{\theta}_{-1,S}$ and $\hat{\theta}_{-1,Y}$
will be rate efficient.

######  Remark 4.9.

From our original model equation 4.1, we have:

$\tilde{Y}=Y-\mathbb{E}[Y\mid\eta]=\tilde{S}\alpha_{0}+\tilde{X}\beta_{0}\,,$

which implies

$\theta^{*}_{Y}=\left(\Sigma_{\tau_{-1,-1}}\right)^{-1}\mathbb{E}[\tilde{W}_{-1}\tilde{Y}\mathds{1}_{|\eta|\leq\tau}]=\theta^{*}_{S}\alpha_{0}+\theta_{0,-1}^{*}$

where $\theta_{0,-1}^{*}=(\beta_{0},0)$. Note that $\theta^{*}_{0,-1}$ is
already sparse with $\left\|\theta_{0,-1}^{*}\right\|_{0}=s_{\beta}$ (recall
that we have defined $s_{\beta}=\|\beta_{0}\|_{0}$). So sparsity assumption on
$\theta^{*}_{S}$ automatically ensures sparsity of $\theta^{*}_{Y}$, in other
words we have $s_{0}\leq s_{1}+s_{\beta}$.

Under the above assumptions we prove the following theorem:

######  Theorem 4.10.

Under the above assumptions and certain conditions on the sparsities:
$s_{0},s_{1},s_{\beta},s_{\gamma}$ (see subsection A.4 for detailed
discussion), we have:

$\frac{\sigma^{2}_{n,2}}{\sigma_{n,1}}\sqrt{n_{3}}\left(\hat{\alpha}-\alpha_{0}\right)\overset{\mathscr{L}}{\implies}\mathcal{N}(0,1)\,.$

where the values of $\sigma_{n,1}$ is as defined in Assumption 4.6 and
$\sigma_{n,2}$ is:

$\displaystyle\sigma_{n,2}$
$\displaystyle=\sqrt{\mathbb{E}\left[\left(\tilde{S}-\tilde{W}_{-1}\theta^{*}_{S}\right)^{2}\mathds{1}_{|\eta|\leq\tau}\right]}=\frac{1}{\sqrt{\left(\Sigma_{\tau}^{-1}\right)_{1,1}}}\,.$

######  Remark 4.11.

Similar to our analysis for the fixed dimensional covariates, we can gain
efficiency (in terms of asymptotic variance) here as well, by rotating the
datasets and taking average of the estimates of $\alpha_{0}$). However, as the
analysis is already quite involved, we do not pursue this extension here.

A roadmap of the proof: We now present a basic roadmap of the proof of Theorem
4.10 for the ease of the readers, the details can be found in Appendix A and
in the supplementary document. As immediate from the definition, the
estimation of $\hat{\alpha}$ consists of three key steps:

  1. 1.

LASSO regression using $\mathbf{Y}^{\perp}$ on
$\breve{\mathbf{W}}_{-1}^{\perp}$.

  2. 2.

LASSO regression using $\mathbf{S}^{\perp}$ on
$\breve{\mathbf{W}}_{-1}^{\perp}$.

  3. 3.

Finally, regression of the residual of the first LASSO regression on the
residual of the second LASSO regression.

Recall that the last $p_{2}$ columns of $\breve{\mathbf{W}}$ involves
$\hat{b}^{\prime}(\hat{\eta})$ as a coefficient of $Z$, which, asymptotically
should be close to $b^{\prime}(\eta)$, i.e. the random vector $\breve{W}$
should be asymptotically close to the random vector $W$. Consequently, if we
consider $\breve{W}^{\perp}$, it should be asymptotically close to
$\tilde{W}$, as projecting out the span of the b-spline bases evaluated at
$\hat{\eta}$ is expected to be asymptotically equivalent to centering around
$\eta$. Using this intuition, we expect that our LASSO estimates
$\hat{\theta}_{-1,Y}$ (first LASSO) should converge to $\theta^{*}_{Y}$ and
the estimator $\hat{\theta}_{-1,S}$ (second LASSO) should converge to
$\theta^{*}_{S}$ under our assumptions. This is what we establish using Lemma
A.4 and Proposition A.5. Lemma A.3, which establishes the _asymptotic
closeness_ of $\breve{\mathbf{W}}^{\perp}$ to $\tilde{\mathbf{W}}$ via the
following approximation:

$\breve{\mathbf{W}}^{\perp}\longrightarrow\mathbf{W}^{\perp}\longrightarrow\widehat{\tilde{\mathbf{W}}}\longrightarrow\tilde{\mathbf{W}}\,,$

is an auxiliary lemma that is used in the proof of Lemma A.4. Furthermore, for
the optimal rate of LASSO estimator, we need restricted eigenvalue condition
of the covariate matrix $\breve{\mathbf{W}}_{-1}^{\perp}$. In Assumption 4.7
we have assumed the matrix $\tilde{\mathbf{W}}_{-1}$ satisfied RE condition.
In Proposition A.5 we establish that, if $\tilde{\mathbf{W}}_{-1}$ satisfies
RE, then $\breve{\mathbf{W}}_{-1}^{\perp}$ also satisfies it with high
probability. To provide further insights of the proof, we expand
$\sqrt{n}(\hat{\alpha}-\alpha_{0})$ as follows:

$\displaystyle\hat{\alpha}$
$\displaystyle=\frac{\left(\mathbf{Y}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,Y}\right)^{\top}\left(\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,S}\right)}{\left(\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,S}\right)^{\top}\left(\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,S}\right)}$
$\displaystyle=\frac{\left(\mathbf{S}^{\perp}\alpha_{0}-\breve{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,S}\alpha_{0}+\breve{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,S}\alpha_{0}-\breve{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,Y}+\breve{\mathbf{W}}^{\perp}_{-1}\theta_{0,-1}+\mathbf{R}^{\perp}+{\bm{\epsilon}}^{\perp}\right)^{\top}\left(\mathbf{S}^{\perp}-\tilde{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,S}\right)}{\left(\mathbf{S}^{\perp}-\tilde{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,S}\right)^{\top}\left(\mathbf{S}^{\perp}-\tilde{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,S}\right)}$
$\displaystyle=\alpha_{0}+\frac{\left(\theta_{0,-1}+\hat{\theta}_{-1,S}\alpha_{0}-\hat{\theta}_{-1,Y}\right)^{\top}\breve{\mathbf{W}}_{-1}^{\perp^{\top}}\left(\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,S}\right)}{\left(\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,S}\right)^{\top}\left(\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,S}\right)}$
$\displaystyle\qquad\qquad+\frac{{\bm{\epsilon}}^{\perp^{\top}}\left(\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,S}\right)}{\left(\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,S}\right)^{\top}\left(\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,S}\right)}+\frac{\tilde{\mathbf{R}}^{\perp^{\top}}\left(\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,S}\right)}{\left(\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,S}\right)^{\top}\left(\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,S}\right)}$

This implies we have:

$\displaystyle\sqrt{n}\left(\hat{\alpha}-\alpha_{0}\right)$
$\displaystyle=\frac{\frac{1}{\sqrt{n}}\left(\theta_{0,-1}+\hat{\theta}_{-1,S}\alpha_{0}-\hat{\theta}_{-1,Y}\right)^{\top}\breve{\mathbf{W}}_{-1}^{\perp^{\top}}\left(\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,S}\right)}{\frac{1}{n}\left(\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,S}\right)^{\top}\left(\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,S}\right)}$
$\displaystyle\qquad\qquad+\frac{\frac{1}{\sqrt{n}}\tilde{\mathbf{R}}^{\perp^{\top}}\left(\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,S}\right)}{\frac{1}{n}\left(\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,S}\right)^{\top}\left(\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,S}\right)}$
$\displaystyle\qquad\qquad\qquad\qquad+\frac{\frac{1}{\sqrt{n}}{\bm{\epsilon}}^{\perp^{\top}}\left(\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,S}\right)}{\frac{1}{n}\left(\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,S}\right)^{\top}\left(\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,S}\right)}$
(4.9)

Using the consistency of the lasso estimates, we first establish the stability
of the common denominator of equation (4.9). More specifically we show in
Proposition A.6 that:

$\frac{1}{n}\left\|\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}\hat{\theta}_{-1,S}\right\|^{2}=\mathbb{E}\left[\left(\tilde{S}-\tilde{W}_{-1}^{\top}\theta^{*}_{S}\right)^{2}\mathds{1}_{|\eta|\leq\tau}\right]+o_{p}(1)=\frac{1}{\left(\Sigma_{\tau}^{-1}\right)_{1,1}}+o_{p}(1)\,.$

By Assumption 4.6, we conclude that the common denominator of equation (4.9)
is $O_{p}(1)$. The next step is to show that the numerators of the first two
terms of the RHS of equation (4.9) is $o_{p}(1)$. The basic intuition behind
this asymptotic negligibility is that the numerator of the first term largely
is basically product of the prediction error of the two LASSO regressions, and
consequently $o_{p}(1)$, even after scaling by $\sqrt{n}$ under certain
condition on sparsity required for debiased LASSO. The numerator of the second
term of the RHS of equation (4.9) is the inner product the residuals of
original model equation (4.4) and the lasso residuals of regression
$\mathbf{S}^{\perp}$ on $\breve{\mathbf{W}}_{-1}^{\perp}$. As we have the
already established the lasso residuals stabilizes, the asymptotic
negligibility of this numerator primarily stems from the asymptotic
negligibility of the residual vector. Details can be found in the Appendix.
The final step is to show that the numerator of the third term of equation
(4.9) is asymptotically normal which follows from an application of
Lindeberg’s central limit theorem. This completes the roadmap of the proof.

##  5 Real data analysis

In this section we illustrate our method by analyzing two real datasets. We
divide our analysis into two subsections, one for each dataset. We first
present a brief description of the data, then present our analysis and compare
our results with the existing one.

###  5.1 Effect of Islamic party on women’s education in Turkey

In this subsection we study the effect of Islamic party rule in Turkey on
women empowerment in terms of their high school education. In the 1994
municipality elections, Islamic parties won several municipal mayor seats in
Turkey. We are interested in investigating whether this winning had any effect
on the education of women, i.e. to determine, statistically, whether the
concern that Islamic control may be inimical towards gender equality is
supported by the data. The dataset we analyze here was collected by Turkish
Statistical Institute and was first analyzed in [29]. Since then it has been
used by several authors, appearing for example, as one of the core
illustrations in [13]444We have downloaded the dataset from
https://github.com/rdpackages-
replication/CIT_2019_CUP/blob/master/CIT_2019_Cambridge_polecon.csv.. The
dataset consists of $n=2629$ rows where the rows represent municipalities, the
units of our analysis. The main target/response variable $Y$ is the percentage
of women in the 15-20 year age-group who were recorded to have completed their
high school education in the 2000 census. As mentioned in [29], this is the
group that should have been most affected by the decisions made by the winners
of the 1994 election. The treatment determining variable $Q$ is the difference
in vote share between the largest Islamic party (i.e. the Islamic party which
got maximum votes among all Islamic parties) and the largest secular party
(i.e. the non-Islamic party which got the maximum votes among all non-Islamic
parties). Hence, the cutoff is $0$: i.e. if $Q_{i}>0$, then the $i^{th}$
municipality unit elected an Islamic party and if $Q_{i}<0$, a secular party.
The description of the available co-variates is presented in Table 4 of the
supplementary document. For $X$ and $Z$, we use all the co-variates presented
in that table except i89 because almost $1/3$’rd of the observations (729
many) are missing for this variable. We estimate $\alpha_{0}$ by
$\overline{\hat{\alpha}}$ as described in Section 2. To test:

$H_{0}:\alpha_{0}=0\ \ \ \ vs\ \ \ \ H_{1}:\alpha_{0}\neq 0$

we construct an Efron-bootstrap based confidence interval over 500 iterations.
We present our finding in Table 1. From the table, it is clear that, we don’t
have enough evidence to reject $H_{0}$ at $5\%$ level as the confidence
interval contains $0$. Hence we conclude that, there is no significant effect
of Islamic ruling party on women’s education.

We next compare our result to that of [29] and [13]. [29] implemented a
simpler model for RDD:

$Y_{i}=\beta_{0}+\alpha_{0}S_{i}+f(Q_{i})+{\epsilon}_{i}\,.$

with $f$ being a polynomial function and only those observations were used
where $Q_{i}\in(-h,h)$ for some optimal choice of the bandwidth $h$ (chosen
according to the prescription of [21]). The authors found that Islamic party
rule has a significantly positive effect on women’s education at $1\%$ level
test. On the other hand, [13] implemented the model based on (1.2). We
replicate their result using the R-package _rdrobust_ as advocated in [13].
The function rdrobust inside the package _rdrobust_ , takes input $Y,Q,X$ and
splits out three different types of estimate of $\alpha_{0}$ (along with
$95\%$ confidence interval), namely conventional, bias-corrected and Robust.
As mentioned in [9], conventional presents when the conventional RD estimates
(i.e. solving equation (1.2) via local polynomial regression) with
conventional standard errors, bias-corrected implies bias-corrected RD
estimate with conventional standard errors and robust indicates bias-corrected
RD estimates with robust standard errors (see [11]). We use the parameters
kernel = ’triangular’, scaleregul = 1, p = 1, bwselect = ’mserd’ of rdrobust
function to run the analysis. As evident from Table 2, all these estimates
reject $H_{0}$ at $5\%$ level stipulating a strictly positive effect of
Islamic party on the education of women, while our method fails to reject the
null at the same level.

Table 1: Summary Statistics of data on Islamic party based on our method

Point Estimate 0.4071513 Bootstrap mean. 0.5760144 Bootstrap s.e. 0.48115
Bootstrap $95\%$ C.I. $(-0.4234894,1.41942)$

Table 2: Summary Statistics of data on Islamic party based on [13]

Name of methods Coeff CI Lower CI Upper Conventional 3.005951 0.9622239
5.049678 Bias-Corrected 3.204837 1.1611103 5.248564 Robust 3.204837 0.8266720
5.583003

###  5.2 Effect of probation on subsequent GPA

We next analyze an educational dataset, originally collected and analyzed in
[25]555We have collected the dataset from
https://www.openicpsr.org/openicpsr/project/113751/version/V1/view;jsessionid=A6C09FD5CD7DB8E18EAA77B75BD893B2.where
we investigate whether putting students on academic probation due to grades
below a pre-determined cutoff has any effect on their subsequent GPA. The data
are based on students from 3 independent campuses of a large Canadian
university – a major campus and other satellite campuses. The acceptance rate
in the major campus is around $55\%$ and in the satellite campuses around
$77\%$. The data were collected over $8$ cohorts of students till the end of
the 2005 academic year. To observe the students for at least two years, only
those who entered the university prior to the beginning of the 2004 academic
year were considered. After being put on academic probation in their first
year, some students left the university. We, therefore, only have access to
GPA during the second year for those students who stayed. Thus, our $Y$
variable is the GPA of the first academic term in the second year and the
treatment $S$ is whether the student was put on probation. The treatment
determining variable $Q$ is the difference between the first year GPA and the
cutoff for academic probation: if $Q<0$, the student is put on academic
probation, otherwise not. The covariates we consider here ($X=Z$) are
presented in Table 5 of the supplementary document. In Table 3 we summarize
our result. It is immediate from the bootstrap confidence interval from Table
3 that for testing $H_{0}:\alpha=0$ vs $H_{1}:\alpha\neq 0$, we have enough
evidence to Reject $H_{0}$ at the $5\%$ level and conclude that the students
who are put on academic probation and continue with their education, tend to
improve their performance (note that the estimated $\alpha_{0}$ is positive)
in the subsequent academic year. This makes sense, as the students who did not
leave university after being put on academic probation must have a strong
incentive to work harder so that they are not expelled from the university.
Our findings are in harmony with those obtained in [25], where the author also
found the treatment effect to be significant at the $5\%$ level.

######  Remark 5.1.

Note that in Table 1 and Table 3, we presented bootstrap confidence intervals
for the treatment effect instead of asymptotic confidence intervals. This is
because consistent estimation of the the asymptotic variance of our estimator
is not straightforward. Recall that from Theorem 3.6, the asymptotic variance
of our estimator is
$e_{1}^{\top}\Omega_{\tau}^{-1}\Omega_{\tau}^{*}\Omega_{\tau}^{-1}e_{1}$. As
mentioned in Step 1 of the sketch of the proof of Theorem 3.6,
$\mathbf{W}^{\top}{\sf proj}^{\perp}_{\tilde{\mathbf{N}}_{K,a}}\mathbf{W}/n$
is a consistent estimator of $\Omega_{\tau}$ but it is hard to estimate
$\Omega^{*}_{\tau}$ consistently, which forces us to resort to the bootstrap
confidence interval.

Table 3: Summary Statistics of data on GPA data

Point Estimate 0.2733371 Bootstrap mean. 0.2817404 Bootstrap s.e. 0.016672
Bootstrap $95\%$ C.I. $(0.248588,0.3158934)$

##  6 Conclusion and possible extensions

In this paper, we proposed a new approach to estimate an non-randomized
treatment effect at the $\sqrt{n}$ rate and showed that under homoscedastic
normal errors our method is semiparametrically efficient. We also pointed out
in Remark 3.8 that one can use weighted least squares instead of ordinary
least squares to take care of heterogenous errors. However, the normality
assumption is necessary for semiparametric efficiency as we use least squares
for estimating the treatment effect. We now discuss some natural extensions of
our models that are worth analyzing as potential future research problems.

###  6.1 Non-constant treatment effect

Consider the following extension of our model with non-constant treatment
effect:

$\displaystyle Y_{i}$
$\displaystyle=\alpha(\eta_{i})\mathds{1}_{Q_{i}>0}+X_{i}^{\top}\beta_{0}+\nu_{i}=\alpha(\eta_{i})\mathds{1}_{Q_{i}>0}+X_{i}^{\top}\beta_{0}+b(\eta_{i})+{\epsilon}_{i}$
(6.1) $\displaystyle Q_{i}$ $\displaystyle=Z_{i}^{\top}\gamma_{0}+\eta_{i}\,.$
(6.2)

where as before $b(\eta_{i})=\mathbb{E}[\nu\mid\eta_{i}]$. This generalization
assumes that the response of a treated candidate with higher abilities is
boosted in comparison to another treated candidate with lower ability. As an
example, a more capable student upon entering into a prestigious grad school,
will most likely get a better mentor, resulting in an amplification of their
academic prowess. The random variable $\alpha(\eta)$ represents the
_conditional treatment effect_ as:

$\mathbb{E}[Y\mid X,\eta,Q>0]-\mathbb{E}[Y\mid X,\eta,Q<0]=\alpha(\eta)\,.$

i.e. conditioning on $(X,\eta)$ (background information and innate ability),
$\alpha(\eta)$ quantifies the difference between the responses of treated and
untreated samples. It can be shown that the estimator proposed in our
manuscript based on the model with constant $\alpha_{0}$, estimates, in the
newly proposed model, a _weighted average_ of $\alpha(\eta)$, i.e.
$\mathbb{E}[p(\eta)\alpha(\eta)]/\mathbb{E}[p(\eta)]$ where $p(\cdot)$ is a
non-negative integrable function depending on other parameters. This weight
function is basically a function of the information for $\alpha$ in each
$\eta$ \- stratum. So our current method is applicable when the parameter of
interest is a _weighted average treatment effect_.

However if the parameter of interest is _unweighted average treatment effect_
$\mathbb{E}[\alpha(\eta)]$, then our method can be slightly modified as
follows to yield a $\sqrt{n}$-consistent estimator:

  1. 1.

Split the data into three (almost) equal parts say
$\mathcal{D}_{1},\mathcal{D}_{2},\mathcal{D}_{3}$.

  2. 2.

From $\mathcal{D}_{1}$, impute $\hat{\eta}$ from equation (6.2) by regressing
$Q$ on $X$.

  3. 3.

Estimate $\alpha(\cdot)$ and $b(\cdot)$ (and their derivatives) from
$\mathcal{D}_{2}$ using a B-spline series expansion.

  4. 4.

Note that equation (6.1) can be expanded as:

$\displaystyle Y_{i}$
$\displaystyle=\alpha(\eta_{i})\mathds{1}_{Q_{i}>0}+X_{i}^{\top}\beta_{0}+b(\eta_{i})+{\epsilon}_{i}$
$\displaystyle=\alpha(\hat{\eta}_{i})\mathds{1}_{Q_{i}>0}+b(\hat{\eta}_{i})+X_{i}^{\top}\beta_{0}+(\hat{\gamma}-\gamma_{0})^{\top}Z_{i}\left(\alpha^{\prime}(\hat{\eta}_{i})\mathds{1}_{Q_{i}>0}+b^{\prime}(\hat{\eta}_{i})\right)+{\epsilon}_{i}+R_{i}$
$\displaystyle\approx\alpha(\hat{\eta}_{i})\mathds{1}_{Q_{i}>0}+b(\hat{\eta}_{i})+X_{i}^{\top}\beta_{0}+(\hat{\gamma}-\gamma_{0})^{\top}Z_{i}\left(\hat{\alpha}^{\prime}(\hat{\eta}_{i})\mathds{1}_{Q_{i}>0}+\hat{b}^{\prime}(\hat{\eta}_{i})\right)+{\epsilon}_{i}+R_{i}$

where in the last line we replaced $\alpha^{\prime}$ and $b^{\prime}$ by their
estimates obtained from $\mathcal{D}_{2}$ in the previous step. Finally, we
can re-estimate $\alpha$ from the above equation via a non-parametric method
(i.e. expanding via a B-spline basis and regressing $Y$ on an appropriate set
of covariates). Our final estimator becomes:

$\hat{\mathbb{E}}[\alpha(\eta)]=\frac{1}{n}\sum_{i=1}^{n}\hat{\alpha}(\hat{\eta}_{i})\,.$

An analysis similar to that in our paper indicates that this estimator has
$\sqrt{n}$ rate of convergence and is asymptotically normal. However whether
this is semi-parametrically efficient is not clear at this moment and a
potentially interesting problem for future research.

###  6.2 Bootstrap consistency

As noted in Remark 5.1, we use a bootstrap confidence interval instead of the
asymptotic one owing to the intricate form of the asymptotic variance of our
estimator, which makes it hard to estimate from the data. Therefore, an
immediate question of interest is to investigate whether the bootstrap is
consistent under our model assumptions. Although empirical evidence suggests
that this is the case, a rigorous theoretical undertaking is essential to
establish the claim.

##  Appendix A Proof of Theorem 4.10

###### Proof.

For the ease of notation, we assume $\dim(X)=p_{1}\asymp\dim(Z)=p_{2}\asymp
p$. One can extend our proof quite easily for general $p_{1},p_{2}$ (i.e. when
they are not of same order) with a careful booking for the dimension factor.
The entire proof is quite long and tedious, therefore in this appendix we will
state the key steps and provide proofs of the main parts. Proofs of all
supplementary lemmas and propositions can be found in Appendix.

Recall that our estimation procedure consists of three parts:

  * •

Estimate $\gamma_{0}$ from $\mathcal{D}_{1}$.

  * •

Estimate $b^{\prime}$ from $\mathcal{D}_{2}$.

  * •

Estimate $\alpha_{0}$ from $\mathcal{D}_{3}$.

###  A.1 Estimation of $\hat{\gamma}_{0}$ from $\mathcal{D}_{1}$

This part is the easiest among all the three parts. We estimate $\gamma_{0}$
by doing LASSO of $\mathbf{Y}$ in $\mathbf{Z}$. Note that, by sub-gaussianity
assumption, $\mathbf{Z}$ satisfies restricted eigenvalue condition with high
probability. Therefore, by standard lasso calculation with the tuning
parameter $\lambda\asymp\sqrt{\log{p}/n}$ we have:

$\displaystyle\left\|\hat{\gamma}-\gamma_{0}\right\|_{2}^{2}\lesssim_{\mathbb{P}}\frac{s_{\gamma}\log{p}}{n}\,,$
$\displaystyle\left\|\hat{\gamma}-\gamma_{0}\right\|_{1}\lesssim_{\mathbb{P}}s_{\gamma}\sqrt{\frac{\log{p}}{n}}\,.$

We will use this estimator in the subsequent analysis.

###  A.2 Estimation of $b^{\prime}$ from $\mathcal{D}_{2}$

In this subsection, we present the estimation error for $b^{\prime}$. An in
our fixed dimensional analysis, we expand the model equation as:

$\mathbf{Y}=\mathbf{X}\beta_{0}+\mathbf{N}_{k}(\hat{\eta})\omega_{b}+\mathbf{R}_{1}+\mathbf{R}_{2}+{\bm{\epsilon}}\,.$

where $\mathbf{N}_{K}(\hat{\eta})\in{\mathbb{R}}^{n\times K}$ with
$\mathbf{N}_{K}(\hat{\eta})_{i,j}=\tilde{N}_{K,j}(\hat{\eta}_{i})$, $N_{K,j}$
being the scaled $K^{th}$ B-spline basis (see Section F of the supplementary
document for details). Here $\omega_{b}$ is the coefficient of _best B-spline
approximator_ of $b$ with respect to $K$ basis. Our aim is to estimate
$\omega_{b}$ from the data, as upon obtaining $\hat{\omega}_{b}$ one can
define $\hat{b}^{\prime}(t)=\nabla N_{K}(t)^{\top}\hat{\omega}_{b}$, where
$\nabla\tilde{N}_{K}(t)$ is the vector of derivatives of B-spline basis.
Therefore the estimation error $\hat{b}^{\prime}$ can be bounded as:

$\sup_{|t|\leq\tau}\left|b^{\prime}(t)-\hat{b}^{\prime}(t)\right|\leq\sup_{|t|\leq\tau}\left\|\nabla\tilde{N}_{k}(t)\right\|\times\left\|\hat{\omega}_{b}-\omega_{b}\right\|+\sup_{|t|\leq\tau}\left|b^{\prime}(t)-\tilde{N}_{K}(t)^{\top}\omega_{b}\right|\,.$

The second term is the B-spline approximation error, which is bounded by the
order of $K^{-(\upsilon-1)}$ (see Theorem F.1 of the supplementary document).
For the first term or estimation error, as mentioned in Section F, we have
$\sup_{|t|\leq\tau}\|\nabla\tilde{N}_{k}(t)\|\lesssim K^{3/2}$. This further
implies:

$\sup_{|t|\leq\tau}\left|b^{\prime}(t)-\hat{b}^{\prime}(t)\right|\lesssim
K^{3/2}\left\|\hat{\omega}_{b}-\omega_{b}\right\|+K^{-(\upsilon-1)}\,,$

and consequently, the estimation error $\hat{b}^{\prime}$ completely depends
on the estimation error of $\hat{\omega}_{b}$. In the proposition we present a
bound on the estimation error on $\hat{\omega}_{b}$ (and consequently on
$\hat{b}^{\prime}$):

######  Proposition A.1.

Under the Assumptions stated in Section 4, we have:

$\left\|\hat{\omega}_{b}-\omega_{b}\right\|\lesssim_{\mathbb{P}}\sqrt{K}s_{\beta}\left(\sqrt{\frac{s_{\gamma}\log{p}}{n}}+K^{-2\upsilon}\right)+K^{-\upsilon}+\sqrt{\frac{s_{\gamma}\log{p}}{n}}+\sqrt{\frac{K}{n}}$

and consequently:

$\sup_{|t|\leq\tau}\left|b^{\prime}(t)-\hat{b}^{\prime}(t)\right|\lesssim_{\mathbb{P}}\left(\frac{n}{s^{2}_{\beta}s_{\gamma}\log{p}}\right)^{\frac{-(\upsilon-\frac{3}{2})}{2\upsilon+1}}\,.$

###### Proof of Proposition A.1.

The method is, first we estimate $\beta_{0}$ using LASSO as follows:

$\hat{\beta}=\mathop{\rm
argmin}_{\beta}\left\\{\frac{1}{2n}\left\|P_{\mathbf{N}_{k}}^{\perp}\left(\mathbf{Y}-\mathbf{X}\beta\right)\right\|^{2}+\lambda\|\beta\|_{1}\right\\}$

Next we estimate $\omega_{b}$ is:

$\displaystyle\hat{\omega}_{b}=\left(\frac{\mathbf{N}_{k}^{\top}\mathbf{N}_{k}}{n}\right)^{-1}\frac{\mathbf{N}_{k}^{\top}\left(\mathbf{Y}-\mathbf{X}\hat{\beta}\right)}{n}$
(A.1)

Therefore, to establish a bound on $\|\hat{\omega}_{b}-\omega_{b}\|$, a bound
on the estimation error $\hat{\beta}-\beta_{0}$ is necessary, which is
established in the the following lemma:

######  Lemma A.2.

Under our assumptions, we have the the following bound on the estimation error
of $\beta_{0}$:

$\left\|\hat{\beta}-\beta_{0}\right\|_{2}^{2}\lesssim_{\mathbb{P}}s_{\beta}\left(\sqrt{\frac{s_{\gamma}\log{p}}{n}}+K^{-2\upsilon}\right)^{2}\lesssim_{\mathbb{P}}\frac{s_{\beta}s_{\gamma}\log{p}}{n}+K^{-4\upsilon}\triangleq
r_{n,\beta_{0}}^{2}\,,$

where $K$ is the number of B-spline basis used for extending $b$.

###### Proof of Lemma A.2.

By the basic LASSO inequality, we have:

$\displaystyle\frac{1}{2n}\left\|P_{\mathbf{N}_{k}}^{\perp}\left(\mathbf{Y}-\mathbf{X}\hat{\beta}\right)\right\|^{2}+\lambda\|\hat{\beta}\|_{1}$
$\displaystyle\leq\frac{1}{2n}\left\|P_{\mathbf{N}_{k}}^{\perp}\left(\mathbf{Y}-\mathbf{X}\beta_{0}\right)\right\|^{2}+\lambda\|\beta_{0}\|_{1}$

Some algebraic manipulation yields:

$\displaystyle\frac{1}{2n}\left\|P_{\mathbf{N}_{k}}^{\perp}\mathbf{X}\left(\hat{\beta}-\beta_{0}\right)\right\|^{2}+\lambda\|\hat{\beta}\|_{1}$
$\displaystyle\leq\frac{1}{n}\left\|\left({\bm{\epsilon}}+\mathbf{R}\right)^{\top}P_{\mathbf{N}_{k}}^{\perp}\mathbf{X}\right\|_{\infty}\left\|\hat{\beta}-\beta_{0}\right\|_{1}+\lambda\|\beta_{0}\|_{1}$

The matrix $\mathbb{P}_{\mathbf{N}_{K}}^{\perp}\mathbf{X}$ satisfies RE
condition with high probability due to Assumption 4.5 and Proposition A.5. For
the optimal value of $\lambda$, we need a bound on
$(1/n)\|\left({\bm{\epsilon}}+\mathbf{R}\right)^{\top}P_{\mathbf{N}_{k}}^{\perp}\mathbf{X}\|_{\infty}$,
for which we bound
$(1/n)\|{\bm{\epsilon}}^{\top}P_{\mathbf{N}_{k}}^{\perp}\mathbf{X}\|_{\infty}$
and $(1/n)\|\mathbf{R}^{\top}P_{\mathbf{N}_{k}}^{\perp}\mathbf{X}\|_{\infty}$
separately. To bound the term with ${\bm{\epsilon}}$, we use the sub-gaussian
concentration inequality:

$\displaystyle\mathbb{P}\left(\frac{1}{n}\left|{\bm{\epsilon}}^{\top}P_{\mathbf{N}_{k}}^{\perp}\mathbf{X}_{*,j}\right|>t\mid\sigma(X,Z,\eta,\mathcal{D}_{1})\right)$
$\displaystyle\leq 2{\sf
exp}{\left(-C\frac{nt^{2}}{\frac{1}{n}\mathbf{X}_{*,j}^{\top}P_{\mathbf{N}_{k}}^{\perp}\mathbf{X}_{*,j}}\right)}$
$\displaystyle\leq 2{\sf
exp}{\left(-C\frac{nt^{2}}{\frac{1}{n}\mathbf{X}_{*,j}^{\top}\mathbf{X}_{*,j}}\right)}$

Now from the sub-gaussianity of $X_{j}$ (with sub-gaussianity constant
$\sigma_{W}$) applying Lemma I.2, we have probability going to $1$:

$\max_{1\leq j\leq n}\frac{\|\mathbf{X}_{*,j}\|^{2}}{n}\leq 3\sigma_{W}\,.$

Define the above event as $\Omega_{n}$. Using this we have:

$\displaystyle\mathbb{P}\left(\frac{1}{n}\left|{\bm{\epsilon}}^{\top}P_{\mathbf{N}_{k}}^{\perp}\mathbf{X}_{*,j}\right|>t\right)$
$\displaystyle\leq\mathbb{P}\left(\frac{1}{n}\left|{\bm{\epsilon}}^{\top}P_{\mathbf{N}_{k}}^{\perp}\mathbf{X}_{*,j}\right|>t,\Omega_{n}\right)+o(1)$
$\displaystyle=\mathbb{E}\left[\mathds{1}_{\Omega_{n}}\mathbb{P}\left(\frac{1}{n}\left|{\bm{\epsilon}}^{\top}P_{\mathbf{N}_{k}}^{\perp}\mathbf{X}_{*,j}\right|>t\mid\sigma(X,Z,\eta,\mathcal{D}_{1})\right)\right]+o(1)$
$\displaystyle\leq 2\mathbb{E}\left[\mathds{1}_{\Omega_{n}}{\sf
exp}{\left(-C\frac{nt^{2}}{\frac{1}{n}\mathbf{X}_{*,j}^{\top}\mathbf{X}_{*,j}}\right)}\right]+o(1)$
$\displaystyle\leq 2{\sf
exp}{\left(-C\frac{nt^{2}}{3\sigma_{W}^{2}}\right)}+o(1)$

Therefore an appropriate choice of $t$ yields:

$\displaystyle\frac{1}{n}\left\|{\bm{\epsilon}}^{\top}P_{\mathbf{N}_{K}}^{\top}\mathbf{X}\right\|_{\infty}\lesssim_{\mathbb{P}}\sqrt{\frac{\log{p}}{n}}\,.$
(A.2)

Now for the remainder term, recall that the remainder term $\mathbf{R}$
consists of $\mathbf{R}=\mathbf{R}_{1}+\mathbf{R}_{2}$ where

$\displaystyle\mathbf{R}_{1,i}$ $\displaystyle=b(\eta_{i})-b(\hat{\eta}_{i})$
$\displaystyle\mathbf{R}_{2,i}$
$\displaystyle=b(\hat{\eta}_{i})-\mathbf{N}_{K}(\hat{\eta}_{i})^{\top}\omega_{b}\,.$

We now bound these two remainder terms separately. For the first remainder
term:

$\displaystyle\frac{1}{n}\left\|\mathbf{R}_{1}^{\top}P_{\mathbf{N}_{K}}^{\perp}\mathbf{X}\right\|_{\infty}$
$\displaystyle=\max_{1\leq j\leq
p}\frac{1}{n}\left|\mathbf{R}_{1}^{\top}P_{\mathbf{N}_{K}}^{\perp}\mathbf{X}_{*,j}\right|$
$\displaystyle\leq\sqrt{\frac{1}{n}\mathbf{R}_{1}^{\top}\mathbf{R}}\times\sqrt{\max_{1\leq
j\leq p}\frac{1}{n}\|\mathbf{X}_{*,j}\|^{2}}$ (A.3)

Again it follows from Lemma I.2 that the second term of the above inequality
is $O_{p}(1)$ based on the subgaussianity of $\mathbf{X}_{*,j}$. For the first
term, note that as $b$ is Lipschitz (Assumption 4.1), we have:

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\mathbf{R}_{1,i}^{2}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\left(b(\hat{\eta}_{i})-b(\eta_{i})\right)^{2}$
$\displaystyle\lesssim\frac{1}{n}\sum_{i=1}^{n}\left(\hat{\eta}_{i}-\eta_{i}\right)^{2}\hskip
14.45377pt[\text{As }b\text{ is Lipschitz}]$
$\displaystyle\leq\left\|\hat{\gamma}_{n}-\gamma_{0}\right\|^{2}\times\frac{1}{n}\sum_{i=1}^{n}\left(Z_{i}^{\top}a_{n}\right)^{2}$
$\displaystyle\lesssim_{\mathbb{P}}\frac{s_{\gamma}\log{p}}{n}\,.$

where the last line follows from the fact that
$\|\hat{\gamma}_{n}-\gamma_{0}\|^{2}\lesssim_{\mathbb{P}}(s_{\gamma}\log{p})/n$
(from LASSO on $\mathcal{D}_{1}$) and
$(1/n)\sum_{i=1}^{n}\left(Z_{i}^{\top}a_{n}\right)^{2}\lesssim_{\mathbb{P}}1$
follows from subgaussianity of $Z_{i}$ along with Lemma I.2. Using this bound
in equation (A.3) we conclude:

$\displaystyle\frac{1}{n}\left\|\mathbf{R}_{1}^{\top}P_{\mathbf{N}_{K}}^{\perp}\mathbf{X}\right\|_{\infty}$
$\displaystyle\lesssim_{\mathbb{P}}\sqrt{\frac{s_{\gamma}\log{p}}{n}}\,.$
(A.4)

For the other remainder term, we expand $\mathbf{X}_{*,j}$ as:

$\mathbf{X}_{*,j}=\mathbf{X}_{*,j}-\check{m}_{j}(\hat{\mathbf{\eta}})+\check{m}_{j}(\hat{\mathbf{\eta}})=\mathbf{X}_{*,j}-\check{m}_{j}(\hat{\mathbf{\eta}})+\mathbf{N}_{K}(\hat{\eta})w_{j}+\mathbf{R}_{j}$

This implies:

$\mathbf{R}_{2}^{\top}P_{\mathbf{N}_{K}}^{\perp}\mathbf{X}_{*,j}=\mathbf{R}_{2}^{\top}P_{\mathbf{N}_{K}}^{\perp}\left(\mathbf{X}_{*,j}-\check{m}_{j}(\hat{\mathbf{\eta}})\right)+\mathbf{R}_{2}^{\top}P_{\mathbf{N}_{K}}^{\perp}\mathbf{R}_{j}\,.$

As both $\mathbf{R}_{2}$ and $P_{\mathbf{N}_{K}}^{\perp}$ are function of
$\hat{\eta}$, we have from the sub-gaussian concentration bound:

$\displaystyle\mathbb{P}\left(\frac{1}{n}\left|\mathbf{R}_{2}^{\top}P_{\mathbf{N}_{K}}^{\perp}\left(\mathbf{X}_{*,j}-\check{m}_{j}(\hat{\mathbf{\eta}})\right)\right|>t\mid\hat{\eta}\right)$
$\displaystyle\leq 2{\sf
exp}{\left(-C\frac{nt^{2}}{\frac{1}{n}\mathbf{R}_{2}^{\top}P_{\mathbf{N}_{K}}^{\perp}\mathbf{R}_{2}}\right)}$
$\displaystyle\leq 2{\sf
exp}{\left(-C\frac{nt^{2}}{\frac{1}{n}\mathbf{R}_{2}^{\top}\mathbf{R}_{2}}\right)}\,.$

From the spline approximation error,
$(1/n)\mathbf{R}_{2}^{\top}\mathbf{R}_{2}\leq K^{-2\upsilon}$ (where $K$ is
the number of basis and $\upsilon$ is the smoothness index of $b$ (see
Assumption 4.1). Using this we have:

$\mathbb{P}\left(\frac{1}{n}\left|\mathbf{R}_{2}^{\top}P_{\mathbf{N}_{K}}^{\perp}\left(\mathbf{X}_{*,j}-\check{m}_{j}(\hat{\mathbf{\eta}})\right)\right|>t\mid\hat{\eta}\right)\leq
2{\sf exp}{\left(-C\frac{nt^{2}}{K^{-2\upsilon}}\right)}$

which, upon unconditioning, taking union bound and choosing a suitable value
of $t$ yields:

$\max_{1\leq j\leq
p}\frac{1}{n}\left|\mathbf{R}_{2}^{\top}P_{\mathbf{N}_{K}}^{\perp}\left(\mathbf{X}_{*,j}-\check{m}_{j}(\hat{\mathbf{\eta}})\right)\right|\lesssim_{\mathbb{P}}K^{-\upsilon}\sqrt{\frac{\log{p}}{n}}\,.$

For the other remainder term, i.e.
$(1/n)\mathbf{R}_{2}^{\top}P_{\mathbf{N}_{K}}^{\perp}\mathbf{R}_{j}$, an
application of Cauchy-Schwarz inequality yields:

$\frac{1}{n}\left|\mathbf{R}_{2}^{\top}P_{\mathbf{N}_{K}}^{\perp}\mathbf{R}_{j}\right|\leq\sqrt{\frac{1}{n}\mathbf{R}_{2}^{\top}\mathbf{R}_{2}}\times\sqrt{\frac{1}{n}\mathbf{R}_{j}^{\top}\mathbf{R}_{j}}\leq
K^{-2\upsilon}\,.$

Therefore we have:

$\displaystyle\frac{1}{n}\left\|\mathbf{R}_{2}^{\top}P_{\mathbf{N}_{K}}^{\perp}\mathbf{X}\right\|_{\infty}\lesssim_{\mathbb{P}}K^{-\upsilon}\sqrt{\frac{\log{p}}{n}}+K^{-2\upsilon}\,.$
(A.5)

Combining the bounds in equation (A.2), (A.5) and (A.5) we have:

$\displaystyle\frac{1}{n}\left\|\left({\bm{\epsilon}}+\mathbf{R}\right)^{\top}P_{\mathbf{N}_{K}}^{\perp}\mathbf{X}\right\|_{\infty}\lesssim_{\mathbb{P}}\sqrt{\frac{s_{\gamma}\log{p}}{n}}+K^{-2\upsilon}\asymp\lambda\,.$
(A.6)

With this choice of $\lambda$, standard LASSO analysis completes the proof. ∎

Going back to the definition of $\hat{\omega}_{b}$ (equation (A.1)), we expand
the its estimation error as follows:

$\displaystyle\hat{\omega}_{b}-\omega_{b}$
$\displaystyle=\underbrace{\left(\frac{\mathbf{N}_{k}^{\top}\mathbf{N}_{k}}{n}\right)^{-1}\frac{\mathbf{N}_{k}^{\top}\mathbf{X}\left(\hat{\beta}-\beta_{0}\right)}{n}}_{T_{1}}+\underbrace{\left(\frac{\mathbf{N}_{k}^{\top}\mathbf{N}_{k}}{n}\right)^{-1}\frac{\mathbf{N}_{k}^{\top}(\mathbf{R}_{1}+\mathbf{R}_{2})}{n}}_{T_{2}}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad+\underbrace{\left(\frac{\mathbf{N}_{k}^{\top}\mathbf{N}_{k}}{n}\right)^{-1}\frac{\mathbf{N}_{k}^{\top}{\bm{\epsilon}}}{n}}_{T_{3}}$

We next bound each $T_{i}$ separately. Same argument as for the fixed
dimensional analysis (See the proof of Proposition 2.1) establishes:

$\left\|\left(\frac{\mathbf{N}_{k}^{\top}\mathbf{N}_{k}}{n}\right)^{-1}\right\|_{op}\lesssim_{\mathbb{P}}1\,.$

Therefore we can bound $T_{1}$ as:

$\displaystyle
T_{1}\lesssim_{\mathbb{P}}\left\|\frac{\mathbf{N}_{k}^{\top}\mathbf{X}\left(\hat{\beta}-\beta_{0}\right)}{n}\right\|_{2}$
$\displaystyle\leq\frac{\sqrt{K}}{n}\max_{1\leq j\leq
k}\left|\mathbf{N}_{K,*j}^{\top}\mathbf{X}\left(\hat{\beta}-\beta_{0}\right)\right|$
$\displaystyle\leq\frac{\sqrt{K}}{n}\left\|\hat{\beta}-\beta_{0}\right\|_{1}\times\max_{1\leq
j\leq k}\max_{1\leq l\leq
p}\left|\mathbf{N}_{K,*j}^{\top}\mathbf{X}_{*,l}\right|$
$\displaystyle\lesssim_{\mathbb{P}}\sqrt{K}\left\|\hat{\beta}-\beta_{0}\right\|_{1}\lesssim_{\mathbb{P}}\sqrt{K}s_{\beta}\left(\sqrt{\frac{s_{\gamma}\log{p}}{n}}+K^{-2\upsilon}\right)\,.$

For $T_{2}$ and $T_{3}$, the term containing residuals:

$\displaystyle
T_{2}\leq\left(\frac{\mathbf{N}_{k}^{\top}\mathbf{N}_{k}}{n}\right)^{-1}\frac{\mathbf{N}_{k}^{\top}(\mathbf{R}_{1}+\mathbf{R}_{2})}{n}$
$\displaystyle\leq\left\|\left(\frac{\mathbf{N}_{k}^{\top}\mathbf{N}_{k}}{n}\right)^{-1/2}\right\|_{op}\left\|\frac{\mathbf{R}_{1}+\mathbf{R}_{2}}{\sqrt{n}}\right\|_{2}$
$\displaystyle\lesssim_{\mathbb{P}}\left\|\frac{\mathbf{R}_{1}+\mathbf{R}_{2}}{\sqrt{n}}\right\|_{2}$
$\displaystyle\leq\sqrt{\frac{1}{n}\mathbf{R}_{1}^{\top}\mathbf{R}_{1}}+\sqrt{\frac{1}{n}\mathbf{R}_{2}^{\top}\mathbf{R}_{2}}$
$\displaystyle\lesssim_{\mathbb{P}}\left(K^{-\upsilon}+\sqrt{\frac{s_{\gamma}\log{p}}{n}}\right)\,.$

where the last rate inequality follows from the bounds on $\mathbf{R}_{1}$ and
$\mathbf{R}_{2}$ established in the proof of Lemma A.2. For the error term we
use the sub-gaussian bound using the fact that $\mathbb{E}[{\epsilon}\mid
X,Z,\eta]=0$. For this, note that for any vector $z\in\mathbb{R}^{K}$ we have:

$\|z\|=\sup_{\|v\|=1}|z^{\top}v|\,.$

Define $\mathcal{N}_{1/2}(S^{K-1}$ to be $1/2$-covering number of the sphere
in dimension $K$. Then we know $\mathcal{N}_{1/2}(S^{K-1}\leq 5^{K}$. Also we
have:

$\|z\|\leq\sup_{v_{1}\in\mathcal{N}_{1/2}(S^{K-1})}|z^{\top}v_{1}|\leq\frac{1}{2}\|z\|$

which implies:

$\|z\|\leq 2\sup_{v_{1}\in\mathcal{N}_{1/2}(S^{K-1})}|z^{\top}v_{1}|\,.$

Using this we have:

$\displaystyle\mathbb{P}\left(\sup_{v_{1}\in\mathcal{N}_{1/2}(S^{K-1})}\left|v_{1}^{\top}\left(\mathbf{N}_{K}^{\top}\mathbf{N}_{K}\right)^{-1}\mathbf{N}_{K}^{\top}{\bm{\epsilon}}\right|>t\right)$
$\displaystyle=\mathbb{E}_{\hat{\eta}}\left[\mathbb{P}\left(\sup_{v_{1}\in\mathcal{N}_{1/2}(S^{K-1})}\left|v_{1}^{\top}\left(\mathbf{N}_{K}^{\top}\mathbf{N}_{K}\right)^{-1}\mathbf{N}_{K}^{\top}{\bm{\epsilon}}\right|>t\mid\hat{\eta}\right)\right]$
$\displaystyle\leq\mathbb{E}_{\hat{\eta}}\left[\sum_{v_{1}\in\mathcal{N}_{1/2}(S^{K-1})}\mathbb{P}\left(\left|v_{1}^{\top}\left(\mathbf{N}_{K}^{\top}\mathbf{N}_{K}\right)^{-1}\mathbf{N}_{K}^{\top}{\bm{\epsilon}}\right|>t\mid\hat{\eta}\right)\right]$
$\displaystyle\leq\mathbb{E}_{\hat{\eta}}\left[\sum_{v_{1}\in\mathcal{N}_{1/2}(S^{K-1})}2{\sf
exp}{\left(-c\frac{nt^{2}}{v_{1}^{\top}\left(\frac{\mathbf{N}_{k}^{\top}\mathbf{N}_{k}}{n}\right)^{-1}v_{1}}\right)}\right]$
$\displaystyle\leq\mathbb{E}_{\hat{\eta}}\left[2{\sf
exp}{\left(K\log{5}-c\frac{nt^{2}}{\left\|\left(\frac{\mathbf{N}_{k}^{\top}\mathbf{N}_{k}}{n}\right)^{-1}\right\|_{op}}\right)}\right]$

This implies:

$T_{3}=\left(\frac{\mathbf{N}_{k}^{\top}\mathbf{N}_{k}}{n}\right)^{-1}\frac{\mathbf{N}_{k}^{\top}{\bm{\epsilon}}}{n}\lesssim_{\mathbb{P}}\sqrt{\frac{K}{n}}\,.$

Combining the bounds on $T_{1},T_{2},T_{3}$ we have:

$\left\|\hat{\omega}_{b}-\omega_{b}\right\|\lesssim_{\mathbb{P}}\sqrt{K}s_{\beta}\left(\sqrt{\frac{s_{\gamma}\log{p}}{n}}+K^{-2\upsilon}\right)+K^{-\upsilon}+\sqrt{\frac{s_{\gamma}\log{p}}{n}}+\sqrt{\frac{K}{n}}\,.$

The above bound on the estimation error on $\omega_{b}$ translates to the
estimation error of $b^{\prime}$ as follows:

$\displaystyle\left|b^{\prime}(x)-\hat{b}^{\prime}(x)\right|$
$\displaystyle\leq\left\|\nabla
N_{K}(x)\right\|\left\|\hat{\omega}_{b}-\omega_{b}\right\|+K^{-(\upsilon-1)}$
$\displaystyle\lesssim_{\mathbb{P}}K^{3/2}\left[\sqrt{K}s_{\beta}\left(\sqrt{\frac{s_{\gamma}\log{p}}{n}}+K^{-2\upsilon}\right)+K^{-\upsilon}+\sqrt{\frac{s_{\gamma}\log{p}}{n}}+\sqrt{\frac{K}{n}}\right]+K^{-(\upsilon-1)}$
$\displaystyle\lesssim_{\mathbb{P}}K^{2}s_{\beta}\sqrt{\frac{s_{\gamma}\log{p}}{n}}+K^{-(\upsilon-\frac{3}{2})}+K^{\frac{3}{2}}\sqrt{\frac{s_{\gamma}\log{p}}{n}}+\frac{K^{2}}{\sqrt{n}}$
$\displaystyle\lesssim_{\mathbb{P}}K^{2}s_{\beta}\sqrt{\frac{s_{\gamma}\log{p}}{n}}+K^{-\left(\upsilon-\frac{3}{2}\right)}\,.$

Hence an optimal choice of $K$ would satisfy:

$\displaystyle K^{2}s_{\beta}\sqrt{\frac{s_{\gamma}\log{p}}{n}}\asymp
K^{-\left(\upsilon-\frac{3}{2}\right)}\implies\ $ $\displaystyle
s_{\beta}\sqrt{\frac{s_{\gamma}\log{p}}{n}}\asymp K^{-(\upsilon+\frac{1}{2})}$
$\displaystyle\implies\ $
$\displaystyle\frac{s^{2}_{\beta}s_{\gamma}\log{p}}{n}\asymp
K^{-(2\upsilon+1)}$ $\displaystyle\implies$ $\displaystyle
K\asymp\left(\frac{n}{s^{2}_{\beta}s_{\gamma}\log{p}}\right)^{\frac{1}{2\upsilon+1}}$

Using this we conclude:

$\sup_{|t|\leq\tau}\left|b^{\prime}(t)-\hat{b}^{\prime}(t)\right|\lesssim_{\mathbb{P}}\left(\frac{n}{s^{2}_{\beta}s_{\gamma}\log{p}}\right)^{\frac{-(\upsilon-\frac{3}{2})}{2\upsilon+1}}\,.$

This completes the proof. ∎

###  A.3 Estimation of $\alpha_{0}$ from $\mathcal{D}_{3}$

For notational simplicity, we here use $r_{n}(=K^{-\upsilon})$ to denote the
B-spline approximation of $m_{j}$ and $\check{m}_{j}$ (see the definitions in
Section 4) and $\dot{r}_{n}$ to denote the estimation error of
$\hat{b}^{\prime}$ obtained in Proposition A.1, i.e. we write:

  1. 1.

$\sup_{|t|\leq\tau}\left|b^{\prime}(t)-\hat{b}^{\prime}(t)\right|\lesssim_{\mathbb{P}}\dot{r}_{n}$.

  2. 2.

For $0\leq j\leq 1+p_{1}+p_{2}$,
$\sup_{|t|\leq\tau}\left|m_{j}(t)-\mathbf{N}_{k}(t)^{\top}\omega_{j}\right|\lesssim_{\mathbb{P}}r_{n}$.

  3. 3.

For $0\leq j\leq 1+p_{1}+p_{2}$,
$\sup_{|t|\leq\tau}\left|\check{m}_{j}(t)-\mathbf{N}_{k}(t)^{\top}\check{\omega}_{j}\right|\lesssim_{\mathbb{P}}r_{n}$

where $\omega_{j}$ and $\check{\omega}_{j}$ are the optimal projection vectors
of $m_{j}$ and $\check{m}_{j}$ respectively on the space spanned by $K$
B-spline basis with respect to $\ell_{\infty}$ norm. Henceforth, we will work
on the intersection of these events. First, Consider the LASSO regression of
$\mathbf{S}^{\perp}$ on $\breve{\mathbf{W}}_{-1}^{\perp}$. By basic inequality
we have:

$\displaystyle\frac{1}{2n}\left\|\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,S}\right\|^{2}+\lambda_{1}\|\hat{\theta}_{-1,S}\|_{1}\leq\frac{1}{2n}\left\|\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\theta^{*}_{S}\right\|^{2}+\lambda_{1}\|\theta^{*}_{S}\|_{1}$

Few algebraic manipulations (similar to that of standard LASSO analysis)
yields:

$\displaystyle\frac{1}{2n}\left\|\breve{\mathbf{W}}^{\perp}_{-1}\left(\hat{\theta}_{-1,S}-\theta^{*}_{S}\right)\right\|^{2}+\lambda_{1}\|\hat{\theta}_{-1,S}\|_{1}$
$\displaystyle\qquad\qquad\qquad\qquad\leq\frac{1}{n}\left\|\left(\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\theta^{*}_{S}\right)^{\top}\breve{\mathbf{W}}^{\perp}_{-1}\right\|_{\infty}\left\|\hat{\theta}_{-1,S}-\theta^{*}_{S}\right\|_{1}+\lambda_{1}\|\theta^{*}_{S}\|_{1}$
(A.7)

To find the optimal value of $\lambda_{1}$, we need to bound
$(1/n)\|(\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\theta^{*}_{S})^{\top}\breve{\mathbf{W}}^{\perp}_{-1}\|_{\infty}$.
Similarly, for the LASSO regression of $\mathbf{Y}$ on
$\breve{\mathbf{W}}_{-1}$, we have:

$\displaystyle\frac{1}{2n}\left\|\breve{\mathbf{W}}^{\perp}_{-1}\left(\hat{\theta}_{-1,Y}-\theta^{*}_{Y}\right)\right\|^{2}+\lambda_{0}\|\hat{\theta}_{-1,Y}\|_{1}$
$\displaystyle\qquad\qquad\qquad\qquad\leq\frac{1}{n}\left\|\left(\mathbf{Y}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\theta^{*}_{Y}\right)^{\top}\breve{\mathbf{W}}^{\perp}_{-1}\right\|_{\infty}\left\|\hat{\theta}_{-1,Y}-\theta^{*}_{Y}\right\|_{1}+\lambda_{0}\|\theta^{*}_{Y}\|_{1}$
(A.8)

and to obtain $\lambda_{0}$, we need to bound
$(1/n)\|(\mathbf{Y}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\theta^{*}_{Y})^{\top}\breve{\mathbf{W}}^{\perp}_{-1}\|_{\infty}$.
Towards that direction, we need the following lemma:

######  Lemma A.3.

Define the random variable $\widehat{\tilde{W}}$ (and consequently the matrix
$\widehat{\tilde{\mathbf{W}}}$) and $W-\mathbb{E}[W\mid\hat{\eta}]$. Under our
assumptions, we have for all $1\leq l\leq 1+p_{1}+p_{2}$ we have with
probability going to $1$:

$\max_{1\leq j\leq
1+p_{1}+p_{2}}\frac{1}{n}\left\|\breve{\mathbf{W}}^{\perp}_{*,j}-\widehat{\tilde{\mathbf{W}}}_{*,j}\right\|^{2}\lesssim_{\mathbb{P}}\
\dot{r}_{n}^{2}+\frac{s_{\gamma}\log{p}}{n}+\dot{r}_{n}\sqrt{\frac{\log{p}}{n}}$
$\max_{1\leq j\leq
1+p_{1}+p_{2}}\frac{1}{n}\left\|\breve{\mathbf{W}}^{\perp}_{*,j}-\tilde{\mathbf{W}}_{*,j}\right\|^{2}\lesssim_{\mathbb{P}}\
\dot{r}_{n}^{2}+\frac{s_{\gamma}\log{p}}{n}+\dot{r}_{n}\sqrt{\frac{\log{p}}{n}}\,,$

and

$\max_{1\leq j,j^{\prime}\leq
p}\frac{1}{n}\left|\left(\breve{\mathbf{W}}_{*,j}^{\perp}\right)^{\top}\left(\breve{\mathbf{W}}^{\perp}_{*,j^{\prime}}-\tilde{\mathbf{W}}_{*,j^{\prime}}\right)\right|\lesssim_{\mathbb{P}}\
\dot{r}_{n}+\sqrt{\frac{s_{\gamma}\log{p}}{n}}\,.$

The proof of this lemma can be found in the Supplementary document. An
immediate consequence of Lemma A.3 is the following bound which will be used
subsequently in this proof:

$\displaystyle\max_{1\leq j\neq j^{\prime}\leq
p}\frac{1}{n}\left|\left(\breve{\mathbf{W}}^{\perp}_{*,j}-\tilde{\mathbf{W}}_{*,j}\right)^{\top}\left(\breve{\mathbf{W}}^{\perp}_{*,j^{\prime}}-\tilde{\mathbf{W}}_{*,j^{\prime}}\right)\right|$
$\displaystyle\leq\sqrt{\max_{1\leq j\leq
p}\frac{1}{n}\left\|\breve{\mathbf{W}}^{\perp}_{*,j}-\tilde{\mathbf{W}}_{*,j}\right\|^{2}}\times\sqrt{\max_{1\leq
j^{\prime}\leq
p}\frac{1}{n}\left\|\breve{\mathbf{W}}^{\perp}_{*,j^{\prime}}-\tilde{\mathbf{W}}_{*,j^{\prime}}\right\|^{2}}$
$\displaystyle\lesssim_{\mathbb{P}}\
\dot{r}_{n}^{2}+\frac{s_{\gamma}\log{p}}{n}+\dot{r}_{n}\sqrt{\frac{\log{p}}{n}}\,.$
(A.9)

Based on the bound obtained in Lemma A.3, the optimal choices for
$\lambda_{0}$ and $\lambda_{1}$ are following:

######  Lemma A.4.

Under our assumptions, we can choose $\lambda_{0}$ and $\lambda_{1}$ as:

$\lambda_{1}\asymp(1+\|\theta^{*}_{S}\|_{1})\left[\dot{r}_{n}+\sqrt{\frac{s_{\gamma}\log{p}}{n}}\right]+\sqrt{\frac{s_{\gamma}\log{p}}{n}}\left(\log{\frac{n}{s_{\gamma}\log{p}}}\right)^{3/2}\,,$
$\lambda_{0}\asymp(1+\|\theta^{*}_{Y}\|_{1})\left[\dot{r}_{n}+\sqrt{\frac{s_{\gamma}\log{p}}{n}}\right]+\sqrt{\frac{s_{\gamma}\log{p}}{n}}\left(\log{\frac{n}{s_{\gamma}\log{p}}}\right)^{3/2}\,.$

where $\theta^{*}_{S},\theta^{*}_{Y}$ are same as defined in Section 4.

The proof of Lemma A.4 is also presented in the Supplementary document.
Another important ingredient in obtaining this LASSO-type bounds is the
restricted eigenvalue assumption on the covariate matrix
$\breve{\mathbf{W}}^{\perp}$ which is presented in the next Proposition (proof
is in Supplementary document):

######  Proposition A.5.

Under our assumptions in Section 4, the matrix
$\breve{\mathbf{W}}^{\perp}_{-1}$ satisfies RE condition with high
probability.

We next obtain the estimation error of $\hat{\theta}_{-1,Y}$ and
$\hat{\theta}_{-1,S}$ combining our findings from Lemma A.4 and Proposition
A.5. From (A.3) we have:

$\displaystyle\frac{1}{2n}\left\|\breve{\mathbf{W}}^{\perp}_{-1}\left(\hat{\theta}_{-1,S}-\theta^{*}_{S}\right)\right\|^{2}+\lambda_{1}\|\hat{\theta}_{-1,S}\|_{1}\leq\frac{\lambda_{1}}{2}\left\|\hat{\theta}_{-1,S}-\theta^{*}_{S}\right\|_{1}+\lambda_{1}\|\theta^{*}_{S}\|_{1}$

Therefore by the choice of $\lambda_{1}$ of Lemma A.4 and Proposition A.5 we
have (via standard LASSO bound calculation):

$\displaystyle\left\|\hat{\theta}_{-1,S}-\theta^{*}_{S}\right\|_{1}\lesssim_{\mathbb{P}}\lambda_{1}s_{1}\,,$
$\displaystyle\left\|\hat{\theta}_{-1,S}-\theta^{*}_{S}\right\|^{2}_{2}\lesssim_{\mathbb{P}}\lambda^{2}_{1}s_{1}\,.$

Similarly, for $\hat{\theta}_{-1,Y}$, we have from equation (A.3):

$\displaystyle\frac{1}{2n}\left\|\breve{\mathbf{W}}^{\perp}_{-1}\left(\hat{\theta}_{-1,Y}-\theta^{*}_{Y}\right)\right\|^{2}+\lambda_{0}\|\hat{\theta}_{-1,Y}\|_{1}\leq\frac{\lambda_{0}}{2}\left\|\hat{\theta}_{-1,Y}-\theta^{*}_{Y}\right\|_{1}+\lambda_{0}\|\theta^{*}_{Y}\|_{1}$

As before, the value of $\lambda_{0}$ from Lemma A.4 and Proposition A.5
yields:

$\displaystyle\left\|\hat{\theta}_{-1,Y}-\left(\theta^{*}_{S}\alpha_{0}+\theta^{*}_{0,-1}\right)\right\|_{1}\lesssim_{p}\lambda_{0}s_{0}\,,$
$\displaystyle\left\|\hat{\theta}_{-1,Y}-\left(\theta^{*}_{S}\alpha_{0}+\theta^{*}_{0,-1}\right)\right\|^{2}_{2}\lesssim_{p}\lambda_{0}^{2}s_{0}\,.$

From the above rates of the lasso estimates $\hat{\theta}_{-1,S}$ and
$\hat{\theta}_{-1,Y}$ we can further conclude:

$\displaystyle\left\|\theta_{0,-1}+\hat{\theta}_{-1,S}\alpha_{0}-\hat{\theta}_{-1,Y}\right\|_{1}$
$\displaystyle\leq\left\|\theta_{0,-1}+\theta^{*}_{S}\alpha_{0}-\hat{\theta}_{-1,Y}\right\|_{1}+|\alpha_{0}|\left\|\theta^{*}_{S}-\hat{\theta}_{-1,S}\right\|_{1}$
$\displaystyle\lesssim_{\mathbb{P}}|\alpha_{0}|\lambda_{1}s_{1}+\left\|\theta^{*}_{0,-1}+\hat{\theta}_{-1,S}\alpha_{0}-\hat{\theta}_{-1,Y}\right\|_{1}+\left\|\theta_{0,-1}-\theta^{*}_{0,-1}\right\|_{1}$
$\displaystyle\lesssim_{\mathbb{P}}\lambda_{0}s_{0}+|\alpha_{0}|\lambda_{1}s_{1}+s_{\gamma}\sqrt{\frac{\log{p}}{n}}$
(A.10)

and similarly:

$\displaystyle\left\|\theta^{0}_{-1}+\hat{\theta}_{-1,S}\alpha_{0}-\hat{\theta}_{-1,Y}\right\|^{2}_{2}$
$\displaystyle\lesssim_{p}\lambda^{2}_{0}s_{0}+\alpha^{2}_{0}\lambda^{2}_{1}s_{1}+\frac{s_{\gamma}\log{p}}{n}\,.$
(A.11)

Going back to equation (4.9), we next show in the following Proposition that
the common denominator of the three terms in the RHS stabilizes:

######  Proposition A.6.

Under our assumptions:

$\frac{1}{n_{3}}\left\|\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}\hat{\theta}_{-1,S}\right\|^{2}=\mathbb{E}\left[\left(\tilde{S}-\tilde{W}_{-1}^{\top}\theta^{*}_{S}\right)^{2}\mathds{1}_{|\eta|\leq\tau}\right]+o_{p}(1)$

The proof of Proposition A.6 can be found in the Appendix. We next show that
the numerators of the first and the second term of the RHS of equation (4.9)
are asymptotically negligible and the numerator of the third term contributes
to asymptotic normality.
Numerator of first term: We start with the first term, which is $\sqrt{n_{3}}$
times:

$\frac{1}{n_{3}}\left(\theta_{0,-1}+\hat{\theta}_{-1,S}\alpha_{0}-\hat{\theta}_{-1,Y}\right)^{\top}\breve{\mathbf{W}}_{-1}^{\perp^{\top}}\left(\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,S}\right)$

We first expand it as follows:

$\displaystyle\frac{1}{n}\left(\theta_{0,-1}+\hat{\theta}_{-1,S}\alpha_{0}-\hat{\theta}_{-1,Y}\right)^{\top}\breve{\mathbf{W}}_{-1}^{\perp^{\top}}\left(\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,S}\right)$
$\displaystyle\qquad\qquad=\underbrace{\frac{1}{n}\left(\theta_{0,-1}+\hat{\theta}_{-1,S}\alpha_{0}-\hat{\theta}_{-1,Y}\right)^{\top}\breve{\mathbf{W}}_{-1}^{\perp^{\top}}\left(\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\theta^{*}_{S}\right)}_{T_{1}}$
$\displaystyle\qquad\qquad\qquad+\underbrace{\frac{1}{n}\left(\theta_{0,-1}+\hat{\theta}_{-1,S}\alpha_{0}-\hat{\theta}_{-1,Y}\right)^{\top}\breve{\mathbf{W}}_{-1}^{\perp^{\top}}\breve{\mathbf{W}}_{-1}^{\perp}\left(\hat{\theta}_{-1,S}-\theta^{*}_{S}\right)}_{T_{2}}\,.$

To bound $T_{1}$ we use $\ell_{1}-\ell_{\infty}$ bound:

$\displaystyle\frac{1}{n}\left(\theta_{0,-1}+\hat{\theta}_{-1,S}\alpha_{0}-\hat{\theta}_{-1,Y}\right)^{\top}\breve{\mathbf{W}}_{-1}^{\perp^{\top}}\left(\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\theta^{*}_{S}\right)$
$\displaystyle\qquad\qquad\leq\left\|\theta_{0,-1}+\hat{\theta}_{-1,S}\alpha_{0}-\hat{\theta}_{-1,Y}\right\|_{1}\times\frac{1}{n}\left\|\breve{\mathbf{W}}_{-1}^{\top}\left(\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\theta^{*}_{S}\right)\right\|_{\infty}$

We have already established bound on the first part in equation (A.10) and the
second part is bounded by $\lambda_{1}$ (see the proof of Lemma A.4).
Therefore we have with probability going to 1:

$\displaystyle T_{1}$
$\displaystyle\lesssim_{\mathbb{P}}\lambda_{1}\left(\lambda_{0}s_{0}+|\alpha_{0}|\lambda_{1}s_{1}+s_{\gamma}\sqrt{\frac{\log{p}}{n}}\right)$
$\displaystyle\lesssim_{\mathbb{P}}\left(\lambda_{0}\vee\lambda_{1}\right)^{2}(s_{0}\vee
s_{1})+\lambda_{1}s_{\gamma}\sqrt{\frac{\log{p}}{n}}\,.$

For $T_{2}$ we can use CS inequality to conclude:

$\displaystyle T_{2}$
$\displaystyle\leq\underbrace{\sqrt{\frac{1}{n}\left\|\breve{\mathbf{W}}_{-1}^{\perp}\left(\hat{\theta}_{-1,S}-\theta^{*}_{S}\right)\right\|^{2}}}_{T_{21}}\times\underbrace{\sqrt{\frac{1}{n}\left\|\breve{\mathbf{W}}_{-1}^{\perp}\left(\theta_{0,-1}+\hat{\theta}_{-1,S}\alpha_{0}-\hat{\theta}_{-1,Y}\right)\right\|^{2}}}_{T_{22}}$

The first term $T_{21}$ is the prediction error of Lasso when we regress
$\mathbf{S}^{\perp}$ on $\breve{\mathbf{W}}^{\perp}_{1}$. Therefore from
standard Lasso prediction error bound we have:

$T_{21}\lesssim_{\mathbb{P}}\lambda_{1}\sqrt{s_{1}}\,,$

For $T_{22}$ we need a bit more detailed calculation. First of all note that,
$T_{22}$ can be further bounded as:

$\displaystyle T^{2}_{22}$
$\displaystyle=\frac{1}{n}\left\|\breve{\mathbf{W}}_{-1}^{\perp}\left(\theta_{0,-1}+\hat{\theta}_{-1,S}\alpha_{0}-\hat{\theta}_{-1,Y}\right)\right\|^{2}$
$\displaystyle\lesssim\frac{1}{n}\left\|\breve{\mathbf{W}}_{-1}^{\perp}\left(\theta_{0,-1}+\theta^{0}_{-1,1}\alpha_{0}-\hat{\theta}_{-1,Y}\right)\right\|^{2}+\frac{\alpha_{0}^{2}}{n}\left\|\breve{\mathbf{W}}_{-1}^{\perp}\left(\hat{\theta}_{-1,S}-\theta^{0}_{-1,1}\right)\right\|^{2}$
$\displaystyle\lesssim_{\mathbb{P}}\frac{1}{n}\left\|\breve{\mathbf{W}}_{-1}^{\perp}\left(\theta^{*}_{0,-1}+\theta^{0}_{-1,1}\alpha_{0}-\hat{\theta}_{-1,Y}\right)\right\|^{2}+\frac{1}{n}\left\|\breve{\mathbf{W}}_{-1}^{\perp}\left(\theta^{*}_{0,-1}-\theta^{0}_{-1}\right)\right\|^{2}+\lambda_{1}^{2}s_{1}$
$\displaystyle\lesssim_{\mathbb{P}}\lambda_{0}^{2}s_{0}+\lambda_{1}^{2}s_{1}+\frac{1}{n}\left\|P_{N_{k}}^{\perp}\tilde{\mathbf{Z}}\left(\hat{\gamma}_{n}-\gamma_{0}\right)\right\|^{2}$
$\displaystyle\lesssim_{\mathbb{P}}\lambda_{0}^{2}s_{0}+\lambda_{1}^{2}s_{1}+\frac{s_{\gamma}\log{p}}{n}\times\frac{1}{n}\sum_{i=1}^{n}\left(\tilde{Z}_{i}^{\top}a_{n}\right)^{2}$
$\displaystyle=\lambda_{0}^{2}s_{0}+\lambda_{1}^{2}s_{1}+\frac{s_{\gamma}\log{p}}{n}\times\frac{1}{n}\sum_{i=1}^{n}\left(\hat{b}^{\prime}(\hat{\eta}_{i})\right)^{2}\left(Z_{i}^{\top}a_{n}\right)^{2}$
$\displaystyle=\lambda_{0}^{2}s_{0}+\lambda_{1}^{2}s_{1}+\frac{s_{\gamma}\log{p}}{n}\times\left[\frac{1}{n}\sum_{i=1}^{n}\left(b^{\prime}(\hat{\eta}_{i})\right)^{2}\left(Z_{i}^{\top}a_{n}\right)^{2}+\frac{1}{n}\sum_{i=1}^{n}\left(\hat{b}^{\prime}(\hat{\eta}_{i})-b^{\prime}(\hat{\eta}_{i})\right)^{2}\left(Z_{i}^{\top}a_{n}\right)^{2}\right]$
$\displaystyle\lesssim_{\mathbb{P}}\lambda_{0}^{2}s_{0}+\lambda_{1}^{2}s_{1}+\frac{s_{\gamma}\log{p}}{n}\times\left[1+\dot{r}_{n}\right]\,.$

where in the last inequality, we use the following facts:

$\displaystyle b^{\prime}(\hat{\eta}_{i})$ $\displaystyle\lesssim 1\,,\hskip
14.45377pt[\text{As }b^{\prime}\text{ is bounded, Assumption
\ref{assm:smooth_b}}]$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\left(b^{\prime}(\hat{\eta}_{i})\right)^{2}\left(Z_{i}^{\top}a_{n}\right)^{2}$
$\displaystyle\lesssim_{\mathbb{P}}1\,,\hskip 14.45377pt[b^{\prime}\text{ is
bounded and Lemma \ref{lem:subg_second_moment}}]$
$\displaystyle\left\|\hat{b}^{\prime}-b^{\prime}\right\|_{\infty}$
$\displaystyle\lesssim_{\mathbb{P}}r^{\prime}_{n}\hskip
14.45377pt[\text{Proposition }\ref{prop:b_hat_est}]\,.$

Therefore we have:

$T_{22}\lesssim_{\mathbb{P}}\lambda_{0}\sqrt{s_{0}}+\lambda_{1}\sqrt{s_{1}}+\sqrt{\frac{s_{\gamma}\log{p}}{n}}\,.$

Taking products of the bounds on $T_{21}$ and $T_{22}$ we have:

$\displaystyle T_{2}\leq T_{21}\times T_{22}$
$\displaystyle\lesssim_{\mathbb{P}}\lambda_{1}\sqrt{s_{1}}\left(\lambda_{0}\sqrt{s_{0}}+\lambda_{1}\sqrt{s_{1}}+\sqrt{\frac{s_{\gamma}\log{p}}{n}}\right)$
$\displaystyle\lesssim_{\mathbb{P}}\left(\lambda_{0}\vee\lambda_{1}\right)^{2}(s_{0}\vee
s_{1})+\lambda_{1}\sqrt{\frac{s_{1}s_{\gamma}\log{p}}{n}}$

Combining bounds on $T_{1}$ and $T_{2}$ we obtain that for the second term in
equation (4.9), with probability going to 1:

$\displaystyle\frac{1}{n}\left(\theta^{0}_{-1}+\hat{\theta}_{-1,S}\alpha_{0}-\hat{\theta}_{-1,Y}\right)^{\top}\breve{\mathbf{W}}_{-1}^{\perp^{\top}}\left(\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,S}\right)$
$\displaystyle\qquad\qquad\lesssim_{\mathbb{P}}\left(\lambda_{0}\vee\lambda_{1}\right)^{2}(s_{0}\vee
s_{1})+\lambda_{1}s_{\gamma}\sqrt{\frac{\log{p}}{n}}+\left(\lambda_{0}\vee\lambda_{1}\right)^{2}(s_{0}\vee
s_{1})+\lambda_{1}\sqrt{\frac{s_{1}s_{\gamma}\log{p}}{n}}$
$\displaystyle\qquad\qquad\lesssim_{\mathbb{P}}\left(\lambda_{0}\vee\lambda_{1}\right)^{2}(s_{0}\vee
s_{1})+\lambda_{1}\sqrt{\frac{s_{\gamma}(s_{1}+s_{\gamma})\log{p}}{n}}$

which implies:

$\textbf{Numerator of Term 1
}\lesssim_{\mathbb{P}}\sqrt{n}\left[\left(\lambda_{0}\vee\lambda_{1}\right)^{2}(s_{0}\vee
s_{1})+\lambda_{1}\sqrt{\frac{s_{\gamma}(s_{1}+s_{\gamma})\log{p}}{n}}\right]\,.$
(A.12)

For the above term to be asymptotically negligible we need:

$\sqrt{n}\left[\left(\lambda_{0}\vee\lambda_{1}\right)^{2}(s_{0}\vee
s_{1})+\lambda_{1}\sqrt{\frac{s_{\gamma}(s_{1}+s_{\gamma})\log{p}}{n}}\right]=o(1)\,.$

Numerator of the second term: We next control $\sqrt{n}$ times the numerator
of second term on the RHS of equation (4.9) which involves the remainder
$\tilde{\mathbf{R}}$. Recall that
$\tilde{\mathbf{R}}=\mathbf{R}_{1}+\mathbf{R}_{2}+\mathbf{R}_{3}$ where:

$\displaystyle\mathbf{R}_{1}$
$\displaystyle=b(\hat{\eta})-\mathbf{N}_{k}(\hat{\eta})\omega_{b}$
$\displaystyle\mathbf{R}_{2}$
$\displaystyle=\left(b^{\prime}(\hat{\eta})-\hat{b}^{\prime}(\hat{\eta})\right)\odot
Z(\hat{\gamma}_{n}-\gamma_{0})$ $\displaystyle\mathbf{R}_{3}$
$\displaystyle=(\eta-\hat{\eta})^{2}\odot b^{\prime\prime}(\tilde{\eta})$

Therefore the numerator of the second term can be written as:

$\displaystyle\frac{1}{n}\tilde{\mathbf{R}}^{\perp^{\top}}\left(\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,S}\right)$
$\displaystyle=\underbrace{\frac{1}{n}\mathbf{R}_{1}^{\perp^{\top}}\left(\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,S}\right)}_{T_{1}}+\underbrace{\frac{1}{n}\mathbf{R}_{2}^{\perp^{\top}}\left(\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,S}\right)}_{T_{2}}$
$\displaystyle\qquad\qquad\qquad\qquad+\underbrace{\frac{1}{n}\mathbf{R}_{3}^{\perp^{\top}}\left(\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,S}\right)}_{T_{3}}\,.$

To bound $T_{1}$ note that $\mathbf{R}_{1}$ is measurable function of
$\hat{\eta}$. We further expand $T_{1}$ as:

$\displaystyle T_{1}$
$\displaystyle=\frac{1}{n}\mathbf{R}_{1}^{\perp^{\top}}\left(\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,S}\right)$
$\displaystyle=\underbrace{\frac{1}{n}\mathbf{R}_{1}^{\perp^{\top}}\left(\hat{\tilde{\mathbf{S}}}^{\perp}-\hat{\tilde{\mathbf{W}}}^{\perp}_{-1}\theta^{*}_{S}\right)}_{T_{11}}+\underbrace{\frac{1}{n}\mathbf{R}_{1}^{\perp^{\top}}\left(\mathbf{S}^{\perp}-\hat{\tilde{\mathbf{S}}}\right)}_{T_{12}}$
$\displaystyle\qquad\qquad\qquad+\underbrace{\frac{1}{n}\mathbf{R}_{1}^{\perp^{\top}}\left(\breve{\mathbf{W}}^{\perp}-\hat{\tilde{\mathbf{W}}}_{-1}\right)\theta^{*}_{S}}_{T_{13}}+\underbrace{\frac{1}{n}\mathbf{R}_{1}^{\perp^{\top}}\hat{\tilde{\mathbf{W}}}_{-1}\left(\theta^{*}_{S}-\hat{\theta}_{-1,S}\right)}_{T_{14}}$

Note that given $\hat{\eta}$, $T_{11}$ is a linear combination of centered
subgaussian random variables. Therefore by subgaussian concentration
inequality, we have with probabilities going to 1:

$T_{11}\lesssim_{\mathbb{P}}\frac{r_{n}}{\sqrt{n}}\,.$

Further we have as
$\mathbf{S}^{\perp}-\hat{\tilde{\mathbf{S}}}=P_{\mathbf{N}_{k}}^{\perp}\check{\mathbf{R}}_{1}-P_{\mathbf{N}_{k}}\check{\nu}_{1}$
we have with probability going to:

$T_{12}=\frac{1}{n}\mathbf{R}_{1}^{\perp^{\top}}\check{\mathbf{R}}_{1}\leq
r_{n}^{2}\,.$

For $T_{13}$ we apply a similar analysis along with $\ell_{1}-\ell_{\infty}$
bound:

$\displaystyle T_{13}$
$\displaystyle\leq\left\|\theta^{*}_{S}\right\|_{1}\frac{1}{n}\left\|\mathbf{R}_{1}^{\perp^{\top}}\left(\breve{\mathbf{W}}^{\perp}-\hat{\tilde{\mathbf{W}}}_{-1}\right)\right\|_{\infty}$
$\displaystyle=\left\|\theta^{*}_{S}\right\|_{1}\max_{2\leq j\leq
2p}\frac{1}{n}\left|\mathbf{R}_{1}^{\perp^{\top}}\left(\breve{\mathbf{W}}_{*,j}^{\perp}-\hat{\tilde{\mathbf{W}}}_{*,j}\right)\right|$
$\displaystyle\leq\left\|\theta^{*}_{S}\right\|_{1}\max_{2\leq j\leq
2p}\frac{1}{n}\left[\left|\mathbf{R}_{1}^{\perp^{\top}}\left(\breve{\mathbf{W}}_{*,j}^{\perp}-\mathbf{W}_{*,j}^{\perp}\right)\right|+\left|\mathbf{R}_{1}^{\perp^{\top}}\left(\mathbf{W}_{*,j}^{\perp}-\hat{\tilde{\mathbf{W}}}_{*,j}\right)\right|\right]$
$\displaystyle\leq\left\|\theta^{*}_{S}\right\|_{1}\max_{p+1\leq j\leq
2p}\frac{1}{n}\left|\mathbf{R}_{1}^{\perp^{\top}}\left(\breve{\mathbf{W}}_{*,j}^{\perp}-\mathbf{W}_{*,j}^{\perp}\right)\right|+\left\|\theta^{*}_{S}\right\|_{1}\max_{2\leq
j\leq
2p}\frac{1}{n}\left|\mathbf{R}_{1}^{\perp^{\top}}\check{\mathbf{R}}_{j}\right|$
$\displaystyle\lesssim\left\|\theta^{*}_{S}\right\|_{1}\max_{p+1\leq j\leq
2p}\sqrt{\frac{\mathbf{R}_{1}^{\perp^{\top}}\mathbf{R}_{1}}{n}}\sqrt{\frac{1}{n}\sum_{i}Z_{i,j}^{2}\left(\hat{b}(\hat{\eta}_{i})-b(\eta_{i})\right)^{2}}+\left\|\theta^{*}_{S}\right\|_{1}r_{n}^{2}$
$\displaystyle\lesssim_{\mathbb{P}}\left\|\theta^{*}_{S}\right\|_{1}\left[r_{n}\left(\dot{r}_{n}+\sqrt{\frac{s_{\gamma}\log{p}}{n}}\right)\right]+\left\|\theta^{*}_{S}\right\|_{1}r_{n}^{2}$
$\displaystyle\lesssim_{\mathbb{P}}\left\|\theta^{*}_{S}\right\|_{1}\left[r_{n}\left(\dot{r}_{n}+\sqrt{\frac{s_{\gamma}\log{p}}{n}}\right)\right]\,.$

And for the last term $T_{14}$ we have:

$\displaystyle T_{14}$
$\displaystyle=\frac{1}{n}\mathbf{R}_{1}^{\perp^{\top}}\hat{\tilde{\mathbf{W}}}_{-1}\left(\theta^{*}_{S}-\hat{\theta}_{-1,S}\right)$
$\displaystyle\leq\left\|\theta^{*}_{S}-\hat{\theta}_{-1,S}\right\|_{1}\max_{1\leq
j\leq
p}\frac{1}{n}\left|\mathbf{R}_{1}^{\perp^{\top}}\hat{\tilde{\mathbf{W}}}_{*,j}\right|$
$\displaystyle\lesssim\lambda_{1}s_{1}r_{n}\sqrt{\frac{\log{p}}{n}}\,.$

Combining the bounds on the different components of $T_{1}$ we have with
probability going to 1:

$\displaystyle
T_{1}\lesssim_{\mathbb{P}}\frac{r_{n}}{\sqrt{n}}+r_{n}^{2}+\left\|\theta^{*}_{S}\right\|_{1}\left[r_{n}\left(\dot{r}_{n}+\sqrt{\frac{s_{\gamma}\log{p}}{n}}\right)\right]+\lambda_{1}s_{1}r_{n}\sqrt{\frac{\log{p}}{n}}\,.$

Now we consider the second remainder term $T_{2}$ which involves
$\mathbf{R}_{2}$. An easy bound on this term will be following:

$\displaystyle
T_{2}=\frac{1}{n}\mathbf{R}_{2}^{\perp^{\top}}\left(\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,S}\right)$
$\displaystyle\leq\sqrt{\frac{1}{n}\mathbf{R}_{2}^{\top}\mathbf{R}_{2}}\sqrt{\frac{1}{n}\left\|\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,S}\right\|^{2}}\lesssim_{\mathbb{P}}\dot{r}_{n}\sqrt{\frac{s_{\gamma}\log{p}}{n}}\,.$

Finally, for the third term $T_{3}$ an easy CS bound is sufficient:

$T_{3}=\frac{1}{n}\mathbf{R}_{3}^{\perp^{\top}}\left(\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,S}\right)\leq\sqrt{\frac{1}{n}\mathbf{R}_{3}^{\top}\mathbf{R}_{3}}\sqrt{\frac{1}{n}\left\|\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,S}\right\|^{2}}\lesssim\frac{s_{\gamma}\log{p}}{n}\,.$

Combining the bounds on $T_{1},T_{2},T_{3}$ we have with probability going to
1:

$\displaystyle T_{1}+T_{2}+T_{3}$
$\displaystyle\lesssim_{\mathbb{P}}\frac{r_{n}}{\sqrt{n}}+r_{n}^{2}+\left\|\theta^{*}_{S}\right\|_{1}\left[r_{n}\left(\dot{r}_{n}+\sqrt{\frac{s_{\gamma}\log{p}}{n}}\right)\right]+\lambda_{1}s_{1}r_{n}\sqrt{\frac{\log{p}}{n}}$
$\displaystyle\qquad\qquad+\dot{r}_{n}\sqrt{\frac{s_{\gamma}\log{p}}{n}}+\frac{s_{\gamma}\log{p}}{n}\,.$

which implies:

Numerator 2nd
$\displaystyle\lesssim_{\mathbb{P}}\sqrt{n}\left[\frac{r_{n}}{\sqrt{n}}+r_{n}^{2}+\left\|\theta^{*}_{S}\right\|_{1}\left[r_{n}\left(\dot{r}_{n}+\sqrt{\frac{s_{\gamma}\log{p}}{n}}\right)\right]+\lambda_{1}s_{1}r_{n}\sqrt{\frac{\log{p}}{n}}\right.$
$\displaystyle\qquad\qquad\left.+\dot{r}_{n}\sqrt{\frac{s_{\gamma}\log{p}}{n}}+\frac{s_{\gamma}\log{p}}{n}\right]\,.$
(A.13)

Numerator of the third term: Last but not the least, we need to analyze the
numerator of the error term, i.e. third term on the RHS of equation (4.9)
First we show that:

$\frac{1}{\sqrt{n}}{\bm{\epsilon}}^{\perp^{\top}}\left(\mathbf{S}^{\perp}-\breve{\mathbf{W}}^{\perp}_{-1}\hat{\theta}_{-1,S}\right)=\frac{1}{\sqrt{n}}{\bm{\epsilon}}^{\top}\left(\tilde{\mathbf{S}}-\tilde{\mathbf{W}}_{-1}\theta^{*}_{S}\right)+o_{p}(1)\,.$

That the first term on the RHS is $O_{p}(1)$ is from the sub-exponential
concentration inequality which will later be shown to be asymptotically
normal. To show that the remainder is asymptotically negligible, first observe
that: