Semiparametric Difference-in-Differences with Potentially Many Control Variables
Neng-Chieh Chang





Semiparametric Difference-in-Differences with Potentially Many Control Variables
    Neng-Chieh Chang[Department of Economics, University of California Los Angeles, 315 Portola Plaza, Los Angeles, CA 90095, USA. email: nengchiehchang@g.ucla.edu]
    December 30, 2023
====================================================================================================================================================================


This paper discusses difference-in-differences (DID) estimation when
there exist many control variables, potentially more than the sample size. In this case, traditional
estimation methods, which require a limited number of variables, do not work. One may consider using statistical or machine learning
(ML) methods.  However, by the well-known theory of inference of ML methods proposed in <cit.>, directly applying ML methods to
the conventional semiparametric DID estimators will cause significant bias and make these
DID estimators fail to be √(N)-consistent. This article proposes three new DID estimators for three different data structures,
which are able to shrink the bias and achieve √(N)-consistency and asymptotic normality with mean zero when applying ML methods. This leads to straightforward inferential procedures.
In addition, I show that these new estimators have the small bias
property (SBP), meaning that their
bias will converge to zero faster than the pointwise bias of the nonparametric
estimator on which it is based. 




Keyword: difference-in-differences, causal inference, high-dimensional data, Neyman orthogonality, √(N)-consistency, undersmoothing





JEL Classification: C13, C14






§ INTRODUCTION


The difference-in-differences (DID) estimator has been widely used in empirical economics to evaluate causal effects when there exists a natural experiment with a treated group and an untreated group. By comparing the variation over time in an outcome variable between  the treated group and the untreated group, the DID estimator can be used to calculate the effect of treatment on the outcome variable. Applications of DID include but are not limited to studies of the effects
of immigration on labor markets <cit.>, the effects of minimum
wage law on wages <cit.>, the effect of tariffs liberalization
on corruption <cit.>, the effect of household income on children's
personalities <cit.>, and the effect of corporate tax on
wages <cit.>. 

The traditional linear DID estimator depends on a parallel trend
assumption that in the absence of treatment, the difference of outcomes between treated and untreated groups remains constant over time. In many situations, however, this assumption
may not hold because there are other individual characteristics that
may be associated with the variations of the outcomes. The treatment
may be taken as exogenous only after controlling these characteristics.
To address this problem, <cit.> proposed the semiparametric
DID estimators. Compared to the traditional linear DID estimators, the advantages of Abadie's estimators are threefold.
First, the characteristics are treated nonparametrically so that any estimation
error caused by functional specification is avoided. Second,
the effect of treatment is allowed to vary among individuals,
while the traditional linear DID estimator does not allow this heterogeneity.
Third, the estimation framework proposed in <cit.> allows
researchers to estimate how the effect of treatment varies with changes in the 
characteristics. 

This paper is an extension of <cit.>.  <cit.> considered
the case where the number of control variables has to be limited. A practical difficulty empirical researchers encounter is choosing what variables to
include when there is a rich data set. Although economic intuition can help us narrow down the choice set, it will not completely select all the important variables. This variable selection problem may lead to the chance of omitted variables in practice. 
In this paper, I consider the DID
estimation with many control variables, potentially more
than the sample size.  The classical estimation methods which
require a fixed number of variables do not work in this situation. One has to consider
using ML methods such as Lasso, Logit Lasso, random forests, boosted
trees, or various hybrids. However, by the well-known theory of inference of ML methods developed in <cit.>,  if one directly applies ML methods to the conventional semiparametric DID estimators proposed in <cit.>, the result will lead to significant bias and invalid inference. In particular,
the regularization bias embedded in ML methods will result in the
conventional semiparametric DID estimators failing to be √(N)-consistent. 

I contribute to the literature by proposing three new DID estimators
for three different data structures: repeated outcomes, repeated cross-sections,
and multilevel treatment. These new estimators can relieve the impact
of the regularization bias of ML methods and achieve √(N)-consistency. The key is to find the so-called Neyman-orthogonal
scores <cit.> of <cit.>'s estimands. The
Neyman-orthogonal score is a function that identifies the parameter
of interest, and its derivatives with respect to the nuisance parameters
are zero. This property helps us remove the first-order bias caused
by ML methods so that only the second-order bias remains, which
is much smaller and easier to control than the first-order bias as
in the conventional semiparametric DID estimators. Using the cross-fitting algorithm in 
<cit.>, I show that the new DID estimators can
be √(N)-consistent and asymptotically normal when using ML
methods. Figure 1 presents a Monte Carlo simulation that illustrates the negative effect of directly combining ML methods with Abadie's estimator and the benefit of using the newly proposed DID estimator.  



The second contribution is concerned with the conventional semiparametric DID estimators
with a limited number of control variables considered in <cit.>. In
this case, the conventional semiparametric DID estimators are able to achieve √(N)-consistency using kernel estimators, but they will require undersmoothing. Undersmoothing is a condition
that requires the pointwise bias of the kernel estimators to converge
to zero faster than the pointwise standard deviation. This condition will be violated if researchers use standard data-driven methods, such as cross-validation (CV), to choose the bandwidths of kernel estimators because those methods do not undersmooth. 

In this paper, I show that the new estimators do not require undersmoothing
to achieve √(N)-consistency. Specifically, I will show that the new estimators
have the small bias property (SBP), in terms of <cit.>, meaning that the bias of the new estimators will
converge to zero faster than the pointwise bias of the nonparametric
estimator on which it is based. The SBP, as shown in <cit.>, is a sufficient condition to remove the undersmoothing requirement. Figure 2 shows the Monte Carlo simulation results of Abadie's estimator and the new estimator with bandwidths chosen by CV. We can observe that Abadie's estimator is biased since CV does not undersmooth, and the newly proposed estimator can correct this bias.  


As an empirical example, I study the effect of tariff reduction on corruption behavior using the trade data between South Africa and Mozambique during 2006 and 2014. The treatment is the large tariff reduction on certain commodities occurring in 2008. This natural experiment was previously studied by <cit.> using the traditional linear DID estimator. I apply my proposed semiparametric DID estimator and  <cit.>'s semiparemetric DID estimator on the same data set (Table 9 of <cit.>). In comparison to <cit.> that a decrease in tariff rate will decrease corruption behavior, the two semiparametric estimators consistently suggest that the effect is actually substantially larger than previously reported by <cit.>. A potential explanation for this difference is that the true data generating process violates the linear specification assumed in the traditional linear DID estimator. In addition, when compared to <cit.>'s estimator, my proposed estimator shows that the effect is even larger.

The new estimators proposed in this paper heavily rely on the recent high-dimensional and
ML literature: <cit.>, <cit.>, <cit.>, <cit.>, and <cit.>; and the literature of the SBP in semiparametric estimation:   <cit.> and <cit.>.
 
Plan of the paper. Section 2 describes the conventional semiparametric
DID estimators and discusses their limitations when applying ML methods.
Section 3 presents the new DID estimators and discusses their theoretical
properties. Section 4 conducts Monte Carlo simulation to shed some light on the finite sample performance of the proposed estimators. Section
5 provides an application, and Section 6 concludes the paper. 



§ THE CONVENTIONAL SEMIPARAMETRIC DID ESTIMATORS


Let Y_i(t) be the outcome of interest for individual
i at time t and D_i(t)∈{ 0,1}
the treatment status. The population is observed in a pre-treatment period
t=0, and in a post-treatment period t=1. With potential outcome
notations <cit.>, we have Y_i(t)=Y_i^0(t)+(Y_i^1(t)-Y_i^0(t))D_i(t),
where Y_i^0(t) is the outcome that individual i
would attain at time t in the absence of the treatment, and Y_i^1(t)
represents the outcome that individual i would attain at time t
if exposed to the treatment. Since individuals are only exposed to
treatment at t=1, we have D_i(0)=0 for all i. To reduce
notation, I define D_i D_i(1). Also, let  X_i∈ℝ^d
be a vector of control variables with dimension d potentially larger than the
sample size N. 

The traditional linear DID estimator is the parameter α in the following linear model

    Y_i(t)=μ+X_i'π(t)+τ· D_i+δ· t+α· D_i(t)+ε_i(t),

where ε_i(t) is an exogenous shock that has mean
zero and (μ, π(t), τ, δ) are the corresponding parameters. Clearly, the linear specification assumed here
is a strong assumption since the true data generating process may be nonlinear. In addition, <cit.> noticed that including control variables in this linear form may not be appropriate if the treatment has different effects for different groups in the population. To deal with these problems, <cit.> proposed the semiparametric DID estimators which can
identify average treatment effect on the treated (ATT)

    θ_0 E[Y_i^1(1)-Y_i^0(1)| D_i=1].

According to the data, there are three particular cases. 



  §.§.§ Case 1: Random sample with repeated outcomes


Consider the case that researchers can observe both pre-treatment
and post-treatment outcomes for each individual of interest. That
is, researchers observe { Y_i(0),Y_i(1),D_i,X_i} _i=1^N. In this case, the ATT can be identified under the following assumptions <cit.>:


  Assumption 2.1. E[Y_i^0(1)-Y_i^0(0)| X_i,D_i=1]=E[Y_i^0(1)-Y_i^0(0)| X_i,D_i=0].

  Assumption 2.2. P(D_i=1)>0 and with probability
one P(D_i=1| X_i)<1.

Assumption (2.1) is the conditional parallel trend assumption. It states
that conditional on individual's characteristics, the average outcomes for treated
and untreated groups would have followed parallel paths in the absence of
treatment. With these two assumptions, the ATT is identified
<cit.> as 

    θ_0=E[Y_i(1)-Y_i(0)/P(D_i=1)D_i-P(D_i=1| X_i)/1-P(D_i=1| X_i)].2.1





  §.§.§ Case 2: Random sample with repeated cross sections


Often times, researchers may not be able to observe both pre-treatment
and post-treatment outcomes of the same individual. Instead, they
observe repeated cross-section data sets. Let T_i be a time indicator
that takes value one if the observation belongs to the post-treatment
sample. Researchers observe { Y_i,D_i,T_i,X_i} _i=1^N,
where Y_i=Y_i(0)+T_i(Y_i(1)-Y_i(0)).



  Assumption 2.3. Conditional on T=0, the data are i.i.d. from
the distribution of (Y(0),D,X); conditional
on T=1, the data are i.i.d. from the distribution of (Y(1),D,X).

Suppose Assumptions (2.1)-(2.3)
hold, the ATT is identified <cit.> as

    θ_0=E[T_i-λ_0/λ_0(1-λ_0)Y_i/P(D_i=1)D_i-P(D_i=1| X_i)/1-P(D_i=1| X_i)],2.2

where λ_0 P(T_i=1).


  §.§.§ Case 3: Multilevel treatments


In many cases, individuals can be exposed to different levels of 
treatment. Let W∈{0, w_1,...,w_J} be the level
of treatment, where W=0 denotes the untreated individuals. Researchers observe  { Y_i(0),Y_i(1),W_i,X_i} _i=1^N.

For w∈{ 0,w_1,...,w_J} and t∈{ 0,1},
let Y^w(t) be the potential outcome for treatment
level w at period t. Denote the ATT for each level of treatment w by 

    θ_0^w E[Y^w(1)-Y^0(1)| W=w].

Suppose that Assumptions (2.1) and (2.2) hold for each level of treatment:

    E[Y_i^0(1)-Y_i^0(0)| X_i,W_i=w]=E[Y_i^0(1)-Y_i^0(0)| X_i,W_i=0]

 for w∈{ w_1,...,w_J} and P(W_i=w)>0
and with probability one P(W_i=w| X_i)<1 for
w∈{ w_1,...,w_J}. Then we have <cit.>

    θ_0^w=E[Y(1)-Y(0)/P(W=w)I(W=w)· P(W=0| X)-I(W=0)· P(W=w| X)/P(W=0| X)],2.3

where I(·) is an indicator function. 


Let us focus on Case 1 in which researchers confront repeated
outcomes data. To use the identification result (2.1), the first step
is to estimate the two nuisance parameters: P(D_i=1) p_0
and P(D_i=1| X_i) g_0(X_i).
The estimator of p_0 is just a sample average p̂=N^-1∑_i=1^ND_i,
while the propensity score  g_0 is infinite-dimensional and needs to be estimated nonparametrically.
Denote by ĝ the estimator of  g_0,
then the plug-in estimator based on equation (2.1) is 

    θ̂=1/N∑_i=1^NY_i(1)-Y_i(0)/p̂D_i-ĝ(X_i)/1-ĝ(X_i).

When  ĝ is estimated using  classical nonparametric
methods such as kernel or series estimators, the estimator θ̂ can
be √(N)-consistent and asymptotically normal under certain
conditions provided in the semiparametric estimation literature <cit.>. 

When ĝ is an ML estimator, however, the estimator θ̂
will fail to be √(N)-consistent in general. By the general theory of inference of ML methods developed in <cit.>, the reason is twofold : 
(1) the score function based on (2.1), φ(W,θ_0,p_0,g_0)Y(1)-Y(0)/P(D=1)D-g_0(X)/1-g_0(X)-θ_0,
has a non-zero directional (Gateaux) derivative with respect to the propensity score g_0:

    ∂_g E[φ(W,θ_0,p_0,g_0)][g-g_0]≠0,
 where the directional (Gateaux) derivative is formally defined in Section 3; (2) ML estimators usually have a convergence
rate slower than N^-1/2 due to regularization bias. Similarly, the estimators obtained by directly plugging ML estimators into (2.2) and (2.3) will not be √(N)-consistent in general. The Monte Carlo simulation in Section 4 supports this theoretical insight and reveals significant bias on the estimators based on (2.1)-(2.3) when using ML estimators in the first-stage nonparametric estimation. 

The next section proposes three new score functions to relieve the regularization bias of the first-stage ML estimators. These three new score functions are derived under the same
identification assumptions as those in <cit.>, so that no extra assumption
is made. Heuristically, a distinctive feature of the new score functions
is that their derivatives with respect to their infinite-dimensional
nuisance parameters are zero. This property can help us  remove
the first-order bias of the first-stage estimation so that the bias
of the estimators based on these new score functions will be much smaller. In addition, I will use
the cross-fitting algorithm to improve the over-fitting phenomena
that frequently arise when using highly adaptive ML methods
<cit.>. 



§ THE NEW DID ESTIMATORS




 §.§ The Main Algorithm


Supposing Assumptions (2.1)-(2.3) hold, consider the following three
new score functions.



  §.§.§ Case 1: Random sample with repeated outcomes


The new score function for repeated outcomes is


    ψ_1(W,θ_0,p_0,η_10)     =Y(1)-Y(0)/P(D=1)D-P(D=1| X)/1-P(D=1| X)-θ_0
        -D-P(D=1| X)/P(D=1)(1-P(D=1| X))E[Y(1)-Y(0)| X,D=0]_c_1,3.1

with the unknown constant p_0 and the infinite-dimensional nuisance
parameter 

    η_10=(P(D=1| X),E[Y(1)-Y(0)| X,D=0])(g_0,ℓ_10).





  §.§.§ Case 2: Random sample with repeated cross sections


The new score function for repeated cross sections is

    ψ_2(W,θ_0,p_0,λ_0,η_20)=T-λ_0/λ_0(1-λ_0)Y/P(D=1)D-P(D=1| X)/1-P(D=1| X)-θ_0-c_2,3.2

where the adjustment term is 

    c_2=D-P(D=1| X)/λ_0(1-λ_0)· P(D=1)·(1-P(D=1| X))× E[(T-λ_0)Y| X,D=0].

The nuisance parameters are the unknown constants p_0 and λ_0,
and the infinite-dimensional parameter

    η_20=(P(D=1| X),E[(T-λ)Y| X,D=0])(g_0,ℓ_20).





  §.§.§ Case 3: Multilevel treatment


For each w∈{ w_1,...,w_J}, the new score function
for multilevel treatment is

    ψ_w(W,θ_w0,p_w0,η_w0)=    Y(1)-Y(0)/P(W=w)I(W=w)· P(W=0| X)-I(W=0)· P(W=w| X)/P(W=0| X)
        -θ_w0-c_w,3.3

where the adjustment term is 

    c_w=    (I(W=w)· P(W=0| X)-I(W=0)· P(W=w| X)/P(W=w)· P(W=0| X))×
        E[Y(1)-Y(0)| X,I(W=0)=1].

The nuisance parameters are the unknown constant p_w0 P(W=w)
and the infinite-dimensional parameter

    η_w0=(P(W=w| X),P(W=0| X),E[Y(1)-Y(0)| X,I(W=0)=1])(g_0w,g_0z,ℓ_30).


Notice that the above three new functions are equal to the original score
functions (2.1)-(2.3) plus the adjustment terms, (c_1,c_2,c_w),
which have zero expectations. Thus, the new score functions (3.1)-(3.3) still identify the ATT in each case. I will use these new scores to construct new DID estimators. 

To avoid repetition, I will focus on the estimation of ATT when data
belongs to repeated outcomes and repeated cross sections. The estimation
of multilevel treatment is provided in appendix. Now I combine the
score functions described above with the cross-fitting estimation
algorithm of <cit.>.


  Algorithm 1  



  * Take a K-fold random partition (I_k)_k=1^K
of observation indices [N]={ 1,...,N}. For simplicity, assume that each fold I_k has the same size  n=N/K. For each k∈[K]={ 1,...,K},
define the auxiliary sample I_k^c{ 1,...,N}∖ I_k. 


  * For each k, construct the intermediate ATT estimators 

    θ̃_k=1/n∑_i∈ I_kD_i-ĝ_k(X_i)/p̂_k(1-ĝ_k(X_i))×(Y_i(1)-Y_i(0)-ℓ̂_1k(X_i)) (repeated outcomes)


    θ̃_k=1/n∑_i∈ I_kD_i-ĝ_k(X_i)/p̂_kλ̂_k(1-λ̂_k)(1-ĝ_k(X_i))×((T_i-λ̂_k)Y_i-ℓ̂_2k(X_i)) (repeated cross sections)

where p̂_k=1/n∑_i∈ I_k^cD_i
, λ̂_k=1/n∑_i∈ I_k^cT_i, and (ĝ_k,ℓ̂_1k,ℓ̂_2k) are the estimators of (g_0,ℓ_10,ℓ_20) constructed using the auxiliary sample I_k^c. 

  * Construct the final ATT estimator θ̃=1/K∑_k=1^Kθ̃_k.

The estimators (ĝ_k,ℓ̂_1k,ℓ̂_2k)
can be constructed using any ML methods or classical estimators such
as kernel or series estimators. For completeness, I present the Logit
Lasso and Lasso estimators here. 

Consider a class of approximating functions of X_i,

    q_i(q_i1(X_i),...,q_ip(X_i))'.

For example, q_i can be polynomials or
B-splines. Let Λ(u)1/(1+exp(-u)) be
the cumulative distribution function of the standard Logistic distribution, construct
the estimator of the propensity score g_0 by 

    ĝ_k(x_i)Λ(q_i'β̂_k),3.4

where 

    β̂_kmin_β∈ℝ^p1/M∑_i∈ I_k^c{ -D_i(q_i'β)+log(1+exp(q_i'β))} +λ_k∥β∥_1

is the Logit Lasso estimator and M=N-n is the sample size of the auxiliary sample  I_k^c.
Next, define I_kz^c I_k^c∩{ i:D_i=0},
M_k the sample size of I_kz^c. Construct
the estimators of ℓ_10 and ℓ_20 by 

    ℓ̂_1k(x_i) q_i'β̂_1k,


    ℓ̂_2k(x_i) q_i'β̂_2k,

where 

    β̂_1k∈min_β∈ℝ^p[1/M_k∑_i∈ I_kz^c(Y_i(1)-Y_i(0)-q_i'β)^2]+λ_1k/M_k∥Υ̂_1kβ∥_1

and 

    β̂_2k∈min_β∈ℝ^p[1/M_k∑_i∈ I_kz^c((T_i-λ̂_k)Y_i-q_i'β)^2]+λ_2k/M_k∥Υ̂_2kβ∥_1

 are the modified Lasso estimators proposed in <cit.>.
The choices of the penalty levels and loadings (λ_1k,λ_2k,Υ̂_1k,Υ̂_2k)
suggested by <cit.> are provided in appendix.  



 §.§ Theoretical Properties


In this section, I discuss the theoretical properties of the new DID
estimator θ̃. In particular, I will show that the estimator
θ̃ can achieve √(N)-consistency and asymptotic
normality as long as the first-stage estimators converge at rates
faster than N^-1/4. This rate of convergence can be achieved
by many ML methods, including Lasso and Logit Lasso. Further, I will
show that when using kernel estimators in the first-stage estimation,
the estimator θ̃ has the SBP while the conventional semiparametric
DID estimators do not. 



  §.§.§ The Neyman Orthogonality 


The differences between the new DID estimators and the conventional semiparametric DID estimators in <cit.> are the score functions on which they are based. The key property of the new score functions (3.1)-(3.3) is that their directional (or the Gateaux) derivatives with
respect to their infinite-dimensional nuisance parameters are zero, while the scores based on (2.1)-(2.3) do not have this property.
This property is the so-called Neyman orthogonality in <cit.>.
The Neyman orthogonality enables us to remove the first-order bias of the first-stage estimation so that the estimators based on these
Neyman-orthogonal scores can achieve √(N)-consistency under
less restrictive conditions. 

The definition of the Neyman-orthogonal score provided here is slightly
different from <cit.> that instead of being orthogonal
against all nuisance parameters, the Neyman-orthogonal score defined
here is orthogonal against only those infinite-dimensional nuisance
parameters. Formally, let θ_0∈Θ be the low-dimensional
parameter of interest, ρ_0 be the true value of the finite-dimensional
nuisance parameter ρ, and η_0 the true value of the
infinite-dimensional nuisance parameter η∈𝒯 . Suppose
that W is a random element taking values in a measurable space
(𝒲,𝒜_𝒲) with probability
measure P. Define the directional (or the Gateaux) derivative against
the infinite-dimensional nuisance parameter D_r:𝒯̃→ℝ,
where 𝒯̃={η-η_0:η∈𝒯},

    D_r[η-η_0]∂_r{ E_P[ψ(W,θ_0,ρ_0,η_0+r(η-η_0))]} ,η∈𝒯,

for all r∈[0,1). For convenience,  denote

    ∂_ηE_Pψ(W,θ_0,ρ_0,η_0)[η-η_0] D_0[η-η_0],η∈𝒯.

In addition, let 𝒯_N⊂𝒯 be a nuisance
realization set such that the estimator of η_0 take values
in this set with high probability.


  Definition (The Neyman Orthogonality)

The score ψ obeys the Neyman orthogonality condition at
(θ_0,ρ_0,η_0) with respect to the nuisance
parameter realization set 𝒯_N⊂𝒯 if
the directional derivative map D_r[η-η_0]
exists for all r∈[0,1) and η∈𝒯_N and vanishes
at r=0:

    ∂_ηE_Pψ(W,θ_0,ρ_0,η_0)[η-η_0]=0,for all η∈𝒯_N.




  Lemma 1 The new score functions (3.1)-(3.3) obey the Neyman
orthogonality.


This property embedded in (3.1)-(3.3) will play the key role to make less restrictive assumptions in the following proofs of asymptotic distribution and the SBP. 



  §.§.§ Asymptotic Distribution


In the following, I will discuss the theoretical properties of the
new estimator θ̃ when data belongs to repeated outcomes
and repeated cross sections. The results of multilevel treatment can be proven using the same arguments. Let κ and C be
strictly positive constants, K≥2 be a fixed integer, and ε_N
be a sequence of positive constants approaching zero. Denote by ∥·∥_P,q the L^q norm of some probability measure P:  ∥ f ∥_P,q(∫| f(w) |^qdP(w))^1/q and ∥ f ∥_P,∞sup_w| f(w) |.


  Assumption 3.1 (Regularity Conditions for Repeated Outcomes)

Let P be the probability law for (Y(0),Y(1),D,X).
Let D=g_0(X)+U and Y(1)-Y(0)=ℓ_10(X)+V_1
with E_P[U| X]=0 and E_P[V_1| X, D=0]=0.
Define G_1p0 E_P[∂_pψ_1(W,θ_0,p_0,η_10)]
and Σ_10 E_P[(ψ_1(W,θ_0,p_0,η_10)+G_1p0(D-p_0))^2].
Suppose the following conditions hold: (a) Pr(κ≤ g_0(X)≤1-κ)=1;
(b) ∥ UV_1∥_P,4≤ C; (c) E[U^2| X]≤ C;
(d) E[V_1^2| X]≤ C; (e) Σ_10>0;
and (f) given the auxiliary sample I_k^c, the estimator η̂_1k=(ĝ_k,ℓ̂_1k)
obeys the following conditions. With probability 1-o(1),
∥η̂_1k-η_10∥_P,2≤ε_N,
∥ĝ_k-1/2∥_P,∞≤1/2-κ,
and ∥ĝ_k-g_0∥_P,2^2+∥ĝ_k-g_0∥_P,2×∥ℓ̂_1k-ℓ_10∥_P,2≤(ε_N)^2. 


  Assumption 3.2 (Regularity Conditions for Repeated Cross Sections)

Let P be the probability law for (Y,T,D,X). Let
D=g_0(X)+U and (T-λ_0)Y=ℓ_20(X)+V_2
with E_p[U| X]=0 and E_p[V_2| X, D=0]=0.
Define G_2p0 E_P[∂_pψ_2(W,θ_0,p_0,λ_0,η_20)],
G_2λ0 E_P[∂_λψ_2(W,θ_0,p_0,λ_0,η_20)],
and Σ_20 E_P[(ψ_1(W,θ_0,p_0,η_10)+G_2p0(D-p_0)+G_2λ0(T-λ_0))^2].
Suppose the following conditions hold: (a) Pr(κ≤ g_0(X)≤1-κ)=1;
(b) ∥ UV_2∥_P,4≤ C; (c) E[U^2| X]≤ C;
(d) E[V_2^2| X]≤ C; (e) E_P[Y^2| X]≤ C;
(f) | E_P[YU]|≤ C; (g) Σ_20>0;
and (h) given the auxiliary sample I_k^c, the estimators η̂_2k=(ĝ_k,ℓ̂_2k)
obeys the following conditions. With probability 1-o(1),
∥η̂_2k-η_20∥_P,2≤ε_N,
∥ĝ_k-1/2∥_P,∞≤1/2-κ,
and ∥ĝ_k-g_0∥_P,2^2+∥ĝ_k-g_0∥_P,2×∥ℓ̂_2k-ℓ_20∥_P,2≤(ε_N)^2. 


  Theorem 1  

For repeated outcomes, suppose Assumptions (2.1), (2.2) and (3.1) hold. For repeated cross sections, suppose Assumptions  (2.1)-(2.3) and (3.2) hold. If ε_N=o(N^-1/4),
then the new ATT estimator θ̃ satisfies

    √(N)(θ̃-θ_0)→ N(0,Σ)

with Σ=Σ_10 for repeated outcomes and Σ=Σ_20
for repeated cross sections. 


  Theorem 2 (Variance Estimator) 

Construct the estimators of the asymptotic variances as

    Σ̂_1=1/K∑_k=1^K𝔼_n,k[(ψ_1(W,θ̃,p̂_k,η̂_1k)+Ĝ_1p(D-p̂_k))^2]repeated outcomes


    Σ̂_2=1/K∑_k=1^K𝔼_n,k[(ψ_2(W,θ̃,p̂_k,λ̂_k,η̂_2k)+Ĝ_2p(D-p̂_k)+Ĝ_2λ(T-λ̂_k))^2]repeated cross sections

where 𝔼_n,k[f(W)]=n^-1∑_i∈ I_kf(W_i),
Ĝ_1p=Ĝ_2p=-θ̃/p̂_k, and Ĝ_2λ
is a consistent estimator of G_2λ0. If the assumptions
of Theorem 1 hold, then Σ̂_1=Σ_10+o_P(1)
and Σ̂_2=Σ_20+o_P(1).


The interpretation of Theorem 1 and 2 is that the new DID  estimator θ̃
can achieve √(N)-consistency and asymptotic normality provided that
the first-stage estimators of the infinite dimensional nuisance parameters
converge at a rate faster than N^-1/4. This rate of convergence
can be achieved by many ML methods. In particular, <cit.>
and <cit.> provided detail conditions for Logit Lasso
and the modified Lasso estimators to satisfy this rate of convergence. It is also
worth noting that even when the first-stage estimators do not converge
as fast as N^-1/4, the new estimator θ̃ still
has smaller bias than the original estimator because the Neyman orthogonality
removes the first-order bias of the first-stage estimators. 



  §.§.§ The Small Bias Property


Consider the conventional semiparametric DID estimators with a limited number
of control variables studied in <cit.>. Let g_h
be the kernel estimator of g_0 with bandwidth h→0 in (2.1) and (2.2). Under the standard assumptions of kernel estimation (Assumption (3.3) below), one can show that the pointwise bias of ĝ_h is of order O(h^m), where m can be interpreted as the minimum number of derivatives of g_0; and the pointwise standard deviation is sd(ĝ_h(x))=O((Nh^d+2s)^-1/2). By Theorem 8.11 of <cit.>, one can show that the √(N)-consistency of the plug-in estimators based on (2.1) and (2.2) requires √(N)h^m→ 0. That is, the pointwise bias of the kernel estimator has to converge to zero faster than N^-1/2.   Since the pointwise standard deviation converges to zero slower than N^-1/2, undersmoothing is required. In this case, standard data-driven bandwidth selection methods which do not undersmooth, such as cross-validation, are invalid.


To avoid undersmoothing, by the analysis of SBP in <cit.>,  the estimator of the parameter of interest needs to have smaller bias than the pointwise bias of the first-stage nonparametric estimators. That is, the SBP requires that the bias of the estimator of θ_0 converges to zero faster than h^m.

In the following, I will show that the new DID estimator θ̃
has the SBP. Let (ĝ_kh,ℓ̂_1kh,ℓ̂_2kh)
be the kernel estimators of (g_0,ℓ_10,ℓ_20)
using auxiliary sample I_k^c. I assume here that they have
the same bandwidth h and kernel K(u) for convenience.



  Assumption 3.3 <cit.>



  * K(u) is differentiable of order s, the derivatives
of order s are bounded, K(u) is zero outside a bounded
set, ∫ K(u)du=1, there is a positive m such
that for all j<m, ∫ K(u)[⊗_ℓ=1^ju]du=0.

  * Define γ_0(x)=f_0(x)E(z| x),
where z∈(1,D,Y(1)-Y(0)| D=0,(T-λ_0)Y| D=0)
and f_0(x) is the true density of x. Assume that
γ_0(x) is continuously differentiable to order
s with bounded derivatives on an open set containing  𝒳, where 𝒳 is the support of x.

  * There is α≥4 such that E[| z|^α]<∞
and E[| z|^α| x]f_0(x)
is bounded.



  Theorem 3 

For repeated outcomes, suppose Assumptions (2.1), (2.2),  (3.1), and (3.3) hold. For repeated cross sections, suppose Assumptions (2.1)-(2.3), (3.2), and (3.3) hold. Suppose
that  inf_x∈𝒳f_0(x)≠ 0, h=h(N) with  log N/(√(N)h^d+2s)→ 0. If √(N)h^2m→0,
then

    √(N)(θ̃-θ_0)→ N(0,Σ)

with Σ=Σ_10 for repeated outcomes and Σ=Σ_20
for repeated cross sections. 


The interpretation of Theorem 3 is that the new estimator θ̃
only requires √(N)h^2m→0 to achieve √(N)-consistency,
while the conventional semiparametric DID estimators require √(N)h^m→0
under the same assumptions. With the Neyman orthogonality,
the bias of θ̃ is only of the second-order of the pointwise
bias of the first-stage kernel estimators. The bias of θ̃
is h^2m instead of h^m. Hence, θ̃ satisfies  the SBP. In particular, the bandwidth h such that  log N/(√(N)h^d+2s)→ 0 and √(N)h^2m→0 exists only if 2m>d+2s. Under this condition, the optimal bandwidth selected by  minimizing  mean-square errors (CV),  h=N^-1/(d+2s+2m),  satisfies the conditions for √(N)-consistency. 


  Theorem 4  

Construct the estimators of the asymptotic variances as

    Σ̂_1=1/K∑_k=1^K𝔼_n,k[(ψ_1(W,θ̃,p̂_k,η̂_1kh)+Ĝ_1p(D-p̂_k))^2] repeated outcomes


    Σ̂_2=1/K∑_k=1^K𝔼_n,k[(ψ_2(W,θ̃,p̂_k,λ̂_k,η̂_2kh)+Ĝ_2p(D-p̂_k)+Ĝ_2λ(T-λ̂_k))^2] repeated cross sections

where Ĝ_1p=Ĝ_2p=-θ̃/p̂_k and
Ĝ_2λ is a consistent estimator of G_2λ0.
If the assumptions of Theorem 3 hold, then Σ̂_1=Σ_10+o_P(1)
and Σ̂_2=Σ_20+o_P(1).





§ SIMULATION


In this section, I present Monte Carlo simulation results of the conventional semiparametric DID estimators and the new DID estimator
θ̃ in three different data structures: repeated outcomes, repeated cross
sections, and multilevel treatment. I use both ML methods and kernel estimators in the first-stage estimation. For ML estimation, I generate high-dimensional (HD) data and estimate the propensity score by Logit Lasso (Multi-Logit Lasso for multilevel treatment). To choose the penalty parameter for Logit Lasso (Multi-Logit Lasso), I use K-fold CV (as recommended by <cit.>) with K=10. Alternatively, one could use a method developed in <cit.>. The other infinite-dimensional nuisance parameters are estimated by random forests with 500 regression trees. For kernel estimation, all the infinite-dimensional nuisance parameters are estimated using the standard Gaussian kernel. 

Figure 3-20 in appendix show the simulation results. I find that the conventional semiparametric DID estimators are biased when using ML methods, while the new DID estimator θ̃ can correct the bias. For kernel estimation, the conventional DID estimator with bandwidth selected by CV is biased, while the new DID estimator θ̃ is centered at the true value. The data generating processes are presented in the following. 



 §.§ Repeated Outcomes




  §.§.§ ML Estimation


Let N∈{200,500} be the sample size and p∈{100,300} the dimension of control variables,
X_i∼ N(0,I_p× p). Also, let γ_0=(1,1/2,1/3,1/4,1/5,0,...,0)∈ℝ^p
and D_i is generated by the propensity score 

    P(D=1| X)=1/1+exp(-X'γ_0)(Logistic).

At t=0, the potential outcome is generated 

    Y_i^0(0)=X_i'β_0+ε_1,

and at t=1,

    Y_i^0(1)=Y_i^0(0)+1+ε_2,


    Y_i^1(1)=θ_0+Y_i^0(1)+ε_3,

where β_0=γ_0+0.5 and θ_0=3, and all error
terms follow N(0,0.1). Researchers observe { Y_i(0),Y_i(1),D_i,X_i} for i=1,...,N,
where Y_i(0)=Y_i^0(0) and Y_i(1)=Y_i^0(1)(1-D_i)+Y_i^1(1)D_i. Figure 3-6 present the results.  



  §.§.§ Kernel Estimation


Let N∈{200,500} be the sample size, D_i∼ Bernoulli(0.5), and X_i| D_i∼ N(D_i,1). At t=0, the potential outcome is generated 

    Y_i^0(0)=ε_1,

and at t=1,

    Y_i^0(1)=Y_i^0(0)+X_i+ε_2,


    Y_i^1(1)=θ_0+Y_i^0(1)+ε_3,

where θ_0=3 and all error
terms follow N(0,0.1). Researchers observe { Y_i(0),Y_i(1),D_i,X_i} for i=1,...,N,
where Y_i(0)=Y_i^0(0) and Y_i(1)=Y_i^0(1)(1-D_i)+Y_i^1(1)D_i. Figure 7-8 present the results.



 §.§ Repeated Cross Sections




  §.§.§ ML Estimation


Let N∈{200,500} be the sample size and p∈{100,300} the dimension of control variables,
X_i∼ N(0.3,I_p× p). Also, let γ_0=(1,1/2,1/3,1/4,1/5,0,...,0)∈ℝ^p
and D is generated by the propensity score 

    P(D=1| X)=1/1+exp(-X'γ_0)(Logistic).

At t=0, the potential outcome is generated 

    Y_i^0(0)=1+ε_1,

and at t=1,

    Y_i^0(1)=Y_i^0(0)+1+ε_2,


    Y_i^1(1)=θ_0+Y_i^0(1)+ε_3,

where β_0=γ_0+0.5 and θ_0=3, and all error
terms follow N(0,0.1). Define Y_i(0)=Y_i^0(0)
and Y_i(1)=Y_i^0(1)(1-D_i)+Y_i^1(1)D_i.
Let T_i follow a Bernoulli distribution with parameter 0.5.
Researchers observe { Y_i,T_i,D_i,X_i} for i=1,...,N,
where Y_i=Y_i(0)+T_i(Y_i(1)-Y_i(0)). Figure 9-12 present the results. 



  §.§.§ Kernel Estimation


Let N∈{200,500} be the sample size, D_i∼ Bernoulli(0.5), and X_i| D_i∼ N(D_i,1). At t=0, the potential outcome is generated 

    Y_i^0(0)=ε_1,

and at t=1,

    Y_i^0(1)=Y_i^0(0)+X_i+ε_2,


    Y_i^1(1)=θ_0+Y_i^0(1)+ε_3,

where θ_0=3 and all error
terms follow N(0,0.1). Let Y_i(0)=Y_i^0(0) and Y_i(1)=Y_i^0(1)(1-D_i)+Y_i^1(1)D_i.
Let T_i∼ Bernoulli(0.5).
Researchers observe { Y_i,T_i,D_i,X_i} for i=1,...,N,
where Y_i=Y_i(0)+T_i(Y_i(1)-Y_i(0)). Figure 13-14 present the results.



 §.§ Multilevel Treatment




  §.§.§ ML Estimation


Suppose there are two levels of treatment so that W∈{ 0,1,2}. Let N∈{200,500} be the sample size and p∈{100,300} the dimension of control variables, X_i∼ N(0,I_p× p).
Also, let γ_0=(1,1/2,1/3,1/4,1/5,0,...,0)∈ℝ^p and

    (P(W=0),P(W=1),P(W=2))=(0.3,0.3,0.4)

At t=0, the potential outcome is generated 

    Y_i^0(0)=X'β_0+ε_1,

and at t=1,

    Y_i^0(1)=Y_i^0(0)+1+ε_2,


    Y_i^1(1)=θ_10+Y_i^0(1)+ε_3,


    Y_i^2(1)=θ_20+Y_i^0(1)+ε_4,

where β_0=γ_0+0.5 and θ_10=3 and θ_20=6,
and all error terms follow N(0,0.1). Researchers observe
{ Y_i(0),Y_i(1),W_i,X_i} for i=1,...,N,
where Y_i(0)=Y_i^0(0) and Y_i(1)=Y_i^0(1)I(W_i=0)+Y_i^1(1)I(W_i=1)+Y_i^2(1)I(W_i=2).
I focus on the estimation of the second level ATT θ_20. Figure 15-18 present the results


  §.§.§ Kernel Estimation


Suppose there are two levels of treatment so that W∈{ 0,1,2}. Let N be the sample
size, X_i| W_i∼ N(W_i,1), and


    P(W_i=0)=P(W_i=1)=P(W_i=2)=1/3.

At t=0, the potential outcome is generated 

    Y_i^0(0)=ε_1,

and at t=1,

    Y_i^0(1)=Y_i^0(0)+X_i+ε_2,


    Y_i^1(1)=θ_10+Y_i^0(1)+ε_3,


    Y_i^2(1)=θ_20+Y_i^0(1)+ε_4,

where  θ_10=3, θ_20=6,
and all error terms follow N(0,0.1). 
Let Y_i(0)=Y_i^0(0) and Y_i(1)=Y_i^0(1)I(W_i=0)+Y_i^1(1)I(W_i=1)+Y_i^2(1)I(W_i=2). Researchers observe
{ Y_i(0),Y_i(1),W_i,X_i} for i=1,...,N.
I focus on the estimation of the second level ATT θ_20. Figure 19-20 present the results. 



§ EMPIRICAL EXAMPLE


In this example, I analyze the effect of tariffs reduction on corruption behaviors using the bribe payment data collected by <cit.> between South Africa and Mozambique. There have been theoretical and empirical debates on whether higher tariff rates increase incentives for corruption to occur <cit.> or lower tariffs encourage agents to pay higher bribes through an income effect <cit.>. The former argues that an increase in the tariff rate makes it more profitable to evade taxes on the margin. The latter argues that an increased tariff rate makes the tax payer less wealthy and this, under the decreasing risk aversion of being penalized, tend to reduce evasion <cit.>. 

<cit.> collected primary data on the bribed payments between the ports in Mozambique and South Africa from 2007 to 2013. The treatment is the large reduction in the average nominal tariff rate (of 5 percent) occurring in 2008.    Since not all products were on the tariff reduction list, a credible control group of products is available. This allows for a DID estimation.  

This natural experiment between South Africa and Mozambique was previously studied by <cit.> by pooling the cross section data between 2007 and 2013, with sample size N=1084, and estimating the effect of treatment using the traditional linear DID. Here I focus on the specification of one of the main results (Table 9 of  <cit.>):

    y_it=    γ_1TariffChangeCategory_i× POST
        +μ POST+γ_2TariffChangeCategory_i
        +β_2BaselineTariff_i+Γ_i
        +p_i+w_t+δ_i+ϵ_it
 where y_it is the natural log of the amount of bribe paid for shipment i in period t, conditional on paying a bribe. TariffChangeCategory∈{ 0,1}  denotes the treatment status of commodities, POST∈{ 0,1}  is an indicator for the years following 2008, and BaselineTariff is the tariff rate before the tariff reduction. The specification also includes a vector of characteristics Γ_i, and time and individual fixed effects p_i, w_t, and δ_i. The parameter γ_1 is the parameter of interest in the traditional linear DID estimation. <cit.> found that the amount of bribes paid dropped after the tariff reduction (γ̂_1=-2.928^**).

I use the same data set but instead of using the traditional linear DID estimation, I estimate the ATT by <cit.>'s DID estimator and my proposed DID estimator θ̃. Since the data is repeated cross sections, I construct the estimators based on (2.2) and (3.2), respectively. The estimators with the first-stage kernel estimation contain one individual characteristic (the natural log of shipment value per ton), which is a significant characteristic in Table 9 of <cit.>. The estimators with the first-stage Lasso estimation contain a list of the significant characteristics in Table 9 of <cit.>, which includes product, shipment, firm-level characteristics,  and their interaction terms. Table 1 below shows the results. I find that all these estimators consistently suggest that a decrease in tariff rate will lead to less bribes payment, but the effect of treatment may be actually substantially larger than previously reported by <cit.>.






§ CONCLUSION

In this article, I have introduced three new DID estimators based on the newly-derived Neyman-orthogonal scores. These new scores do not require any additional conditions other than the original conditions made in <cit.>. The new DID estimators will be particularly appropriate when researchers would like to use ML methods in the first-stage nonparametric estimation. When using kernel estimators in the first-stage estimation , the new DID estimators do not require undersmoothing to achieve √(N)-consistency. Hence, researchers can use standard data-driven methods, such as CV, to select bandwidths. 

apacite






APPENDIX







Multilevel Treatment






Similarly, I use the cross-fitting algorithm <cit.>. 


  * Take a K-fold random partition (I_k)_k=1^K
of observation indices [N]={ 1,...,N}
such that the size of each fold I_k is n=N/K. For each k∈[K]={ 1,...,K},
define the auxiliary sample I_k^c{ 1,...,N}∖ I_k. 

  * For each k∈[K], construct the estimator of p_0
and λ_0 by p̂_w=1/n∑_i∈ I_k^cD_i.
Also, construct the estimators of g_w, g_z, and ℓ_30
using the auxiliary sample I_k^c: ĝ_wk=ĝ_w((W_i)_i∈ I_k^c),
ĝ_zk=ĝ_z((W_i)_i∈ I_k^c),
and ℓ̂_3k=ℓ̂_3((W_i)_i∈ I_k^c).

  * For each k, construct the intermediate ATT estimators 

    θ̃_wk=1/n∑_i∈ I_kI(W_i=w)·ĝ_zk(X_i)-I(W_i=0)·ĝ_wk(X_i)/p̂_wĝ_zk(X_i)×(Y(1)-Y(0)-ℓ̂_3k(X_i)),



  * Construct the final ATT estimators θ̃=1/K∑_k=1^Kθ̃_k.

The estimator ĝ_wk and ĝ_zk can be constructed
by Multi-Logit Lasso. 




The Lasso Penalty





The following is suggested by <cit.>. Let y_i denote Y_i(1)-Y_i(0) or (T_i-λ̂_k), λ_k denote λ_1k or λ_2k, and Υ̂_k denote Υ̂_1k or Υ̂_2k. For k∈[K], the loading Υ̂_k is a diagonal matrix with entries γ̂_kj, j=1,...,p, constructed by the following steps:


    Initial γ̂_kj=√(1/M_k∑_i∈ I_kz^cq_ij^2(y_i-y̅_k)^2),λ_k=2c√(M_k)Φ^-1(1-γ/2p),


    Refined γ̂_kj=√(1/M_k∑_i∈ I_kz^cq_ij^2ε̂_i^2),λ_k=2c√(M_k)Φ^-1(1-γ/2p),

where y̅_k=M^-1∑_i∈ I_k^cy_i, c>1 and γ→0. The empirical residual ε̂_i is calculated by the modified Lasso estimator β_k^* in the previous step: ε̂_i=y_i-q_i'β_k^*. Repeat the
second step B>0 times to obtain the final loading. 




PROOFS







  Proof of Lemma 1 

  Repeated outcomes: 

The Gateaux derivative of (3.1) in the direction η_1-η_10=(g-g_0,ℓ_1-ℓ_10)
is 

    ∂_η_1E_P[ψ_1(W,θ_0,p_0,η_10)](η_1-η_10)=     E_P[D-1/p_0(1-g_0(X))^2(Y(1)-Y(0)-ℓ_10)(g(X)-g_0(X))]
        -E_P[D-g_0(X)/p_0(1-g_0(X))(ℓ_1(X)-ℓ_10(X))]
    
    =     -E_P[g(X)-g_0(X)/p_0(1-g_0(X))E[Y(1)-Y(0)-ℓ_10(X)| X, D=0]]
        -E_P[(ℓ_1(X)-ℓ_10(X))/p_0(1-g_0(X))E_P[D-g_0(X)| X]]
    
    =     -E_P[g(X)-g_0(X)/p_0(1-g_0(X))(ℓ_10(X)-ℓ_10(X))]-0
    
    =     0,

where the second inequality follows from the law of iterated expectations, the third from the definition of ℓ_10(X) and E_P[D-g_0(X)| X]=0.




  Repeated cross sections: 

Define ∂_η_2E_P[ψ_20](η_2-η_20)∂_η_2E_P[ψ_2(W,θ_0,p_0,λ_0,η_20)](η_2-η_20). Similar to the proof of repeated outcomes,  the Gateaux derivative of (3.2) in the direction η_2-η_20=(g-g_0,ℓ_2-ℓ_20)
is

    ∂_η_2E_P[ψ_20](η_2-η_20)=     E_P[D-1/p_0'(1-g_0(X))^2[(T-λ_0)Y-ℓ_20(X)](g(X)-g_0(X))]
        -E_P[D-g_0(X)/p_0'(1-g_0(X))(ℓ_2(X)-ℓ_20(X))]
    
    =     -E_P[g(X)-g_0(X)/p_0'(1-g_0(X))(ℓ_20(X)-ℓ_20(X))]
        -E_P[ℓ_2(X)-ℓ_20(X)/pλ(1-λ)(1-g(X))E_P[D-g_0(X)| X]]
    
    =     0,

where p_0' p_0λ_0(1-λ_0).


  Multilevel treatment: 

Let Δ_w=g_w-g_w0, Δ_z=g_z-g_z0, and Δ_ℓ3=ℓ_3-ℓ_30. The Gateaux derivative of (3.3) in the direction η_w-η_w0=(g_w-g_w0,g_z-g_z0,ℓ_3-ℓ_30) is


    ∂_η_wE_P[ψ_w(W,θ_0,p_w0,η_w0)](η_w-η_w0)=   	E_P[I(W=0)g_w0(X)/p_w0g_z0(X)^2(Y(1)-Y(0)-ℓ_30)Δ_w]
       	-E_P[I(W=0)/p_w0g_z0(X)(Y(1)-Y(0)-ℓ_30)Δ_z]
       	+E_P[I(W=0)g_w0(X)-I(W=w)g_z0(X)/p_w0g_z0(X)Δ_ℓ3]
    
    =    	0

by the law of iterated expectation on each terms. 

The proofs of Theorem 1 and 2 follow the general framework proposed in <cit.>.


  Proof of Theorem 1 

  Repeated Outcomes: 

The proof proceeds in five steps. In Step 1, I show the main result
using the auxiliary results (A.1)-(A.4). In Step 2-5, I prove the
auxiliary results. 

    sup_η_1∈𝒯_N(E[∥ψ_1(W,θ_0,p_0,η_1)-ψ_1(W,θ_0,p_0,η_10)∥^2])^1/2≤ε_N,A.1



    sup_r∈(0,1),η_1∈𝒯_N∥∂_r^2E[ψ_1(W,θ_0,p_0,η_10+r(η_1-η_10))]∥≤(ε_N)^2,A.2



    sup_η_1∈𝒯_N(E_P[∥∂_pψ_1(W,θ_0,p_0,η_1)-∂_pψ_1(W,θ_0,p_0,η_10)∥^2])^1/2≤ε_N,A.3


    sup_p∈𝒫_N,η_1∈𝒯_N(E_P[∥∂_p^2ψ_1(W,θ_0,p,η_1)-∂_p^2ψ_1(W,θ_0,p_0,η_10)∥^2])^1/2≤ε_N,A.4

where  𝒯_N
is the set of all η_1=(g,ℓ_1) consisting of
P-square-integrable functions g and ℓ_1 such that 

    ∥η_1-η_10∥_P,2≤ε_N,


    ∥ g-1/2∥_P,∞≤1/2-κ,


    ∥ g-g_0∥_P,2^2+∥ g-g_0∥_P,2×∥ℓ_1-ℓ_10∥_P,2≤(ε_N)^2,

and 𝒫_N is the set of all p>0 such that | p-p_0|≤ N^-1/2. By the regularity condition (3.1) and |p̂_k-p_0|=O_P(N^-1/2), η̂_1k∈𝒯_N and p̂_k∈𝒫_N with probability 1-o(1).

Step 1. Observe that we have the decomposition


    √(N)(θ̃-θ_0)=    √(N)(1/K∑_k=1^Kθ̃_k-θ_0)
    
    =    √(N)1/K∑_k=1^K𝔼_n,k[ψ_1(W,θ_0,p̂_k,η̂_1k)]
    
    =    √(N)1/K∑_k=1^K𝔼_n,k[ψ_1(W,θ_0,p_0,η̂_1k)]
        +√(N)1/K∑_k=1^K𝔼_n,k[∂_pψ_1(W,θ_0,p_0,η̂_1k)](p̂_k-p_0)_a
        +√(N)1/K∑_k=1^K𝔼_n,k[∂_p^2ψ_1(W,θ_0,p̅_k,η̂_1k)](p̂_k-p_0)^2_b,

where p̅_k∈(p̂_k,p_0). For term (a),
by the triangle inequality, we have 

    |𝔼_n,k[∂_pψ_1(W,θ_0,p_0,η̂_1k)]-E_p[∂_pψ_1(W,θ_0,p_0,η_10)]|≤ J_1,k+J_2,k,

where 

    J_1,k=|𝔼_n,k[∂_pψ_1(W,θ_0,p_0,η̂_1k)]-𝔼_n,k[∂_pψ_1(W,θ_0,p_0,η_10)]|,


    J_2,k=|𝔼_n,k[∂_pψ_1(W,θ_0,p_0,η_10)]-E_p[∂_pψ_1(W,θ_0,p_0,η_10)]|.

To bound J_2,k, we have  

    E_P[J_2,k^2]≤    n^-1 E_P[( ∂_pψ_1(W,θ_0,p_0,η_10)^2)]
    
    =     n^-1E_P[1/p_0^4U^2V_1^2/(1-g_0)^2]
    ≤    n^-1(C^2/p_0^4κ^2),

where the last inequality follows from the regularity condition (3.1). By Chebyshev's inequality, J_2,k=O_P(n^-1/2)=o_P(1). Next, we bound J_1,k. Conditional on the auxiliary sample
I_k^c, η̂_1k can be treated as fixed. Under the
event that η̂_1k∈𝒯_N, we have

    E_P[J_1,k^2|(W_i)_i∈ I_k^c]=     E_P[∥∂_pψ_1(W,θ_0,p_0,η̂_1k)-∂_pψ_1(W,θ_0,p_0,η_10)∥^2|(W_i)_i∈ I_k^c]
    ≤   sup_η_1∈𝒯_NE_P[∥∂_pψ_1(W,θ_0,p_0,η_1)-∂_pψ_1(W,θ_0,p_0,η_10)∥^2|(W_i)_i∈ I_k^c]
    
    =    sup_η_1∈𝒯_NE_P[∥∂_pψ_1(W,θ_0,p_0,η_1)-∂_pψ_1(W,θ_0,p_0,η_10)∥^2]
    
    =    ε_N^2

by (A.3). By Lemma A.1, J_1,k=O_P(ε_N)=o_P(1).
Together, we have 

    𝔼_n,k[∂_pψ_1(W,θ_0,p_0,η̂_1k)]p→E_p[∂_pψ_1(W,θ_0,p_0,η_10)]=G_1p0.

For term (b), by the triangle inequality, we have

    |𝔼_n,k[∂_p^2ψ_1(W,θ_0,p̅_k,η̂_1k)]-E_p[∂_p^2ψ_1(W,θ_0,p_0,η_10)]|≤ J_3,k+J_4,k,

where 

    J_3,k=|𝔼_n,k[∂_p^2ψ_1(W,θ_0,p̅_k,η̂_1k)]-𝔼_n,k[∂_p^2ψ_1(W,θ_0,p_0,η_10)]|,


    J_4,k=|𝔼_n,k[∂_p^2ψ_1(W,θ_0,p_0,η_10)]-E_p[∂_p^2ψ_1(W,θ_0,p_0,η_10)]|.

To bound J_4,k, we have  

    E_P[J_4,k^2]≤    n^-1 E_P[( ∂_p^2ψ_1(W,θ_0,p_0,η_10)^2)]
    
    =     n^-1E_P[4/p_0^6U^2V_1^2/(1-g_0)^2]
    ≤    n^-1(4C^2/p_0^6κ^2),

where the last inequality follows from the regularity conditions. By Chebyshev's inequality, J_4,k=O_P(n^-1/2)=o_P(1).
Conditional on I_k^c, both p̅_k and η̂_1k
can be treated as fixed. Under the event that p̂_k∈𝒫_N
(thus p̅_k∈𝒫_N) and η̂_1k∈𝒯_N
, we have 

    E_P[J_3,k^2|(W_i)_i∈ I_k^c]=     E_P[∥∂_p^2ψ_1(W,θ_0,p̅_k,η̂_1k)-∂_p^2ψ_1(W,θ_0,p_0,η_10)∥^2|(W_i)_i∈ I_k^c]
    ≤   sup_p∈𝒫_N,η_1∈𝒯_NE_P[∥∂_pψ_1(W,θ_0,p,η_1)-∂_pψ_1(W,θ_0,p_0,η_10)∥^2|(W_i)_i∈ I_k^c]
    
    =    sup_p∈𝒫_N,η_1∈𝒯_NE_P[∥∂_pψ_1(W,θ_0,p,η_1)-∂_pψ_1(W,θ_0,p_0,η_10)∥^2]
    ≤   ε_N^2

by (A.4). By Lemma A.1, J_3,k=O_P(ε_N)=o_P(1).
Hence, 𝔼_n,k[∂_p^2ψ_1(W,θ_0,p̅_k,η̂_1k)]=O_P(1). 

Combine the above results with p̂_k-p_0=𝔼_n,k[D-p_0]
and (p̂_k-p_0)^2=O_P(N^-1),
the decomposition of θ̃ becomes

    √(N)(θ̃-θ_0)=    √(N)1/K∑_k=1^K𝔼_n,k[ψ_1(W,θ_0,p_0,η̂_1k)]
        +√(N)1/K∑_k=1^K𝔼_n,k[G_1p0(D-p_0)]+o_P(1)
    
    =    √(N)1/K∑_k=1^K𝔼_n,k[ψ_1(W,θ_0,p_0,η̂_1k)+G_1p0(D-p_0)]+o_P(1)
    
    =    1/√(N)∑_i=1^N[ψ_1(W_i,θ_0,p_0,η_10)+G_1p0(D_i-p_0)]+√(N)R_N+o_P(1),

where 

    R_N=    1/K∑_k=1^K𝔼_n,k[ψ_1(W,θ_0,p_0,η̂_1k)+G_1p0(D-p_0)]-1/N∑_i=1^N[ψ_1(W_i,θ_0,p_0,η_10)+G_1p0(D_i-p_0)]
    
    =    1/K∑_k=1^K𝔼_n,k[ψ_1(W,θ_0,p_0,η̂_1k)]-1/N∑_i=1^Nψ_1(W_i,θ_0,p_0,η_10).

If we can show that √(N)R_N=o_P(1), then we
are done.

This part is essentially identical to Step 3 in the proof of Theorem 3.1 (DML2) in <cit.>. I reproduce it here for reader's convenience. Since K is a fixed integer, which is independent of N, it suffices
to show that for any k∈[K],

    𝔼_n,k[ψ_1(W,θ_0,p_0,η̂_1k)]-1/n∑_i∈ I_kψ_1(W_i,θ_0,p_0,η_10)=o_P(N^-1/2).

Define the empirical process notation:

    𝔾_n,k[ϕ(W)]=1/√(n)∑_i∈ I_k(ϕ(W_i)-∫ϕ(w)dP),

where ϕ is any P-integrable function on 𝒲. By
the triangle inequality, we have

    ∥𝔼_n,k[ψ_1(W,θ_0,p_0,η̂_1k)]-1/n∑_i∈ I_kψ_1(W_i,θ_0,p_0,η_10)∥≤I_1,k+I_2,k/√(n),

where

    I_1,k∥𝔾_n,k[ψ_1(W,θ_0,p_0,η̂_1k)]-𝔾_n,k[ψ_1(W,θ_0,p_0,η_10)]∥,


    I_2,k√(n)∥ E_P[ψ_1(W,θ_0,p_0,η̂_1k)|(W_i)_i∈ I_k^c]-E_P[ψ_1(W,θ_0,p_0,η_10)]∥.

To bound I_1,k, note that conditional on (W_i)_i∈ I_k^c
the estimator η̂_1k is nonstochastic. Under the event
that η̂_1k∈𝒯_N, we have 

    E_P[I_1,k^2|(W_i)_i∈ I_k^c]=     E_P[∥ψ_1(W,θ_0,p_0,η̂_1k)-ψ_1(W,θ_0,p_0,η_10)∥^2|(W_i)_i∈ I_k^c]
    ≤   sup_η_1∈𝒯_NE_P[∥ψ_1(W,θ_0,p_0,η_1)-ψ_1(W,θ_0,p_0,η_10)∥^2|(W_i)_i∈ I_k^c]
    
    =    sup_η_1∈𝒯_NE_P[∥ψ_1(W,θ_0,p_0,η_1)-ψ_1(W,θ_0,p_0,η_10)∥^2]
    
    =    (ε_N)^2

by (A.1). Hence, I_1,k=O_P(ε_N) by Lemma
A.1. To bound I_2,k, define the following function

    f_k(r)=E_P[ψ_1(W,θ_0,p_0,η_10+r(η̂_1k-η_10))|(W_i)_i∈ I_k^c]-E[ψ_1(W,θ_0,p_0,η_10)],r∈[0,1).

By Taylor series expansion, we have

    f_k(1)=f_k(0)+f'_k(0)+f”_k(r̃)/2,for some r̃∈(0,1).

Note that f_k(0)=0 since E[ψ_1(W,θ_0,p_0,η_10)|(W_i)_i∈ I_k^c]=E[ψ_1(W,θ_0,p_0,η_10)].
Further, on the event η̂_1k∈𝒯_N,

    ∥ f_k'(0)∥=∥∂_η_1E_Pψ_1(W,θ_0,p_0,η_10)[η̂_1k-η_10]∥=0

by the orthogonality of ψ_1. Also, on the event η̂_1k∈𝒯_N,

    ∥ f_k”(r̃)∥≤sup_r∈(0,1)∥ f_k”(r)∥≤(ε_N)^2

by (A.2). Thus,

    I_2,k=√(n)∥ f_k(1)∥=O_P(√(n)(ε_N)^2).

Together with the result on I_1,k, we have 

    𝔼_n,k[ψ_1(W,θ_0,p_0,η̂_1k)]-1/n∑_i∈ I_kψ_1(W_i,θ_0,p_0,η_10)≤   I_1,k+I_2,k/√(n)
    
    =     O_P(n^-1/2ε_N+(ε_N)^2)
    
    =     o_P(N^-1/2)

by n=O(N) and ε_N=o(N^-1/4). Hence, √(N)R_N=o_P(1).

Step 2. In this step, I present the proof of (A.1).  We have the following decomposition: 

    ψ_1(W,θ_0,p_0,η_1)-ψ_1(W,θ_0,p_0,η_10)=    D-g(X)/p_0(1-g(X))(Y(1)-Y(0)-ℓ_1(X))
        -D-g_0(X)/p_0(1-g_0(X))(Y(1)-Y(0)-ℓ_10(X))
    
    =    U+g_0(X)-g(X)/p_0(1-g(X))(V_1+ℓ_10(X)-ℓ_1(X))
        -UV_1/p_0(1-g_0(X)).

Thus, we have

    ψ_1(W,θ_0,p_0,η_1)-ψ_1(W,θ_0,p_0,η_10)=    UV_1/p_0(1-g(X))+U(ℓ_10(X)-ℓ_1(X))/p_0(1-g(X))+(g_0(X)-g(X))V_1/p_0(1-g(X))
        +(g_0(X)-g(X))(ℓ_10(X)-ℓ_1(X))/p_0(1-g(X))-UV_1/p_0(1-g_0(X)).

Given κ≤ g_0(X)≤1-κ
and κ≤ g(X)≤1-κ, 

    ∥ψ_1(W,θ_0,p_0,η_1)-ψ_1(W,θ_0,p_0,η_10)∥_P,2≤   1/p_0κ^2∥ UV_1(1-g_0(X))
        +U(ℓ_10(X)-ℓ_1(X))(1-g_0(X))
        +V_1(g_0(X)-g(X))(1-g_0(X))
        +(g_0-g)(ℓ_10-ℓ_1)(1-g_0(X))
        -UV_1(1-g(X))∥_P,2.

By κ≤ g_0(X)≤1-κ and κ≤ g(X)≤1-κ
again, we can obtain

    ∥ψ_1(W,θ_0,p_0,η_1)-ψ_1(W,θ_0,p_0,η_10)∥_P,2≤   1-κ/p_0κ^2∥ UV_1+U(ℓ_10(X)-ℓ_1(X))
        +V_1(g_0(X)-g(X))
        +(g_0(X)-g(X))(ℓ_10(X)-ℓ_1(X))
        -UV_1∥_P,2.

Thus, by E_P[U^2| X]≤ C and E_P[V_1^2| X]≤ C,


    ∥ψ_1(W,θ_0,p_0,η_1)-ψ_1(W,θ_0,p_0,η_10)∥_P,2≤   (1-κ)√(C)/p_0κ^2∥ℓ_10-ℓ_1∥_P,2
        +(1-κ)√(C)/p_0κ^2∥ g_0-g∥_P,2
        +(1-κ)/p_0κ^2∥ g_0-g∥_P,2∥ℓ_10-ℓ_1∥_P,2
    ≤    O(ε_N+ε_N+(ε_N)^2)
    
    =     O(ε_N).


Step 3. In this step, I present the proof of (A.2). Define 
    f(r)=E_P[ψ_1(W,θ_0,p_0,η_10+r(η_1-η_10))].

Then its second-order derivative is

    ∂_r^2f(r)=    2/p_0E_P[(D-1)(g-g_0)^2/(1-g_0-r(g-g_0))^3(Y(1)-Y(0)-ℓ_10-r(ℓ_1-ℓ_10))]
        -2/p_0E_P[D-1/(1-g_0-r(g-g_0))^2(ℓ_1-ℓ_10)(g-g_0)].

It follows that 

    |∂_r^2f(r)|≤ O(∥(g-g_0)∥_P,2^2+∥(g-g_0)∥_P,2×∥(ℓ_1-ℓ_10)∥_P,2)≤(ε_N)^2.


Step 4. Notice that 

    ∂_pψ_1(W,θ,p,η_1)=     -1/pD-g(X)/1-g(X)(Y(1)-Y(0)-ℓ_1(X))
    
    =     -1/p(ψ_1(W,θ,p,η_1)+θ),

then we have

    ∥∂_pψ_1(W,θ_0,p_0,η_1)-∂_pψ_1(W,θ_0,p_0,η_10)∥_P,2=    1/p_0∥ψ_1(W,θ_0,p_0,η_1)-ψ_1(W,θ_0,p_0,η_10)∥_P,2
    
    =     O(ε_N)

by Step 2. 

Step 5. Notice that 

    ∂_p^2ψ_1(W,θ,p,η_1)=    2/p^3D-g(X)/1-g(X)(Y(1)-Y(0)-ℓ_1(X))
    
    =    2/p^2(ψ_1(W,θ,p,η_1)+θ),

then we have 

    ∂_p^2ψ_1(W,θ_0,p,η_1)-∂_p^2ψ_1(W,θ_0,p_0,η_10)=    ∂_p^2ψ_1(W,θ_0,p_0,η_1)-∂_p^2ψ_1(W,θ_0,p_0,η_10)
        +∂_p^3^3ψ_1(W,θ_0,p̅,η_1)(p-p_0)
    
    =    2/p_0^2(ψ_1(W,θ_0,p_0,η_1)-ψ_1(W,θ_0,p_0,η_10))
        -6/p̅^4D-g(X)/1-g(X)(Y(1)-Y(0)-ℓ_1(X))(p-p_0),

where p̅∈(p,p_0). Thus,

    ∥∂_p^2ψ_1(W,θ_0,p,η_1)-∂_p^2ψ_1(W,θ_0,p_0,η_10)∥_P,2≤   2/p_0^2∥ψ_1(W,θ_0,p_0,η_1)-ψ_1(W,θ_0,p_0,η_10)∥_P,2
        +∥6/p̅^4D-g(X)/1-g(X)(Y(1)-Y(0)-ℓ_1(X))∥_P,2
       ×| p-p_0|.

The term in the second line is bounded by 

    6/p̅^4κ∥(U+g_0-g)(V_1+ℓ_10-ℓ_1)∥_P,2≤   6/p̅^4κ∥ UV_1∥_P,2+6/p̅^4κ∥ U(ℓ_10-ℓ_1)∥_P,2
        +6/p̅^4κ∥ V_1(g_0-g)∥_P,2
        +6/p̅^4κ∥ g_0-g∥_P,2∥ℓ_10-ℓ_1∥_P,2
    ≤   6/p̅^4κ(C+√(C)∥ℓ_10-ℓ_1∥_P,2+√(C)∥ g_0-g∥_P,2)
        +6/p̅^4κ∥ g_0-g∥_P,2∥ℓ_10-ℓ_1∥_P,2
    
    =     O(1)

by ∥ UV_1∥_P,2≤∥ UV_1∥_P,4≤ C , E_P[U^2| X]≤ C,
E_P[V_1^2| X]≤ C, and the conditions on
the rates of convergence. Together with Step 2, we obtain

    ∥∂_p^2ψ_1(W,θ_0,p,η_1)-∂_p^2ψ_1(W,θ_0,p_0,η_10)∥_P,2≤    O(ε_N)+O(1)× O(N^-1/2)
    
    =     O(ε_N),

where I assume that ε_N converges to zero no faster
than N^-1/2. 


  Repeated cross sections: 

In step 1, I show the main result with the following auxiliary results:

    sup_η_2∈𝒯_N(E[∥ψ_2(W,θ_0,p_0,λ_0,η_2)-ψ_2(W,θ_0,p_0,λ_0,η_20)∥^2])^1/2≤ε_N,A.5



    sup_r∈(0,1),η_2∈𝒯_N∥∂_r^2E[ψ_2(W,θ_0,p_0,λ_0,η_20+r(η_2-η_20))]∥≤(ε_N)^2.A.6


    sup_η_2∈𝒯_N(E_P[∥∂_pψ_2(W,θ_0,p_0,λ_0,η_2)-∂_pψ_2(W,θ_0,p_0,λ_0,η_20)∥^2])^1/2≤ε_N,A.7


    sup_η_2∈𝒯_N(E_P[∥∂_λψ_2(W,θ_0,p_0,λ_0,η_2)-∂_λψ_2(W,θ_0,p_0,λ_0,η_20)∥^2])^1/2≤ε_N,A.8


    sup_p∈𝒫_N,η_2∈𝒯_N(E_P[∥∂_p^2ψ_2(W,θ_0,p,λ_0,η_2)-∂_p^2ψ_2(W,θ_0,p_0,λ_0,η_20)∥^2])^1/2≤ε_N,A.9


    sup_p∈𝒫_N,λ∈Λ_N,η_2∈𝒯_N(E_P[∥∂_λ^2ψ_2(W,θ_0,p,λ,η_2)-∂_λ^2ψ_2(W,θ_0,p_0,λ_0,η_20)∥^2])^1/2≤ε_N,A.10


    sup_p∈𝒫_N,η_2∈𝒯_N(E_P[∥∂_λ∂_pψ_2(W,θ_0,p,λ_0,η_2)-∂_λ∂_pψ_2(W,θ_0,p_0,λ_0,η_20)∥^2])^1/2≤ε_N,A.11

where 𝒯_N is the set of all η_2=(g,ℓ_2)
consisting of P-square-integrable functions g and ℓ_2
such that 

    ∥η_2-η_20∥_P,2≤ε_N,


    ∥ g-1/2∥_P,∞≤1/2-κ,


    ∥(g-g_0)∥_P,2^2+∥(g-g_0)∥_P,2×∥(ℓ_2-ℓ_20)∥_P,2≤(ε_N)^2,


𝒫_N and Λ_N are the sets consisting
all p>0 and λ>0 such that | p-p_0|≤ N^-1/2
and |λ-λ_0|≤ N^-1/2, respectively. By the regularity condition (3.2), |p̂_k-p_0|=O_P(N^-1/2), and |λ̂_k-λ_0|=O_P(N^-1/2),
we have η̂_2k∈𝒯_N, p̂_k∈𝒫_N, and λ̂_k∈Λ_N
with probability 1-o(1).  

In Step 2-4, I show the above auxiliary results.

Step 1. Notice that 

    √(N)(θ̃-θ_0)=    √(N)(1/K∑_k=1^Kθ̃_k-θ_0)
    
    =    √(N)1/K∑_k=1^K𝔼_n,k[ψ_2(W,θ_0,p̂_k,λ̂_k,η̂_2k)]
    
    =    √(N)1/K∑_k=1^K𝔼_n,k[ψ_2(W,θ_0,p_0,λ_0,η̂_2k)]
        +√(N)1/K∑_k=1^K𝔼_n,k[∂_pψ_2(W,θ_0,p_0,λ_0,η̂_2k)](p̂_k-p_0)
        +√(N)1/K∑_k=1^K𝔼_n,k[∂_λψ_2(W,θ_0,p_0,λ_0,η̂_2k)](λ̂_k-λ_0)+o_P(1),

where the term o_P(1), by the same arguments for the
term b in repeated outcomes and the auxiliary results (A.9)-(A.11),
contains the second-order terms 

    √(N)1/K∑_k=1^K𝔼_n,k[∂_p^2ψ_2(W,θ_0,p̅_k,λ_0,η̂_2k)](p̂_k-p_0)^2,


    √(N)1/K∑_k=1^K𝔼_n,k[∂_λ^2ψ_2(W,θ_0,p̂_k,λ̅_k,η̂_2k)](λ̂_k-λ_0)^2,


    √(N)1/K∑_k=1^K𝔼_n,k[∂_λ∂_pψ_2(W,θ_0,p̅_k,λ_0,η̂_2k)](λ̂_k-λ_0)(p̂_k-p_0),

where p̅_k∈(p̂_k,p_0) and λ̅_k∈(λ̂_k,λ_0). On the other hand, 
by the same arguments for the term a in repeated outcomes and the
auxiliary results (A.7)-(A.8), we have 

    𝔼_n,k[∂_pψ_2(W,θ_0,p_0,λ_0,η̂_2k)]p→E_p[∂_pψ_2(W,θ_0,p_0,λ_0,η_20)]=G_2p0,


    𝔼_n,k[∂_λψ_2(W,θ_0,p_0,λ_0,η̂_2k)]p→E_p[∂_λψ_2(W,θ_0,p_0,λ_0,η_20)]=G_2λ0.

Hence, since p̂_k-p_0=𝔼_n,k[D-p_0] and
λ̂_k-λ_0=𝔼_n,k[T-λ_0],
we have 

    √(N)(θ̃-θ_0)=    √(N)1/K∑_k=1^K𝔼_n,k[ψ_2(W,θ_0,p_0,λ_0,η̂_2k)]
    
    =    √(N)1/K∑_k=1^K𝔼_n,k[ψ_2(W,θ_0,p_0,λ_0,η̂_1k)+G_2p0(D-p_0)+G_2λ0(T-λ_0)]+o_P(1)
    
    =    1/√(N)∑_i=1^N[ψ_2(W_i,θ_0,p_0,λ_0,η_20)+G_2p0(D_i-p_0)+G_2λ0(T_i-λ_0)]
        +√(N)R_N'+o_P(1),

where 

    R_N'=    1/K∑_k=1^K𝔼_n,k[ψ_2(W,θ_0,p_0,λ_0,η̂_2k)+G_2p0(D-p_0)+G_2λ0(T-λ_0)]
        -1/N∑_i=1^N[ψ_2(W_i,θ_0,p_0,λ_0,η_20)+G_2p0(D_i-p_0)+G_2λ0(T_i-λ_0)]
    
    =    1/K∑_k=1^K𝔼_n,k[ψ_2(W,θ_0,p_0,λ_0,η̂_2k)]-1/N∑_i=1^Nψ_2(W_i,θ_0,p_0,λ_0,η_10).

Using (A.5)-(A.6) and the same arguments as the step 1 in repeated
outcomes, one can show that √(N)R_N'=o_P(1).
Hence, it remains to prove the auxiliary results (A.5)-(A.11). 

Step 2. 
Recall that p_0'= p_0λ_0(1-λ_0). For (A.5), notice that

    ψ_2(W,θ_0,p_0,λ_0,η_2)-ψ_2(W,θ_0,p_0,λ_0,η_20)=    D-g(X)/p_0'(1-g(X))((T-λ_0)Y-ℓ_2(X))
        -D-g_0(X)/p_0'(1-g_0(X))((T-λ_0)Y-ℓ_20(X))
    
    =    U+g_0(X)-g(X)/p_0'(1-g(X))(V_2+ℓ_20(X)-ℓ_2(X))
        -UV_2/p_0'(1-g_0(X)).

The decomposition becomes

    ψ_2(W,θ_0,p_0,λ_0,η_2)-ψ_2(W,θ_0,p_0,λ_0,η_20)=    UV_2/p_0'(1-g(X))+U(ℓ_20(X)-ℓ_2(X))/p_0'(1-g(X))
        +(g_0(X)-g(X))V_2/p_0'(1-g(X))
        +(g_0(X)-g(X))(ℓ_20(X)-ℓ_2(X))/p_0'(1-g(X))
        -UV_2/p_0'(1-g_0(X)).

Given that κ≤ g_0(X)≤1-κ,
κ≤ g(X)≤1-κ, we have

    ∥ψ_2(W,θ_0,p_0,λ_0,η_2)-ψ_2(W,θ_0,p_0,λ_0,η_20)∥_P,2≤   1/p_0'κ^2∥ UV_2(1-g_0(X))
        +U(ℓ_20(X)-ℓ_2(X))(1-g_0(X))
        +V_2(g_0(X)-g(X))(1-g_0(X))
        +(g_0-g)(ℓ_20-ℓ_2)(1-g_0(X))
        -UV_2(1-g(X))∥_P,2.

By κ≤ g_0(X)≤1-κ, κ≤ g(X)≤1-κ
again, we obtain

    ∥ψ_2(W,θ_0,p_0,λ_0,η_2)-ψ_2(W,θ_0,p_0,λ_0,η_20)∥_P,2≤   1-κ/p_0'κ^2∥ UV_2
        +U(ℓ_20(X)-ℓ_2(X))
        +V_2(g_0(X)-g(X))
        +(g_0-g)(ℓ_20-ℓ_2)
        -UV_2∥_P,2.

Given E_P[U^2| X]≤ C, E_P[V_2^2| X]≤ C,
and the conditions on the rates of convergence, 

    ∥ψ_2(W,θ_0,p_0,λ_0,η_2)-ψ_2(W,θ_0,p_0,λ_0,η_20)∥_P,2≤   (1-κ)√(C)/p_0'κ^2∥ℓ_20(X)-ℓ_2(X)∥_P,2
        +(1-κ)√(C)/p_0'κ^2∥ g_0(X)-g(X)∥_P,2
        +(1-κ)/p_0'κ^2∥ g_0-g∥_P,2∥ℓ_20-ℓ_2∥_P,2
    ≤    O(ε_N+ε_N+(ε_N)^2)
    
    =     O(ε_N).

For (A.6), let f(r)=E_P[ψ_2(W,θ_0,p_0,λ_0,η_20+r(η_2-η_20))].
Then the second-order derivative is

    ∂_r^2f(r)=    2/p_0'E_P[(D-1)(g-g_0)^2/(1-g_0-r(g-g_0))^3((T-λ_0)Y-ℓ_20-r(ℓ_2-ℓ_20))]
        -2/p_0'E_P[D-1/(1-g_0-r(g-g_0))^2(ℓ_2-ℓ_20)(g-g_0)]

It follows that 

    |∂_r^2f(r)|≤ O(∥(g-g_0)∥_P,2^2+∥(g-g_0)∥_P,2×∥(ℓ_2-ℓ_20)∥_P,2)≤(ε_N)^2.


Step 3. For (A.7), notice that 

    ∂_pψ_2(W,θ,p,λ,η_2)=     -1/p^2λ(1-λ)D-g(X)/1-g(X)((T-λ)Y-ℓ_2(X))
    
    =     -1/p(ψ_2(W,θ,p,λ,η_2)+θ),

then we have

    ∥∂_pψ_2(W,θ_0,p_0,λ_0,η_2)-∂_pψ_2(W,θ_0,p_0,λ_0,η_20)∥_P,2=    1/p_0∥ψ_2(W,θ_0,p_0,λ_0,η_2)
        -ψ_2(W,θ_0,p_0,λ_0,η_20)∥_P,2
    
    =     O(ε_N)

by the proof of (A.5).

For (A.8), notice that 

    ∂_λψ_2(W,θ,p,λ,η_2)=     -1-2λ/λ^2(1-λ)^2D-g(X)/p(1-g(X))((T-λ)Y-ℓ_2(X))
        -Y/pλ(1-λ)D-g(X)/1-g(X).

Define ∂_λψ_20∂_λψ_2(W,θ_0,p_0,λ_0,η_20),
then

    ∥∂_λψ_2(W,θ_0,p_0,λ_0,η_2)-∂_λψ_20∥_P,2=    ∥ψ_2(W,θ_0,p_0,λ_0,η_2)-ψ_2(W,θ_0,p_0,λ_0,η_20)∥_P,2
       ×|1-2λ_0|/λ_0(1-λ_0)+∥Y/p_0'(D-g(X)/1-g(X)-D-g_0(X)/1-g_0(X))∥_P,2
    
    =     O(ε_N)+∥Y/p_0'(D-g(X)/1-g(X)-D-g_0(X)/1-g_0(X))∥_P,2
    ≤    O(ε_N)+1/p_0'κ^2∥ Y(g-g_0)(D-1)∥_P,2
    ≤    O(ε_N)+√(C)/p_0'κ^2∥ g-g_0∥_P,2
    
    =     O(ε_N),

by (A.5) and E_P[Y^2| X]≤ C. 

Step 4. For (A.9), notice that we have 

    ∂_p^2ψ_2(W,θ,p,λ,η_2)=2/p^3λ(1-λ)D-g(X)/1-g(X)((T-λ)Y-ℓ_2(X)).

Define ∂_p^2ψ_20∂_p^2ψ_2(W,θ_0,p_0,λ_0,η_20),
then we have 

    ∂_p^2ψ_2(W,θ_0,p,λ_0,η_2)-∂_p^2ψ_20=    ∂_p^2ψ_2(W,θ_0,p_0,λ_0,η_2)-∂_p^2ψ_20
        +∂_p^3ψ_2(W,θ_0,p̅,λ_0,η_2)(p-p_0)
    
    =    2/p^2(ψ_2(W,θ_0,p_0,λ_0,η_2)-ψ_2(W,θ_0,p_0,λ_0,η_20))
        +∂_p^3ψ_2(W,θ_0,p̅,λ_0,η_2)(p-p_0),

where p̅∈(p,p_0). Hence, we have

    ∥∂_p^2ψ_2(W,θ_0,p,λ_0,η_2)-∂_p^2ψ_20∥_P,2≤   2/p^2∥ψ_2(W,θ_0,p_0,λ_0,η_2)-ψ_2(W,θ_0,p_0,λ_0,η_20)∥_P,2
        +∥D-g(X)/1-g(X)((T-λ_0)Y-ℓ_2(X))∥_P,2
       ×6/p̅^4λ_0(1-λ_0)| p-p_0|.

By (A.5), we have ∥ψ_2(W,θ_0,p_0,λ_0,η_2)-ψ_2(W,θ_0,p_0,λ_0,η_20)∥_P,2=O(ε_N).
The term in the second line is bounded by

    1/κ∥(U+g_0-g)(V_2+ℓ_20-ℓ_2)∥_P,2≤   1/κ∥ UV_2∥_P,2+1/κ∥ U(ℓ_20-ℓ_2)∥_P,2
        +1/κ∥ V_2(g_0-g)∥_P,2+1/κ∥ g_0-g∥_P,2∥ℓ_20-ℓ_2∥_P,2
    ≤   1/κ(C+√(C)∥ℓ_20-ℓ_2∥_P,2+√(C)∥ g_0-g∥_P,2)
        +1/κ∥ g_0-g∥_P,2∥ℓ_20-ℓ_2∥_P,2
    
    =     O(1)

by ∥ UV_2∥_P,2≤∥ UV_2∥_P,4≤ C, E_P[U^2| X]≤ C,
and E_P[V_2^2| X]≤ C. Thus, we obtain 

    ∥∂_p^2ψ_2(W,θ_0,p,λ_0,η_2)-∂_p^2ψ_20∥_P,2≤    O(ε_N)+O(1)× O(N^-1/2)
    
    =     O(ε_N),

where I assume that ε_N converges to zero no faster
than N^-1/2.

For (A.10), notice that we have

    ∂_λ^2ψ_2(W,θ,p,λ,η_2)=    c_1/pλ^3(1-λ)^3D-g(X)/1-g(X)((T-λ)Y-ℓ_2(X))
        +2-4λ/pλ^2(1-λ)^2D-g(X)/1-g(X)Y,

where c_1 is a constant depending on λ. Define ∂_λ^2ψ_20∂_λ^2ψ_2(W,θ_0,p_0,λ_0,η_20),
we have

    ∂_λ^2ψ_2(W,θ_0,p,λ,η_2)-∂_λ^2ψ_20=    ∂_λ^2ψ_2(W,θ_0,p_0,λ_0,η_2)-∂_λ^2ψ_20
        +∂_λ^2∂_pψ_2(W,θ_0,p̅,λ,η_2)(p-p_0)
        +∂_λ^3ψ_2(W,θ_0,p_0,λ̅,η_2)(λ-λ_0)
    
    =    c_1/λ_0^2(1-λ_0)^2(ψ_2(W,θ_0,p_0,λ_0,η_2)-ψ_2(W,θ_0,p_0,λ_0,η_20))
        +2-4λ_0/p_0λ_0^2(1-λ_0)^2(D-g(X)/1-g(X)-D-g_0(X)/1-g_0(X))Y
        +∂_λ^2∂_pψ_2(W,θ_0,p̅,λ,η_2)(p-p_0)
        +∂_λ^3ψ_2(W,θ_0,p_0,λ̅,η_2)(λ-λ_0),

where p̅∈(p,p_0) and λ̅∈(λ,λ_0).
By the triangle inequality, we have

    ∥∂_λ^2ψ_2(W,θ_0,p,λ,η_2)-∂_λ^2ψ_20∥_P,2≤   | c_1|/λ^2(1-λ)^2∥ψ_2(W,θ_0,p_0,λ_0,η_2)-ψ_2(W,θ_0,p_0,λ_0,η_20)∥_P,2
        +|2-4λ_0|/p_0λ_0^2(1-λ_0)^2∥(D-g(X)/1-g(X)-D-g_0(X)/1-g_0(X))Y∥_P,2
        +∥∂_λ^2∂_pψ_2(W,θ_0,p̅,λ,η_2)∥_P,2| p-p_0|
        +∥∂_λ^3ψ_2(W,θ_0,p_0,λ̅,η_2)∥_P,2|λ-λ_0|.

The norm term is the second line is bounded by 

    1/κ^2∥ Y(D-1)(g-g_0)∥_P,2≤   √(C)/κ^2∥ g-g_0∥_P,2
    
    =     O(ε_N),

by E_P[Y^2| X]≤ C and D∈{ 0,1}.
The two high-order terms are bounded by

    ∥∂_λ^2∂_pψ_2(W,θ_0,p̅,λ,η_2)∥_P,2≤   | c_1|/p̅^2λ^3(1-λ)^3∥D-g(X)/1-g(X)((T-λ)Y-ℓ_2(X))∥_P,2
        +|2-4λ|/p̅λ^2(1-λ)^2∥D-g(X)/1-g(X)Y∥_P,2.

and

    ∥∂_λ^3ψ_2(W,θ_0,p_0,λ̅,η_2)∥_P,2≤   | c_2|/p_0λ̅^4(1-λ̅)^4∥D-g(X)/1-g(X)((T-λ̅)Y-ℓ_2(X))∥_P,2
        +| c_3|/p_0λ̅^3(1-λ̅)^3∥D-g(X)/1-g(X)× Y∥_P,2,

where c_2 and c_3 are constants depending on λ.
Using the same arguments in (A.9), one can show that 

    ∥D-g(X)/1-g(X)((T-λ)Y-ℓ_2(X))∥_P,2≤ O(1),


    ∥D-g(X)/1-g(X)((T-λ̅)Y-ℓ_2(X))∥_P,2≤ O(1).

Also, we have 

    ∥D-g(X)/1-g(X)× Y∥_P,2=    ∥U+g_0(X)-g(X)/1-g(X)× Y∥_P,2
    ≤   1/κ(∥ UY∥_P,2+∥(g_0-g)Y∥_P,2)
    ≤   1/κ(C+√(C)∥ g_0-g∥_P,2)
    
    =     O(1)

by ∥ UY∥_P,2≤ C and E_P[Y^2| X]≤ C. 

Finally, we obtain

    ∥∂_λ^2ψ_2(W,θ_0,p,λ,η_2)-∂_λ^2ψ_20∥_P,2≤    O(ε_N)+O(ε_N)+O(1)O(N^-1/2)+O(1)O(N^-1/2)
    
    =     O(ε_N),

where I assume that ε_N converges to zero no faster
than N^-1/2.

For (A.11), notice that the derivative is

    ∂_λ∂_pψ_2(W,θ,p,λ,η_2)=    1-2λ/p^2λ^2(1-λ)^2D-g(X)/1-g(X)((T-λ)Y-ℓ_2(X))
        +Y/p^2λ(1-λ)D-g(X)/1-g(X).

Define ∂_λ∂_pψ_20∂_λ∂_pψ_2(W,θ_0,p_0,λ_0,η_20),
then we have

    ∂_λ∂_pψ_2(W,θ_0,p,λ_0,η_2)-∂_λ∂_pψ_20=    ∂_λ∂_pψ_2(W,θ_0,p_0,λ_0,η_2)-∂_λ∂_pψ_20
        +∂_λ∂_p^2ψ_2(W,θ_0,p̅,λ_0,η_2)(p-p_0),

where p̅∈(p,p_0). By the triangle inequality,
we obtain

    ∥∂_λ∂_pψ_2(W,θ_0,p,λ_0,η_2)-∂_λ∂_pψ_20∥_P,2≤   1/p∥∂_λψ_2(W,θ_0,p_0,λ_0,η_2)-∂_λψ_2(W,θ_0,p_0,λ_0,η_20)∥_P,2
        +∥∂_λ∂_p^2ψ_2(W,θ_0,p̅,λ_0,η_2)∥_P,2| p-p_0|.

Using the same arguments in (A.9) and (A.10), one can show that the high-order
term is bounded by 

    ∥∂_λ∂_p^2ψ_2(W,θ_0,p̅,λ_0,η_2)∥_P,2≤   ∥2-4λ_0/p̅^3λ_0^2(1-λ_0)^2D-g(X)/1-g(X)((T-λ_0)Y-ℓ_2(X))∥_P,2
        +∥2Y/p̅^3λ_0(1-λ_0)D-g(X)/1-g(X)∥_P,2
    ≤    O(1).

Together with (A.8), we obtain 

    ∥∂_λ∂_pψ_2(W,θ_0,p,λ_0,η_2)-∂_λ∂_pψ_20∥_P,2≤    O(ε_N)+O(1)O(N^-1/2)
    
    =     O(ε_N),

where I assume that ε_N converges to zero no faster
than N^-1/2.


  Proof of Theorem 2 

  Repeated outcomes: 

In Step 1, I show the main result using the auxiliary results 

    sup_p∈𝒫_N,η_1∈𝒯_N(E_P[∥ψ̅_1(W,θ_0,p,η_1)-ψ̅_1(W,θ_0,p_0,η_10)∥^2])^1/2≤ε_N,A.12



    (E_P[ψ̅_1(W,θ_0,p_0,η_10)^4])^1/4≤ C_1, A.13

where 𝒫_N and 𝒯_N are specified in the proof of Theorem 1, C_1 is a constant, and

    ψ̅_1(W,θ,p,η_1)1/pD-g(X)/1-g(X)(Y(1)-Y(0)-ℓ_1(X))-Dθ/p.

In fact, we have E_P[(ψ̅_1(W,θ_0,p_0,η_10))^2]=Σ_10.
In Step 2, I show the auxiliary results (A.12) and (A.13).

Step 1. Notice that 

    Σ̂_1=    1/K∑_k=1^K𝔼_n,k[(ψ_1(W,θ̃,p̂_k,η̂_1k)+Ĝ_1p(D-p̂_k))^2]
    
    =    1/K∑_k=1^K𝔼_n,k[(1/p̂_kD-ĝ_k(X)/1-ĝ_k(X)(Y(1)-Y(0)-ℓ̂_1k(X))-Dθ̃/p̂_k)^2]
    
    =    1/K∑_k=1^K𝔼_n,k[ψ̅_1(W,θ̃,p̂_k,η̂_1k)^2],

where the second equality follows from Ĝ_1p=-θ̃/p̂_k.

Since K is fixed, which is independent of N, it suffices to
show that for each k∈[k],

    I_k|𝔼_n,k[ψ̅_1(W,θ̃,p̂_k,η̂_1k)^2]-E_P[ψ̅_1(W,θ_0,p_0,η_10)^2]|=o_P(1).

By the triangle inequality, we have

    I_k≤ I_3,k+I_4,k,

where

    I_3,k|𝔼_n,k[ψ̅_1(W,θ̃,p̂_k,η̂_1k)^2]-𝔼_n,k[ψ̅_1(W,θ_0,p_0,η_10)^2]|,


    I_4,k|𝔼_n,k[ψ̅_1(W,θ_0,p_0,η_10)^2]-E_P[ψ̅_1(W,θ_0,p_0,η_10)^2]|.

To bound I_4,k, we have  

    E_P[I_4,k^2]≤    n^-1 E_P[ ψ̅_1(W,θ_0,p_0,η_10)^4]
    ≤    n^-1C_1^4,

where the last inequality follows from (A.13). By Chebyshev's inequality, I_4,k=O_P(n^1/2). 

Next, we bound I_3,k. This part is essentially identical to the proof of Theorem 3.2 in <cit.>, I reproduce it here for reader's convenience. Observe that for any number a and δ a,

    |(a+δ a)^2-a^2|≤2(δ a)(a+δ a).

Denote ψ_i=ψ̅_1(W_i,θ_0,p_0,η_10)
and ψ̂_i=ψ̅_1(W_i,θ̃,p̂_k,η̂_1k),
and aψ_i, a+δ aψ̂_i. Then

    I_3,k=    |1/n∑_i∈ I_k(ψ̂_i)^2-(ψ_i)^2|≤1/n∑_i∈ I_k|(ψ̂_i)^2-(ψ_i)^2|
    ≤   2/n∑_i∈ I_k|ψ̂_i-ψ_i|×(|ψ_i|+|ψ̂_i-ψ_i|)
    ≤   (2/n∑_i∈ I_k|ψ̂_i-ψ_i|^2)^1/2(2/n∑_i∈ I_k(|ψ_i|+|ψ̂_i-ψ_i|)^2)^1/2
    ≤   (2/n∑_i∈ I_k|ψ̂_i-ψ_i|^2)^1/2[(2/n∑_i∈ I_k|ψ_i|^2)^1/2+(2/n∑_i∈ I_k|ψ̂_i-ψ_i|^2)^1/2].

Thus,

    I_3,k^2≲ S_N×(1/n∑_i∈ I_k∥ψ̅_1(W_i,θ_0,p_0,η_10)∥^2+S_N),

where

    S_N1/n∑_i∈ I_k∥ψ̅_1(W_i,θ̃,p̂_k,η̂_1k)-ψ̅_1(W_i,θ_0,p_0,η_10)∥^2.

Since 1/n∑_i∈ I_k∥ψ̅_̅1̅(W_i,θ_0,p_0,η_0)∥^2=O_P(1),
it suffices to bound S_N. We have the decomposition

    S_N=    1/n∑_i∈ I_k∥ψ̅_1(W_i,θ_0,p̂_k,η̂_1k)+∂_θψ̅_1(W_i,θ̅,p̂_k,η̂_1k)(θ̃-θ_0)-ψ̅_1(W_i,θ_0,p_0,η_10)∥^2
    
    =    1/n∑_i∈ I_k∥ψ̅_1(W_i,θ_0,p̂_k,η̂_1k)+D_i/p̂_k(θ̃-θ_0)-ψ̅_1(W_i,θ_0,p_0,η_10)∥^2
    ≤   1/n∑_i∈ I_k∥D_i/p̂_k(θ̃-θ_0)∥^2+1/n∑_i∈ I_k∥ψ̅_1(W_i,θ_0,p̂_k,η̂_1k)-ψ̅_1(W_i,θ_0,p_0,η_10)∥^2,

where θ̅∈(θ̃-θ_0). The 
first term is bounded by

    1/n∑_i∈ I_k∥D_i/p̂_k(θ̃-θ_0)∥^2≤   (1/n∑_i∈ I_k(D_i/p̂_k)^2)∥θ̃-θ_0∥^2
    
    =    (1/n∑_i∈ I_k(D_i/p_0)^2+o_P(1))∥θ̃-θ_0∥^2
    
    =     O_P(1)× O_P(N^-1).

Also, notice that conditional on (W_i)_i∈ I_k^c,
both p̂_k and η̂_1k can be treated as fixed.
Under the event that p̂_k∈𝒫_N and η̂_1k∈𝒯_N,
we have

    E_P[∥ψ̅_1(W_i,θ_0,p̂_k,η̂_1k)-ψ̅_1(W_i,θ_0,p_0,η_10)∥^2|(W_i)_i∈ I_k^c]
    ≤sup_p∈𝒫_N,η_1∈𝒯_NE_P[∥ψ̅_1(W_i,θ_0,p,η_1)-ψ̅_1(W_i,θ_0,p_0,η_10)∥^2]     =(ε_N)^2

by (A.12). It follows that S_N=O_P(N^-1+(ε_N)^2).
Therefore, we obtain

    I_k=O_P(N^-1/2)+O_P(N^-1/2+ε_N)=o_P(1).


Step 2. It remains to prove (A.12) and (A.13). By Taylor series expansion, 

    ψ̅_1(W,θ_0,p,η_1)-ψ̅_1(W,θ_0,p_0,η_10)=    ψ̅_1(W,θ_0,p_0,η_1)-ψ̅_1(W,θ_0,p_0,η_10)
        +∂_pψ_1(W,θ_0,p̅,η_1)(p-p_0)
    
    =    ψ_1(W,θ_0,p_0,η_1)-ψ_1(W,θ_0,p_0,η_10)
        +∂_pψ_1(W,θ_0,p̅,η_1)(p-p_0),

where p̅∈(p,p_0). Then we have 

    ∥ψ̅_1(W,θ_0,p,η_1)-ψ̅_1(W,θ_0,p_0,η_10)∥_P,2≤   ∥ψ_1(W,θ_0,p_0,η_1)-ψ_1(W,θ_0,p_0,η_10)∥_P,2
        +∥1/p̅^2D-g(X)/1-g(X)(Y(1)-Y(0)-ℓ_1(X))+Dθ_0/p̅^2∥_P,2
       ×| p-p_0|.

By (A.1), we have ∥ψ_1(W,θ_0,p_0,η_1)-ψ_1(W,θ_0,p_0,η_10)∥_P,2=O(ε_N).
The term in the second line is bounded by 

    ∥1/p̅^2U+g_0-g/1-g(U+ℓ_10-ℓ_1)∥_P,2+∥Dθ_0/p̅^2∥_P,2≤   1/p̅^2κ∥ UV_1∥_P,2+1/p̅^2κ∥ U(ℓ_10-ℓ_1)∥_P,2
        +1/p̅^2κ∥ V_1(g_0-g)∥_P,2+1/p̅^2|θ_0|
        +1/p̅^2κ∥ g_0-g_1∥_P,2∥ℓ_10-ℓ_1∥_P,2
    ≤   1/p̅^2κ(C+√(C)∥ℓ_10-ℓ_1∥_P,2+√(C)∥ g_0-g∥_P,2)
        +C/p̅^2p_0κ+1/p̅^2κ∥ g_0-g_1∥_P,2∥ℓ_10-ℓ_1∥_P,2
    
    =     O(1),

where I use ∥ UV_1∥_P,2≤∥ UV_1∥_P,4≤ C , E_P[U^2| X]≤ C,
E_P[V_1^2| X]≤ C, and

    |θ_0|=    | E_P[Y(1)-Y(0)/p_0D-g_0(X)/1-g_0(X)]|
    ≤   1/p_0κ| E_P[(Y(1)-Y(0))U]|
    
    =    1/p_0κ| E_P[(ℓ_10(X)+V_1)U]|
    
    =    1/p_0κ| E_P[UV_1]|
    ≤   C/p_0κ

by | E_P[UV_1]|≤∥ UV_1∥_P,4≤ C.
Thus, we obtain 

    ∥ψ̅_1(W,θ_0,p,η_1)-ψ̅_1(W,θ_0,p_0,η_10)∥_P,2≤    O(ε_N)+O(1)O(N^-1/2)
    
    =     O(ε_N),

where I assume that ε_N converges to zero no faster
than N^-1/2. 

For (A.13), 

    ∥ψ̅_1(W, θ_0, p_0, η_10)∥_P,4=   ∥1/p_0UV_1/1-g_0-Dθ_0/p_0∥_P,4
    ≤   ∥1/p_0UV_1/1-g_0∥_P,4+ ∥Dθ_0/p_0∥_P,4
    ≤   1/p_0κ∥ UV_1∥_P,4+1/p_0|θ_0|
    ≤   C/p_0κ+C/p_0^2κ

since ∥ UV_1∥_P,4≤ C.


  Repeated cross sections: 

In Step 1, I show the main result with the auxiliary results:

    sup_p∈𝒫_N,λ∈Λ_N,η_2∈𝒯_N(E_P[∥ψ̅_2(W,θ_0,p,λ,G_2λ0,η_2)-ψ̅_2(W,θ_0,p_0,λ_0,G_2λ0,η_20)∥^2])^2≤ε_N,A.14



    (E_P[ψ̅_2(W,θ_0,p_0,λ_0,G_2λ0,η_20)^4])^1/4≤ C_2, A.15

where (P_N, Λ_N, 𝒯_N) are specified in the proof of Theorem 1, C_2 is a constant, and

    ψ̅_2(W,θ,p,λ,G_2λ,η_2)1/λ(1-λ)pD-g(X)/1-g(X)((T-λ)Y-ℓ_2(X))-Dθ/p+G_2λ(T-λ).

In fact, we have E_P[(ψ̅_2(W,θ_0,p_0,λ_0,G_2λ_0,η_20))^2]=Σ_20. In Step 2, I prove (A.14) and (A.15).

Step 1. Notice that 

    Σ̂_2=    1/K∑_k=1^K𝔼_n,k[(ψ_2(W,θ̃,p̂_k,η̂_1k)+Ĝ_2p(D-p̂_k)+Ĝ_2λ(T-λ̂_k))^2]
    
    =    1/K∑_k=1^K𝔼_n,k[(1/λ̂_k(1-λ̂_k)p̂_kD-ĝ_k(X)/1-ĝ_k(X)((T-λ̂_k)Y-ℓ̂_2k(X))-Dθ̃/p̂_k+Ĝ_2λ(T-λ̂_k))^2]
    
    =    1/K∑_k=1^K𝔼_n,k[ψ̅_2(W,θ̃,p̂_k,λ̂_k,Ĝ_2λ,η̂_2k)^2],

where the second inequality follows from Ĝ_2p=-θ̃/p̂_k
.

Since K is fixed, which is independent of N, it suffices to
show that 

    J_k|𝔼_n,k[ψ̅_2(W,θ̃,p̂_k,λ̂_k,Ĝ_2λ,η̂_2k)^2]-E_P[ψ̅_2(W,θ_0,p_0,λ_0,G_2λ0,η_20)^2]|=o_P(1).

By the triangle inequality, we have

    J_k≤ J_5,k+J_6,k,

where 

    J_5,k|𝔼_n,k[ψ̅_2(W,θ̃,p̂_k,λ̂_k,Ĝ_2λ,η̂_2k)^2]-𝔼_n,k[ψ̅_2(W,θ_0,p_0,λ_0,G_2λ0,η_20)^2]|,



    J_6,k|𝔼_n,k[ψ̅_2(W,θ_0,p_0,λ_0,G_2λ0,η_20)^2]-E_P[ψ̅_2(W,θ_0,p_0,λ_0,G_2λ0,η_20)^2]|.

By the same arguments for I_4,k in the proof of repeated outcomes and (A.15), we can show J_6,k=o_P(1). Also,
by the same arguments for I_3,k in the proof of repeated outcomes, we have 

    J_5,k^2≲ S_N'×(1/n∑_i∈ I_k∥ψ̅_̅2̅(W,θ_0,p_0,λ_0,G_2λ0,η_20)∥^2+S_N'),

where

    S_N'1/n∑_i∈ I_k∥ψ̅_̅2̅(W,θ̃,p̂_k,λ̂_k,Ĝ_2λ,η̂_2k)-ψ̅_̅2̅(W,θ_0,p_0,λ_0,G_2λ0,η_20)∥^2.

Since 1/n∑_i∈ I_k∥ψ̅_̅2̅(W,θ_0,p_0,λ_0,G_2λ0,η_20)∥^2=O_P(1),
it remains to bound S_N'. Define ψ̅_20ψ̅_̅2̅(W,θ_0,p_0,λ_0,G_2λ0,η_20).
By the triangle inequality, we have 

    S_N'≤   1/n∑_i∈ I_k∥ψ̅_2(W_i,θ_0,p̂_k,λ̂_k,G_2λ0,η̂_2k)-ψ̅_20∥^2
        +1/n∑_i∈ I_k∥∂_θψ̅_2(W_i,θ̅,p̂_k,λ̂_k,Ĝ_2λ,η̂_2k)(θ̃-θ_0)∥^2
        +1/n∑_i∈ I_k∥∂_G_2λψ̅_2(W_i,θ_0,p̂_k,λ̂_k,G_2λ,η̂_2k)(Ĝ_2λ-G_2λ0)∥^2,

where θ̅∈(θ̃,θ_0) and
G_2λ∈(G_2λ,G_2λ0).
Then we have 

    S_N'≤   1/n∑_i∈ I_k∥ψ̅_2(W_i,θ_0,p̂_k,λ̂_k,G_2λ0,η̂_2k)-ψ̅_20∥^2
        +1/n∑_i∈ I_k∥D_i/p̂_k(θ̃-θ_0)∥^2+1/n∑_i∈ I_k∥(T_i-λ̂_k)(Ĝ_2λ-G_2λ0)∥^2.

 The two terms in the last line are bounded by 

    (1/n∑_i∈ I_k(D_i/p̂_k)^2)×∥θ̃-θ_0∥^2=O_P(1)× O_P(N^-1)

and

    (1/n∑_i∈ I_k(T_i-λ̂_k)^2)×∥Ĝ_2λ-G_2λ0∥^2=O_P(1)×(o_P(1))^2.

Conditional on the auxiliary sample I_k^c, (p̂_k,λ̂_k,η̂_2k)
can be treated as fixed. Also, under the event that p̂_k∈𝒫_N,
λ̂_k∈Λ_N, and η̂_2k∈𝒯_N,
we have

    E_P[∥ψ̅_2(W_i,θ_0,p̂_k,λ̂_k,G_2λ0,η̂_2k)-ψ̅_2(W_i,θ_0,p_0,λ_0,G_2λ0,η_20)∥^2|(W_i)_i∈ I_k^c]
    ≤sup_p∈𝒫_N,λ∈Λ_N,η_2∈𝒯_NE_P[∥ψ̅_2(W_i,θ_0,p,λ,G_2λ0,η_2)-ψ̅_2(W_i,θ_0,p_0,λ_0,G_2λ0,η_20)∥^2]     =(ε_N)^2

by (A.14). It follows that S_N'=O_P(N^-1+ε_N^2)+(o_P(1))^2
so that 

    J_k=o_P(1)+O_P(N^-1/2+ε_N)+o_P(1)=o_P(1).


Step 2. It remains to show (A.14) and (A.15). Define ψ̅_20ψ̅_2(W,θ_0,p_0,λ_0,G_2λ0,η_20). By the triangle inequality
and 

    ψ̅_2(W,θ_0,p_0,λ_0,G_2λ0,η_2)-ψ̅_20=ψ_2(W,θ_0,p_0,λ_0,η_2)-ψ_2(W,θ_0,p_0,λ_0,η_20),

we have 

    ∥ψ̅_2(W,θ_0,p,λ,G_2λ0,η_2)-ψ̅_20∥_P,2≤   ∥ψ_2(W,θ_0,p_0,λ_0,η_2)-ψ_2(W,θ_0,p_0,λ_0,η_20)∥_P,2
        +∥∂_λψ̅_2(W_i,θ_0,p_0,λ̅,G_2λ0,η_2)∥_P,2|λ-λ_0|
        +∥∂_pψ̅_2(W_i,θ_0,p̅,λ,G_2λ0,η_2)∥_P,2| p-p_0|,

where p̅∈(p,p_0) and λ̅∈(λ,λ_0).
The term in the second line is bounded by 

    ∥∂_λψ̅_2(W_i,θ_0,p_0,λ̅,G_2λ0,η_2)∥_P,2≤   |1-2λ̅|/p_0λ̅^2(1-λ̅)^2∥D-g(X)/1-g(X)((T-λ̅)Y-ℓ_2(X))∥_P,2
        +1/p_0λ̅(1-λ̅)∥D-g(X)/1-g(X)× Y∥_P,2+| G_2λ0|
    ≤    O(1)

by the same arguments in (A.9)-(A.11) and 

    | G_2λ0|=    | E_P[-1-2λ_0/λ_0^2(1-λ_0)^2p_0D-g_0/1-g_0((T-λ_0)Y-ℓ_20)-Y/λ_0(1-λ_0)p_0D-g_0/1-g_0]|
    ≤   |1-2λ_0|/λ_0^2(1-λ_0)^2p_0κ| E_P[UV_2]|+1/λ_0(1-λ_0)p_0κ| E_P[YU]|
    ≤   |1-2λ_0|/λ_0^2(1-λ_0)^2p_0κC+1/λ_0(1-λ_0)p_0κC
    
    =     O(1)

since | E_P[UV_2]|≤∥ UV_2∥_P,4≤ C
and | E_P[YU]|≤ C. Also, we have 

    ∥∂_pψ̅_2(W_i,θ_0,p̅,λ,G_2λ0,η_2)∥_P,2≤   1/λ(1-λ)p̅^2∥D-g(X)/1-g(X)((T-λ)Y-ℓ_2(X))∥_P,2
        +∥Dθ_0/p̅^2∥_P,2
    ≤    O(1)

by the same arguments in (A.9)-(A.11) and

    |θ_0|=    | E_P[D-g_0(X)/p_0'(1-g_0(X))(T-λ_0)Y]|
    ≤   1/p_0'κ| E_P[(T-λ_0)YU]|
    
    =    1/p_0κ| E_P[(ℓ_20(X)+V_2)U]|
    
    =    1/p_0κ| E_P[UV_2]|
    ≤   C/p_0κ

since | E_P[UV_2]|≤∥ UV_2∥_P,4≤ C.
Together with (A.5), we have 

    ∥ψ̅_2(W,θ_0,p,λ,G_2λ0,η_2)-ψ̅_20∥_P,2≤    O(ε_N)+O(1)O(N^-1/2)+O(1)O(N^-1/2)
    
    =     O(ε_N),

where I assume that ε_N converges to zero no faster
than N^-1/2. 

For (A.15), we have 

    ∥ψ̅_2(W,θ_0, p_0, λ_0, G_2λ0, η_20)∥_P,4=   ∥1/λ_0(1-λ_0)p_0UV_2/1-g_0-Dθ_0/p_0+G_2λ0(T-λ_0)∥_P,4
    ≤   1/λ_0(1-λ_0)p_0κ∥ UV_2∥_P,4+ 1/p_0|θ_0| +| G_2λ0|
    ≤    O(1)

since ∥ UV_2∥_P,4≤ C.


  Proof of Theorem 3: 
 
By the same arguments in the proof of Theorem 1, we can have 

    √(N)(θ̃-θ_0)=1/√(N)∑_i=1^Nψ_1(W_i,θ_0,p_0,η_10)+G_1p0(D-p_0)+O_P(ε_N+√(N)(ε_N)^2)

for repeated outcomes and

    √(N)(θ̃-θ_0)=    1/√(N)∑_i=1^Nψ_2(W_i,θ_0,p_0,λ_0,η_20)+G_2p0(D-p_0)+G_2λ0(T-λ_0)
        +O_P(ε_N+√(N)(ε_N)^2)

for repeated cross sections. The term ε_N is the rate of convergence of the kernel
estimators ĝ_kh,  ℓ̂_1kh, and ℓ̂_2kh.
It remains to show ∥ĝ_kh-g_0∥_P,2=o_P(N^-1/4),
∥ℓ̂_1kh-ℓ_10∥_P,2=o_P(N^-1/4), and
∥ℓ̂_2kh-ℓ_20∥_P,2=o_P(N^-1/4). 


Here I use the standard result of kernel estimation in <cit.>. Let γ̂_kh(x) denote the kernel estimator of γ_0(x)=f_0(x)E_P[z| x] using the auxiliary sample I_k^c, where z∈{1,D,Y(1)-Y(0)| D=0,(T-λ_0)Y| D=0}. By Assumption (3.3) and Lemma 8.10 of <cit.>, we have 


    sup_x∈𝒳|γ̂_kh(x)-γ_0(x)|=    O_P((log N)^1/2(Nh^d+2s)^-1/2+h^m)
    
    =     o_P(N^-1/4)

by the conditions on h. Let f̂_kh(x) denote γ̂_k(x) with z=1 and m̂_kh(x) denote γ̂_k(x) with z=D. Then

    ĝ_kh(x)=m̂_kh(x)/f̂_kh(x)=m̂_kh(x)/f_0(x)/f̂_kh(x)/f_0(x).
 
For the denominator, we have

    sup_x∈𝒳|f̂_kh(x)/f_0(x)-1 |≤sup_x∈𝒳|f̂_kh(x)-f_0(x) |/inf_x∈𝒳f_0(x)=o_P(N^-1/4)

given inf_x∈𝒳f_0(x)≠ 0. For the numerator, let m_0(x) denote γ_0(x) with z=D, we have 

    sup_x∈𝒳|m̂_kh(x)/f_0(x)-g_0(x) |≤sup_x∈𝒳|m̂_kh(x)-m_0(x) |/inf_x∈𝒳f_0(x)=o_P(N^-1/4)

given inf_x∈𝒳f_0(x)≠ 0. The above two inequalities imply that uniformly over x∈𝒳,  


    ĝ_kh(x)=g_0(x)+o_P(N^-1/4)/1+o_P(N^-1/4)= g_0(x)+o_P(N^-1/4).

That is, sup_x∈𝒳|ĝ_kh(x)-g_0(x)|=o_P(N^-1/4). Using the same arguments, one can also show that sup_x∈𝒳|ℓ̂_1kh(x)-ℓ_0(x)|=o_P(N^-1/4) and sup_x∈𝒳|ℓ̂_2kh(x)-ℓ_0(x)|=o_P(N^-1/4). Since uniform convergence implies L^2-norm convergence, we complete the proof. 


  Proof of Theorem 4: 

The proof is the same as the proof in Theorem 2 provided the assumptions in Theorem 3 hold.  


  Lemma A.1 (CONDITIONAL CONVERGENCE IMPLIES UNCONDITIONAL) 

Let
{ X_m} and { Y_m} be sequences
of random vectors. (i) If for ϵ_m→0, Pr(∥ X_m∥>ϵ_m| Y_m)p→0,
then Pr(∥ X_m∥>ϵ_m)→0.
This occurs if E[∥ X_m∥^q/ϵ_m^q| Y_m]p→0
for some q≥1, by Markov's inequality. (ii) Let { A_m}
be a sequence of positive constants. If ∥ X_m∥=O_P(A_m)
conditional on Y_m, namely, that for any ℓ_m→∞,
Pr(∥ X_m∥>ℓ_mA_m| Y_m)p→0,
then ∥ X_m∥=O_P(A_m) unconditionally,
namely, that for any ℓ_m→∞, Pr(∥ X_m∥>ℓ_mA_m)→0. 

PROOF: This lemma is the Lemma 6.1 in <cit.>.




 
SIMULATION
 








































