The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
81,735 characters
Estimation and Inference for Policy Relevant Treatment Effects
\title{Estimation and Inference for Policy Relevant Treatment Effects\thanks{First arXiv date: May 29, 2018.
We benefited from very useful comments by
Serena Ng and Elie Tamer (Editors), an anonymous Associate Editor, anonymous referees, other numerous researchers, and
participants at
UC Davis,
2018 CEME Conference: ``Inference on Nonstandard Problems,''
the 5th Annual Seattle-Vancouver Econometrics Conference, and
2018 California Econometrics Conference. All remaining errors are ours.}}
\author{Yuya Sasaki\thanks{Y. Sasaki: Department of Economics, Vanderbilt University, VU Station B \#351819, 2301 Vanderbilt Place, Nashville, TN 37235-1819. Email: [email removed]}\\ Department of Economics\\ Vanderbilt University
\and
Takuya Ura\thanks{T. Ura: Department of Economics, University of California, Davis, One Shields Avenue, Davis, CA 95616. Email: [email removed]}\\ Department of Economics\\ University of California, Davis}
\date{}
\maketitle
\begin{abstract}
The policy relevant treatment effect (PRTE) measures the average effect of switching from a status-quo policy to a counterfactual policy.
Estimation of the PRTE involves estimation of multiple preliminary parameters, including propensity scores, conditional expectation functions of the outcome and covariates given the propensity score, and marginal treatment effects.
These preliminary estimators can affect the asymptotic distribution of the PRTE estimator in complicated and intractable manners.
In this light, we propose an orthogonal score for double debiased estimation of the PRTE, whereby the asymptotic distribution of the PRTE estimator is obtained without any influence of preliminary parameter estimators as far as they satisfy mild requirements of convergence rates.
To our knowledge, this paper is the first to develop limit distribution theories for inference about the PRTE.
\begin{description}
\item {\bf Keywords:} double debiased estimation, orthogonal score, policy relevant treatment effects
\item {\bf JEL Codes:} C14, C21
\end{description}
\end{abstract}
\section{Introduction}
The policy relevant treatment effect \citep[PRTE,][]{heckman/vytlacil:2001,heckman/vytlacil:2005,heckman/vytlacil:2007} measures the average effect of switching from a status-quo policy to a counterfactual policy.\footnote{Also see \citet{stock:1989} and \citet{ichimura/taber:2000} for parameters related to the PRTE.}
Among a number of measures for treatment effects, the PRTE has the advantage of directly evaluating alternative policy scenarios under consideration.
A rich set of identification results have been established in the literature for this treatment effect parameter \citep{heckman/vytlacil:2001,heckman/vytlacil:2005,heckman/vytlacil:2007} based on the marginal treatment effects \citep{bjorklund/moffitt:1987}.
Estimation of the PRTE involves nonparametric estimation for multiple preliminary parameters.
First, we nonparametrically estimate the propensity score, the conditional probability of being treated given an instrumental variable.
Next, when the outcome model involves additive controls as in \citet[][Section 4.1]{carneiro/lee:2009} and \citet[][Section 8]{carneiro/heckman/vytlacil:2010}, we estimate the conditional expectation functions of the outcome and covariates given the propensity score through the partially linear regression estimation procedure of \citet{robinson1988root}.
Third, with the residual from the partially linear regression, we estimate the marginal treatment effect via a nonparametric derivative estimation.
These first-, second-, and third-stage estimates can affect the asymptotic distribution of the PRTE estimator in complicated and intractable manners, especially when their convergence rates are slower than the parametric rate of $n^{-1/2}$.
Estimation of the PRTE based on an orthogonal score can mitigate and asymptotically vanish the estimation errors of these first, second, and third preliminary parameters.
In this light, we propose an orthogonal score for double debiased estimation of the PRTE, whereby these first-, second-, and third-stage estimators will not affect the first-order asymptotic distribution of the PRTE estimator.
Consequently, we characterize the asymptotic distribution for the PRTE estimator based on the proposed orthogonal score.
Even conventionally in the absence of orthogonal scores, there has been no semi- or non-parametric inference method for the PRTE available in the existing literature, perhaps because of the complicated dependence of this parameter on multiple preliminary functions.
Therefore, to our knowledge, our work is the first to develop limit distribution theories for inference about the PRTE.
Orthogonal scores for double debiased and doubly robust estimation appear in a few contexts in the literature of econometrics and statistics.
One instance is about partial linear models,\footnote{Estimation of partial linear models is studied by \citet{robinson1988root} under exogeneity, and is studied by \citet{okui2012doubly} under endogeneity. For a robust minimum-distance approach under endogeneity, see \citet{ai2007estimation}. See \citet{belloni2013honest}, \citet{belloni2014high}, and \citet{farrell2015robust} for machine learning approaches.} and another is about weighted averages\footnote{Weighted estimation is studied by \citet{rosenbaum1987model}, \citet{robins1992estimating}, \citet{imbens1992efficient}, \citet{robins1994estimation}, \citet{robins1995semiparametric}, \citet{lawless1999}, \citet{wooldridge1999,wooldridge2001asymptotic}, \citet{lee/okui/whang:2017} and \citet{sloczynski_wooldridge_2018} under parametric weight function, is studied by \citet{Hahn:1998}, \citet{hirano2003efficient}, \citet{suzukawa2004unbiased}, \citet{NeweyRuud2005}, \citet{Firpo:2007}, \citet{lewbel2007simple}, \citet{magnac2007identification}, and \citet{chen2008semiparametric}, \citet{graham2011efficiency}, \citet{graham2012inverse}, \citet{graham2016efficient}, \citet{sant2018doubly}, and \citet{rothe2019properties} under nonparametric weight function, and is studied by \citet{BCFH:2017} and \citet{Wager/Athey:2018} with machine learning. Particularly, \citet{rothe2019properties} analyzes the orthogonal scores for the policy effect parameter defined in \citet{stock:1989}.
} such as inverse propensity score weighting.
Estimation of the PRTE consists of both of these two types of estimation procedures: (i) weighted estimation with the propensity score, estimation of a partial linear model to account for additive covariates as in \citet[][Section 4.1]{carneiro/lee:2009} and \citet[][Section 8]{carneiro/heckman/vytlacil:2010}, and (ii) weighted estimation with a counterfactual policy.
As such, it is a natural idea to develop an orthogonal score for estimation and inference for the PRTE by taking advantage of and combining these techniques.
With this said, the existing literature has not proposed an orthogonal score for the PRTE, and hence we aim to contribute to this literature by proposing one in this paper.
There is a large literature on orthogonal scores, and it characterizes desired properties of orthogonal scores such as double robustness against misspecification, local robustness in semiparametric estimation, and even semiparametric efficiency for certain orthogonal scores -- \citet{hitomi2008puzzling} discuss these features of orthogonal scores in general from geometric perspectives.
Among these characteristics, our purpose of developing an orthogonal score is, as mentioned earlier, to remove the influence of the first, second, and third-stage estimators on the asymptotic distribution of the PRTE estimator, which concerns the local robustness.
There are a number of general procedures for estimation and inference based on orthogonal scores which we rely on in the present paper.
The seminal paper by \citet{newey1994asymptotic} proposes orthogonal scores in semiparametric methods and provides forms of adjustment terms to to obtain orthogonal scores from moment functions.
The way in which we derive our orthogonal score is based on his prescription \citep[cf.][Proposition 4]{newey1994asymptotic}.
\citet{belloni2014uniform} and \citet{belloni2018uniformly} propose a generic Z-estimation framework based on orthogonal scores.
\citet{chernozhukov2016locally} propose a general procedure for construction of orthogonal scores from moment restriction models.
\citet{chernozhukov2018double} combine the use of orthogonal scores and cross-fitting as a generic semiparametric strategy.
Our proposed method of estimation and inference takes advantages of this existing body of knowledge.
To the best of our knowledge, no limit distribution has been established for the PRTE in the existing literature.
Hence, this work is the first to develop asymptotic theories for inference about the PRTE.
The closest is \citet{carneiro/lee:2009} who develop point-wise limit distributions for the MTE, although they do not seem to imply or lead to a limit distribution for the PRTE.
Our orthogonal score does not only pave the way for a limit distribution for the PRTE for the first time in the literature, but also allows for a flexibility in types of preliminary estimators, e.g., kernel, sieve, and lasso.
In order to make our method and theory more accessible in practice, we also provide specific estimation procedures and lower-level sufficient conditions.\footnote{We consider a generalized linear model for the propensity score function allowing for possibly high-dimensional covariates and propose primitive conditions based on nonparametric estimation of other preliminary functional parameters. This framework involves both shrinkage and kernel estimation, and is related to a branch of the recent literature including \citet{kennedy:2017}, \citet{fan:2019}, \citet{su/ura/zhang:2019}, \citet{zimmert/lechner:2019} and \citet{colangelo/lee:2020}.}
Our proposed method is applicable to empirical studies that identify and estimate marginal treatment effects and/or policy relevant treatment effects.
Examples include, but are not limited to,
\citet{auld:2005},
\citet{basu/heckman/navarro/urzua:2007},
\citet{doyle:2008},
\citet{moffitt:2008},
\citet{chuang/lai:2010},
\citet{carneiro/heckman/vytlacil:2011},
\citet{galasso/schankerman/serrano:2013},
\citet{basu/jena/goldman/philipson/dubois:2014},
\citet{belskaya/peter/posso:2014},
\citet{johar/maruyama:2014},
\citet{lindquist/santavirta:2014},
\citet{moffitt:2014},
\citet{dobbie/song:2015},
\citet{jensen/nielsen:2016},
\citet{kasahara/liang/rodrigue:2016},
\citet{carneiro/lokshin/umapathi:2017},
\citet{cornelissen/dustmann/raute/schonberg:2018},
\citet{felfe/lalive:2018},
and
\citet{kamhofer/schmitz/westphal:2018}.
The rest of this paper is organized as follows.
Section \ref{sec:model} introduces a model based on \citet{carneiro/lee:2009}.
Section \ref{sec:overview} presents an informal overview of our proposed method.
Section \ref{sec:orthogonal_score} presents an orthogonal score for double debiased estimation of the PRTE.
Section \ref{sec:large_sample} discusses large sample properties of the double debiased estimator for inference.
Section \ref{sec:simulations} presents Monte Carlo simulation studies.
Section \ref{sec:empirical} presents an empirical illustration.
The paper is summarized in Section \ref{sec:summary}.
The appendix contains proofs and additional details that are important but relegated there due to their lengths.
\section{Model}\label{sec:model}
Following the model and notations in \citet[][Section 2]{carneiro/lee:2009}, consider the structure
\begin{align*}
Y=SY_1+\left(1-S\right)Y_0
\end{align*}
of outcome production, where $Y_1$ and $Y_0$ denote the potential outcomes under treatment ($S=1$) and no treatment ($S=0$), respectively.
With covariates $X$, the potential outcomes are modeled by $Y_1=\mu_1\left(X,U_1\right)$ and $Y_0=\mu_0\left(X,U_0\right)$ with unobserved variables $\left(U_0,U_1\right)$.
The binary treatment assignment status $S$ is in turn determined by the threshold crossing model
$S=1\{\mu_S\left(Z\right)-U_S>0\}$, where $Z$ is a vector of exogenous variables and $U_S$ is an error term.
$Z$ can contain a sub-vector of $X$ as included exogenous variables, and the remaining excluded exogenous variables serve as instruments.
In this model, the dependence between $U_S$ and $\left(U_0,U_1\right)$ is the source of endogeneity in the treatment selection.
Researchers observe the random vector $W = \left(Y,S,X',Z'\right)'$, but do not observe $Y_1$, $Y_0$, $U_1$, $U_0$ or $U_S$.
Following \citet[][Section 3]{carneiro/lee:2009}, we further introduce the following notations for convenience.
Let $V=F_{U_S}\left(U_S\right)$ denote unobserved innate propensity to select into treatment, and let
$P=F_{U_S}\left(\mu_S\left(Z\right)\right)$ denote the treatment selection probability or the propensity score.
Under these setup and notations, we adopt the following conventional assumption.
\begin{assumption}[\citeauthor{carneiro/lee:2009}, \citeyear{carneiro/lee:2009}, Assumptions 1 and 2]\label{assn1}
(1)
$\mu_S\left(Z\right)$ is non-degenerate conditional on $X$.
(2)
$\left(U_1,U_S\right)$ and $\left(U_0,U_S\right)$ are independent of $\left(Z,X\right)$.\footnote{We note that Assumption 1 of \citet{carneiro/lee:2009} only assumes that $\left(U_1,U_S\right)$ and $\left(U_0,U_S\right)$ are independent of $Z$ conditionally on $X$ for the purpose of identification, while their Assumption 2 assumes that $\left(U_1,U_S\right)$ and $\left(U_0,U_S\right)$ are independent of $\left(Z,X\right)$ similarly to our assumption for the purpose of semiparametric estimation.}
(3)
The distribution of $U_S$ and the conditional distribution of $\mu_S\left(Z\right)$ given $X$ are absolutely continuous with respect to the Lebesgue measure.
(4)
$Y_1$ and $Y_0$ have finite first moments.
(5)
$0<Pr\left(S=1\mid Z\right)<1$.
(6)
$p \mapsto E[U_1\mid P=p,S=1]$, $p \mapsto E[U_0\mid P=p,S=0]$, $\left(u_1,p\right) \mapsto f_{U_1\mid P,S=1}\left(u_1\mid p\right)$ and $\left(u_0,p\right) \mapsto f_{U_0\mid P,S=0}\left(u_0\mid p\right)$ are continuously differentiable with respect to $p$.
\end{assumption}
We can write the treatment selection model in this setting by
\begin{align*}
S=1\{P>V\}.
\end{align*}
This representation provides the interpretation that those individuals with $V=P$ are at the margin of indifference between the two treatment statuses, and motivates the marginal treatment effect \citep[][]{bjorklund/moffitt:1987} defined by
\begin{align*}
MTE\left(x,p\right)=E[Y_1-Y_0\mid X=x,V=p].
\end{align*}
This treatment parameter serves as a building block for many treatment parameters \citep{heckman/vytlacil:1999,heckman/vytlacil:2001,heckman/vytlacil:2005}.
We consider a counterfactual propensity score $P^\ast=P^\ast\left(P,Z\right)$ with a known function $P^\ast\left(\cdot,\cdot\right)$ satisfying Assumption \ref{assn:pstar} to be stated ahead.\footnote{\citet[][Section 3.2]{carneiro/heckman/vytlacil:2010} consider alternative specifications for counterfactual policies $P^\ast$. One scenario takes the form of $P^\ast=P^\ast\left(P,Z\right)$, in which a policy maker directly affects the probability of being in the treatment group. Another scenario takes the form of $P^\ast=P\left(Z^\ast\left(Z\right)\right)$, in which the counterfactual policy affects the distribution of $Z$. We consider the former scenario in the main text, whereas we consider the latter scenario in Appendix \ref{sec:alternative_P_star}.}
Under a counterfactual propensity score $P^\ast$, the counterfactual treatment and outcome are
\begin{align*}
S^\ast&=1\{P^\ast>V\}
\qquad\text{and}\\
Y^\ast&=S^\ast Y_1+\left(1-S^\ast\right)Y_0.
\end{align*}
The parameter of our interest is the policy relevant treatment effect \citep[PRTE;][]{heckman/vytlacil:1999,heckman/vytlacil:2001,heckman/vytlacil:2005}:
\begin{equation}
PRTE=\frac{E[Y^\ast]-E[Y]}{E[S^\ast]-E[S]}.
\end{equation}
This treatment parameter is of policy interest because it allows to compare alternative virtual policies $P^\ast$ under consideration.
\cite{heckman/vytlacil:1999,heckman/vytlacil:2001,heckman/vytlacil:2005} has demonstrated that $PRTE$ can be represented as the weighted average of $MTE\left(x,p\right)$:
$$
PRTE=\int\int_0^1MTE\left(x,p\right)\frac{F_{P\mid X}\left(p\mid x\right)-F_{P^\ast\mid X}\left(p\mid x\right)}{E[P^\ast]-E[P]}dpf_X\left(x\right)dx.
$$
Note that $MTE\left(x,p\right)$ measures the average treatment effects for the subpopulation of those individuals with $V=p$ and $F_{P|X}\left(p|x\right)-F_{P^\ast|X}\left(p|x\right)$ measures the probability of compliance with the counterfactual policy for this subpopulation.
Therefore, setting aside the denominator $E[P^\ast]-E[P]$, this integral measures the average treatment effect of the counterfactual policy for the population.
A division of quantity by $E[P^\ast]-E[P]$ in turn yields the average treatment effect among those who complied with the policy.
In the rest of this paper, we propose a method of double debiased estimation and inference for the PRTE.
Since the identification of the PRTE hinges on that of the marginal treatment effect, the theories of our estimation and inference methods rely on the prior result of the identification of the marginal treatment effect, which is formally stated in the following theorem.
\begin{theorem}[\citeauthor{heckman/vytlacil:1999}, \citeyear{heckman/vytlacil:1999,heckman/vytlacil:2001,heckman/vytlacil:2005}; and \citeauthor{carneiro/lee:2009}, \citeyear{carneiro/lee:2009}]\label{theorem:MTE_LIV}
If Assumption \ref{assn1} is satisfied, then
\begin{align*}
MTE\left(x,p\right)=\frac{\partial E[Y\mid X=x,P=p]}{\partial p}
\end{align*}
holds provided that $p \mapsto E[Y\mid X,P=p,S=1]$ and $p \mapsto E[Y\mid X,P=p,S=0]$ are continuously differentiable with respect to $p$ almost surely (with respect to $X$).
\end{theorem}
In actual empirical research, researchers often use observed controls $X$.
Without imposing additional structural restrictions, they would suffer from the curse of dimensionality of $X$ in estimation of the marginal treatment effect.
The existing literature proposes alternative suggestions to address this issue.
Following \citet[][Section 4.1]{carneiro/lee:2009} and \citet[][Section 8]{carneiro/heckman/vytlacil:2010}, we employ the following structural restriction on $\mu_1$ and $\mu_0$.
\begin{assumption}[\citeauthor{carneiro/lee:2009}, \citeyear{carneiro/lee:2009}, Section 4.1]\label{assn2}
$Y_1=\mu_1\left(X\right)'\beta_1+U_1$ and $Y_0=\mu_0\left(X\right)'\beta_0+U_0$
for $d_X$-dimensional unknown parameters $\beta_1$ and $\beta_0$ and known functions $\mu_1$ and $\mu_0$.
\end{assumption}
Under Assumptions \ref{assn1} and \ref{assn2}, we can write the conditional mean of $Y$ given $\left(X,P\right)$ as
$$
E[Y\mid X,P]=P\mu_1\left(X\right)'\beta_1+\left(1-P\right)\mu_0\left(X\right)' \beta_0+E[U\mid P],
$$
where $U=SU_1+\left(1-S\right)U_0$.
Thus the marginal treatment effect takes the partial linear form of
\begin{align}\label{eq:partial_linear_MTE}
MTE\left(x,p\right) = \mu_1\left(x\right)' \beta_1 - \mu_0\left(x\right)' \beta_0 + \Delta_{U\mid P}\left(p\right),
\end{align}
where $\Delta_{U\mid P}\left(p\right) = \frac{d}{dp} E[U|P=p]$.
In the rest of the paper, we shall use this form (\ref{eq:partial_linear_MTE}) of the MTE expression for analysis of the PRTE:
\begin{align}
PRTE
=&
\frac{E[\left(P^\ast\left(\boldsymbol{g}_{S\mid Z}\left(Z\right),Z\right)-\boldsymbol{g}_{S\mid Z}\left(Z\right)\right)\mu_1\left(X\right)]'}{E[P^\ast\left(\boldsymbol{g}_{S\mid Z}\left(Z\right),Z\right)-\boldsymbol{g}_{S\mid Z}\left(Z\right)]}\beta_1-\frac{E[\left(P^\ast\left(\boldsymbol{g}_{S\mid Z}\left(Z\right),Z\right)-\boldsymbol{g}_{S\mid Z}\left(Z\right)\right)\mu_0\left(X\right)]}{E[P^\ast\left(\boldsymbol{g}_{S\mid Z}\left(Z\right),Z\right)-\boldsymbol{g}_{S\mid Z}\left(Z\right)]}'\beta_0
\notag\\&+\frac{E\left[\Delta_{U\mid P}\left(\boldsymbol{g}_{S\mid Z}\left(Z\right)\right)\frac{F_{P}\left(\boldsymbol{g}_{S\mid Z}\left(Z\right)\right)-F_{P^\ast}\left(\boldsymbol{g}_{S\mid Z}\left(Z\right)\right)}{f_{P}\left(\boldsymbol{g}_{S\mid Z}\left(Z\right)\right)}\right]}{E[P^\ast\left(\boldsymbol{g}_{S\mid Z}\left(Z\right),Z\right)-\boldsymbol{g}_{S\mid Z}\left(Z\right)]},\label{eq:PRTE_expression}
\end{align}
where $\boldsymbol{g}_{S\mid Z}\left(z\right) = E[S|Z=z]$.
Appendix \ref{sec:eq:partial_linear_PRTE_rewritten} provides a derivation of \eqref{eq:PRTE_expression}.
We close this section by discussing feasible and infeasible extensions and variants of our parameter, $PRTE$, of interest.
First, the PRTE may be defined as $E[\omega\left(X\right)PRTE\left(X\right)]$ with a known weight function $\omega$ for heterogeneous policy designs.
Our method and theory presented ahead extends to this variant of the PRTE by replacing $Y$ by $\omega\left(X\right)Y$.
Second, our framework also extends to analysis of $PRTE\left(x\right)$ when $X$ is discrete by focusing on the subpopulation with $X=x$.
When $X$ is continuous, however, it is generally infeasible to estimate $PRTE\left(x\right)$ at the root-$n$ convergence rate.
\section{An Overview of the Method}\label{sec:overview}
While we present a full-fledged framework and formal asymptotic theories in Sections \ref{sec:orthogonal_score} and \ref{sec:large_sample}, respectively, we first provide an informal overview in this section focusing on a simplified model without covariates.
Define $\theta = \left(\theta_N,\theta_D\right)$ by
\begingroup
\allowdisplaybreaks
\begin{align*}
\theta_N =&
E\left[
\Delta_{U\mid P}\left(\boldsymbol{g}_{S\mid Z}\left(Z\right)\right)\frac{{F}_{P}\left(\boldsymbol{g}_{S\mid Z}\left(Z\right)\right)-{F}_{P^\ast}\left(\boldsymbol{g}_{S\mid Z}\left(Z\right)\right)}{f_{P}\left(\boldsymbol{g}_{S\mid Z}\left(Z\right)\right)}
\right]
\\
\theta_D =&
E\left[
P^\ast\left(\boldsymbol{g}_{S\mid Z}\left(Z\right),Z\right)-\boldsymbol{g}_{S\mid Z}\left(Z\right)
\right],
\end{align*}
\endgroup
where $\Delta_{U\mid P}\left(p\right) = \frac{d}{dp} E[U|P=p]$ and $\boldsymbol{g}_{S\mid Z}\left(z\right) = E[S|Z=z]$ from Section \ref{sec:model}.
In the current setting without $X$, the PRTE consists of the last term in \eqref{eq:PRTE_expression} and can be written as
\begin{align}\label{eq:prte_n_over_d}
PRTE = \theta_N / \theta_D.
\end{align}
We collect the possibly infinite-dimensional nuisance parameters as
\begin{align*}
\gamma=
\left(
\frac{f_{P^\ast}\left({\boldsymbol{g}}_{S\mid Z}\right)}{f_{P}\left({\boldsymbol{g}}_{S\mid Z}\right)},
\boldsymbol{g}_{S\mid Z},
\boldsymbol{g}_{U\mid P},
\Delta_{U\mid P}
\right)
\qquad
\text{where $\boldsymbol{g}_{U\mid Z}\left(z\right) = E[U|Z=z]$.}
\end{align*}
\subsection{A Step-by-Step Procedure}\label{sec:overview_step_by_step}
Let $L>1$ be a predetermined natural number of folds that is fixed as the sample size increases.
Randomly split the sample into sub-samples $I_1,...,I_L$ of (approximately) equal size, i.e., $|I_\ell| = \lfloor n/L \rfloor$ or $\lfloor n/L \rfloor+1$.
For each $\ell \in L$, estimate $\gamma$ by using the sub-sample $I_\ell^c$, and denote the estimator by $\hat\gamma_\ell$.\footnote{See Appendix \ref{sec:step_by_step_procedure_application} for concrete estimators along with tuning parameter choice rules used for them in our empirical application.}
We then define our double debiased estimator $\hat\theta=\left(\hat\theta_N,\hat\theta_D\right)$ of $\theta$ as the solution to
\begin{align*}
\frac{1}{L} \sum_{\ell=1}^L \frac{1}{|I_\ell|} \sum_{i \in I_\ell} m_N\left(Y_i,Z_i;\hat\theta_N,\hat\gamma_\ell\right) &= 0
\qquad\text{and}\\
\frac{1}{L} \sum_{\ell=1}^L \frac{1}{|I_\ell|} \sum_{i \in I_\ell} m_D\left(Y_i;\hat\theta_D,\hat\gamma_\ell\right) &= 0,
\end{align*}
where $\left(m_N,m_D\right)$ is the orthogonal score defined by
\begingroup
\allowdisplaybreaks
\begin{align}
m_N\left(Y,Z;\tilde\theta_N,\tilde\gamma\right)
=&
\tilde{\boldsymbol{g}}_{U\mid P}\left(P^\ast\left(\tilde{\boldsymbol{g}}_{S\mid Z}\left(Z\right),Z\right)\right)-Y-\tilde\theta_N
\label{eq:theta_3_a}\\
&+
\frac{\tilde{f}_{P^\ast}\left(\tilde{\boldsymbol{g}}_{S\mid Z}\left(Z\right)\right)}{\tilde{f}_{P}\left(\tilde{\boldsymbol{g}}_{S\mid Z}\left(Z\right)\right)}\left(Y-\tilde{\boldsymbol{g}}_{U\mid P}\left(\tilde{\boldsymbol{g}}_{S\mid Z}\left(Z\right)\right)\right)
\label{eq:theta_3_b}\\&+
\left(
\tilde{\Delta}_{U\mid P}\left(P^\ast\left(\tilde{\boldsymbol{g}}_{S\mid Z}\left(Z\right),Z\right)\right)\partial P^\ast\left(\tilde{\boldsymbol{g}}_{S\mid Z}\left(Z\right),Z\right)
-
\frac{\tilde{f}_{P^\ast}\left(\tilde{\boldsymbol{g}}_{S\mid Z}\left(Z\right)\right)}{\tilde{f}_{P}\left(\tilde{\boldsymbol{g}}_{S\mid Z}\left(Z\right)\right)}\tilde{\Delta}_{U\mid P}\left(\tilde{\boldsymbol{g}}_{S\mid Z}\left(Z\right)\right)
\right)\nonumber
\\
&\times
\left(S-\tilde{\boldsymbol{g}}_{S\mid Z}\left(Z\right)\right),
\label{eq:theta_3_c}
\\
m_D\left(Z;\tilde\theta_D,\tilde\gamma\right)
=&
P^\ast\left(\tilde{\boldsymbol{g}}_{S\mid Z}\left(Z\right),Z\right)-\tilde{\boldsymbol{g}}_{S\mid Z}\left(Z\right)-\tilde\theta_D
\label{eq:theta_2_a}\\
&+
\left(\partial P^\ast\left(\tilde{\boldsymbol{g}}_{S\mid Z}\left(Z\right),Z\right)-1\right)
\left(S-\tilde{\boldsymbol{g}}_{S\mid Z}\left(Z\right)\right).
\label{eq:theta_2_b}
\end{align}
\endgroup
Plug $\hat\theta$ into \eqref{eq:prte_n_over_d} to in turn obtain an estimator of the PRTE:
$$
\widehat{PRTE} = \hat\theta_N / \hat\theta_D.
$$
This double debiased estimator for the PRTE is root-$n$ asymptotically normal as
$$
\sqrt{n}\left(\widehat{PRTE}-PRTE\right) \rightarrow_d N\left(0,D' \Omega D\right),
$$
where
$
D = \left(1/\theta_D, \ -\theta_N/\theta_D^2\right)'
$
and
$
\Omega = \text{Var}\left(\left( m_N\left(Y,Z;\theta_N,\gamma\right), \ m_D\left(Y;\theta_d,\gamma\right)\right)'\right).
$
\subsection{Intuitions and Discussions}\label{sec:intuitions_discussions}
Having outlined the step-by-step procedure of our proposed method, we next present intuitions behind this proposed method of double debiased estimation and inference for the PRTE as well as some heuristic explanations of why it works.
{\bf On the orthogonal score:}
If we knew the nuisance parameters $(\boldsymbol{g}_{S|Z}, \boldsymbol{g}_{U|P})$, then we could estimate $\theta$ by using the moment conditions $E[M_N\left(Y,Z;\theta,\gamma\right)]=E[M_D\left(Y,Z;\theta,\gamma\right)]=0$, where
\begin{align}
M_N\left(Y,Z;\tilde\theta_N,\tilde\gamma\right)=&
\tilde{\boldsymbol{g}}_{U\mid P}\left(P^\ast\left(\tilde{\boldsymbol{g}}_{S\mid Z}\left(Z\right),Z\right)\right)-Y-\tilde{\theta}_N
\label{eq:momentMN}
\qquad\text{and}\\
M_D\left(Y;\tilde\theta_D,\tilde\gamma\right)=&
P^\ast\left(\tilde{\boldsymbol{g}}_{S\mid Z}\left(Z\right),Z\right)-\tilde{\boldsymbol{g}}_{S\mid Z}\left(Z\right)-\tilde{\theta}_D.
\label{eq:momentMD}
\end{align}
Since we estimate $(\boldsymbol{g}_{S|Z}, \boldsymbol{g}_{U|P})$, however, influence function adjustments should be added as in \eqref{eq:theta_3_b}, \eqref{eq:theta_3_c} and \eqref{eq:theta_2_b}.
Lines \eqref{eq:theta_3_a} and \eqref{eq:theta_2_a} precisely correspond to the na\"ive moment functions \eqref{eq:momentMN} and \eqref{eq:momentMD}, respectively.
An influence function adjustment for the estimation of ${\boldsymbol{g}}_{U\mid P}$ appears in \eqref{eq:theta_3_b}.
An estimation of ${\boldsymbol{g}}_{U\mid P}\left(P\right)$ entails an influence function adjustment by $Y-{\boldsymbol{g}}_{U\mid P}\left({\boldsymbol{g}}_{S\mid Z}\left(Z\right)\right)$.
Since the function ${\boldsymbol{g}}_{U\mid P}$ is evaluated at $P^\ast\left({\boldsymbol{g}}_{S\mid Z}\left(Z\right),Z\right)$ in line \eqref{eq:theta_3_a}, the adjustment is scaled by the coefficient ${{f}_{P^\ast}\left({\boldsymbol{g}}_{S\mid Z}\left(Z\right)\right)}/{{f}_{P}\left({\boldsymbol{g}}_{S\mid Z}\left(Z\right)\right)}$ as in \eqref{eq:theta_3_b}.
Likewise, an estimation of ${\boldsymbol{g}}_{S\mid Z}\left(Z\right)$ entails an influence function adjustment by $S-{\boldsymbol{g}}_{S\mid Z}\left(Z\right)$, and it appears in lines \eqref{eq:theta_3_c} and \eqref{eq:theta_2_b}.
Since the function ${\boldsymbol{g}}_{S\mid Z}\left(Z\right)$ shows up inside other functions in lines \eqref{eq:theta_3_a} and \eqref{eq:theta_2_a}, the adjustments are scaled by their derivatives in lines \eqref{eq:theta_3_c} and \eqref{eq:theta_2_b}.
Thanks to these influence function adjustments, these moment functions $m_N$ and $m_D$ are robust against local perturbations in $\gamma$, i.e., $m_N$ and $m_D$ constitute an orthogonal score.
This property, together with the cross fitting discussed below, allows the effects of an estimation of $\gamma$ on an estimation of $\theta$ to be asymptotically negligible.
{\bf On the cross fitting:}
Using the same sample to estimate both $\gamma$ and $\theta$ would result in an over-fitting bias.
To circumvent such a bias, we use the sub-sample $I_\ell^c$ to estimate $\gamma$ by $\hat\gamma_\ell$ and use the complementary sub-sample $I_\ell$ to evaluate the orthogonal score, $|I_\ell|^{-1}\sum_{i\in I_\ell} m_N\left(Y_i,Z_i;\theta_N,\hat\gamma_\ell\right)$ and $|I_\ell|^{-1}\sum_{i\in I_\ell} m_D\left(Y_i;\theta_D,\hat\gamma_\ell\right)$, for each $\ell = 1,...,L$.
These lead to the double debiased estimator $\hat\theta$ introduced in Section \ref{sec:overview_step_by_step}.
{\bf On the root-$n$ asymptotic normality:}
Because of the orthogonality property and the cross fitting that allow the effects of the estimation error of $\gamma$ on an estimation of $\theta$ to be asymptotically negligible, $\hat\theta$ enjoys the root-$n$ asymptotic normality
\begin{align*}
\sqrt{n}\left(\hat\theta-\theta\right) \rightarrow_d N\left(0,\Omega\right)
\end{align*}
as if $\gamma$ were known.
The delta method thus yields the root-$n$ asymptotic normality for the PRTE estimator:
$$
\sqrt{n}\left(\widehat{PRTE}-PRTE\right) \rightarrow_d N\left(0,D' \Omega D\right),
$$
again as if $\gamma$ were known.
{\bf Orthogonal Score and Double Robustness:}
The orthogonality property and double robustness are closely related concepts, but neither implies the other -- see \citet{chernozhukov2016locally} for example.
A natural question is whether our orthogonal score also satisfies the double robustness with respect to the nuisance parameters, $f_{P^\ast}\left({\boldsymbol{g}}_{S\mid Z}\right)/f_P\left({\boldsymbol{g}}_{S\mid Z}\right)$, ${\boldsymbol{g}}_{S\mid Z}$, and $\left({\boldsymbol{g}}_{U\mid P},\Delta_{U\mid P}\right)$.
It turns out that our orthogonal score does not possess the double robustness property against the propensity score function ${\boldsymbol{g}}_{S\mid Z}$ in general.
This is because ${\boldsymbol{g}}_{S\mid Z}$ appears in our score inside possibly nonlinear functions -- we show a concrete case in point in Appendix \ref{sec:orthogonal_score_double_robustness}.
On the other hand, given a fixed propensity score function ${\boldsymbol{g}}_{S\mid Z}$, our orthogonal score has double robustness between $f_{P^\ast}\left({\boldsymbol{g}}_{S\mid Z}\right)/f_P\left({\boldsymbol{g}}_{S\mid Z}\right)$ and $\left({\boldsymbol{g}}_{U\mid P},\Delta_{U\mid P}\right)$.
This is because the score takes forms of products of affine functions of these nuisance parameters -- see Appendix \ref{sec:orthogonal_score_double_robustness} for details.
{\bf More intuitions:}
Using this discussion on the relation to the double robustness, we now provide more intuitions behind why the orthogonal score works to our goal.
One may wonder why our PRTE estimator converges at the rate of root-$n$ while possibly all the preliminary estimators converge at slower rates.
As demonstrated in Appendix \ref{sec:orthogonal_score_double_robustness}, our orthogonal score after some rewritings takes the form of a product of two estimation errors as in the right-hand side of
\begin{eqnarray*}
E[m_N\left(Y,Z;\theta_N,\tilde\gamma\right)]
=
\int\left(
\frac{\tilde{f}_{P^\ast}\left(p\right)}{\tilde{f}_{P}\left(p\right)}-\frac{{f}_{P^\ast}\left(p\right)}{{f}_{P}\left(p\right)}\right)\left({\boldsymbol{g}}_{U\mid P}\left(p\right)-\tilde{\boldsymbol{g}}_{U\mid P}\left(p\right)\right)f_{P}\left(p\right)dp.
\end{eqnarray*}
See Appendix \ref{sec:orthogonal_score_double_robustness} for its derivation.
This shows that, if each of
$
\tilde{f}_{P^\ast}\left( \cdot \right)/\tilde{f}_{P}\left( \cdot \right)-{f}_{P^\ast}\left( \cdot \right)/{f}_{P}\left( \cdot \right)
$
and
$
{\boldsymbol{g}}_{U\mid P}\left( \cdot \right)-\tilde{\boldsymbol{g}}_{U\mid P}\left( \cdot \right)
$
converges to zero at a rate faster than $n^{1/4}$ (yet slower than $n^{1/2}$), then the product converges at a rate faster than $n^{1/4} \cdot n^{1/4} = n^{1/2}$.
This property of the orthogonality (taking the form of a product in this case) solves the puzzle that the PRTE estimator of our interest can converge at the rate of root-$n$ while convergence rates of the preliminary estimators are often slower.
\section{Orthogonal Score and Double Debiased Estimation}\label{sec:orthogonal_score}
From Theorem \ref{theorem:MTE_LIV}, we can express the PRTE as a function of estimable moments.
Consider the vector, $\theta=\left(\theta_1',\theta_2',\theta_3\right)'$, of estimable moments defined by
\begingroup
\allowdisplaybreaks
\begin{align*}
\theta_1
&=
E\left[\xi_1\left(X,Y,\boldsymbol{g}_{S\mid Z}\left(Z\right)\right)\right],
&&\left(2p\left(2p+1\right)\text{-dimensional}\right)
\\
\theta_2
&=
E\left[
\left(\mu_0\left(X\right)',\mu_1\left(X\right)',1\right)'\left(P^\ast\left(\boldsymbol{g}_{S\mid Z}\left(Z\right),Z\right)-\boldsymbol{g}_{S\mid Z}\left(Z\right)\right)
\right],
&&\left(\left(2p+1\right)\text{-dimensional}\right)
\\
\theta_3
&=
E\left[
\Delta_{U\mid P}\left(\boldsymbol{g}_{S\mid Z}\left(Z\right)\right)\frac{{F}_{P}\left(\boldsymbol{g}_{S\mid Z}\left(Z\right)\right)-{F}_{P^\ast}\left(\boldsymbol{g}_{S\mid Z}\left(Z\right)\right)}{f_{P}\left(\boldsymbol{g}_{S\mid Z}\left(Z\right)\right)}
\right],
&&\left(\text{1-dimensional}\right)
\end{align*}
\endgroup
where $\boldsymbol{g}_{S\mid Z}\left(z\right) = E[S|Z=z]$,
$\Delta_{U\mid P}\left(p\right) = \frac{d}{dp} E[U|P=p]$,
\begin{align*}
\xi_1\left(x,y,p\right)=\mathrm{vec}
\left(
\left(\begin{array}{c}
\left(1-p\right)\left(\mu_0\left(x\right)-\boldsymbol{g}_{\mu_0\left(X\right)\mid P}\left(p\right)\right)\\
p\left(\mu_1\left(x\right)-\boldsymbol{g}_{\mu_1\left(X\right)\mid P}\left(p\right)\right)
\end{array}\right)
\left(\begin{array}{c}
\left(1-p\right)\left(\mu_0\left(x\right)-\boldsymbol{g}_{\mu_0\left(X\right)\mid P}\left(p\right)\right)\\
p\left(\mu_1\left(x\right)-\boldsymbol{g}_{\mu_1\left(X\right)\mid P}\left(p\right)\right)\\
y-\boldsymbol{g}_{Y\mid P}\left(p\right)
\end{array}\right)'
\right),
\end{align*}
$\boldsymbol{g}_{\mu_0\left(X\right)\mid P}\left(p\right) = E[\mu_0\left(X\right) | P=p]$,
$\boldsymbol{g}_{\mu_1\left(X\right)\mid P}\left(p\right) = E[\mu_1\left(X\right) | P=p]$, and
$\boldsymbol{g}_{Y\mid P}\left(p\right) = E[Y | P=p]$.
To express $PRTE$ in (\ref{eq:PRTE_expression}) as a function of $\theta$, we write
$$
PRTE
=
\Lambda\left(\theta\right)
\equiv
\frac{\theta_{2,1}' \boldsymbol{d}_1\left(\theta_1\right) - \theta_{2,0}' \boldsymbol{d}_0\left(\theta_1\right)+ \theta_3}{\theta_{2,2}},
$$
where
$\boldsymbol{d}=\left(\boldsymbol{d}_0',\boldsymbol{d}_1'\right)'$ is the function defined by
$\boldsymbol{d}\left(\mathrm{vec}\left(\mathbf{B},\mathbf{A}\right)\right)={\mathbf{B}^{-1}} {\mathbf{A}}$ for a $2p \times 1$ vector $\mathbf{A}$ and a $2p \times 2p$ matrix $\mathbf{B}$, cf. Equation \eqref{eq:beta_moments} below for the function $\boldsymbol{d}$.
Note that the definition of $\xi_1$ comes from the moment condition for $\left(\beta_0',\beta_1'\right)'$.
Specifically, as in \citet{robinson1988root}, the parameter vector $\left(\beta_0',\beta_1'\right)'$ can be written as
\begin{align}
\left(\beta_0',\beta_1'\right)'
&=
\left(\boldsymbol{d}_0\left(E\left[\xi_1\left(X,Y,\boldsymbol{g}_{S\mid Z}\left(Z\right)\right)\right]\right)',\boldsymbol{d}_1\left(E\left[\xi_1\left(X,Y,\boldsymbol{g}_{S\mid Z}\left(Z\right)\right)\right]\right)'\right)'
\label{eq:beta_moments}
\\
&=
E\left[
\left(\begin{array}{c}
\left(1-P\right)\left(\mu_0\left(X\right)-\boldsymbol{g}_{\mu_0\left(X\right)\mid P}\left(P\right)\right)\\
P\left(\mu_1\left(X\right)-\boldsymbol{g}_{\mu_1\left(X\right)\mid P}\left(P\right)\right)
\end{array}\right)
\left(\begin{array}{c}
\left(1-P\right)\left(\mu_0\left(X\right)-\boldsymbol{g}_{\mu_0\left(X\right)\mid P}\left(P\right)\right)\\
P\left(\mu_1\left(X\right)-\boldsymbol{g}_{\mu_1\left(X\right)\mid P}\left(P\right)\right)
\end{array}\right)'\right]^{-1}
\nonumber
\\& \qquad\times
E\left[\left(\begin{array}{c}
\left(1-P\right)\left(\mu_0\left(X\right)-\boldsymbol{g}_{\mu_0\left(X\right)\mid P}\left(P\right)\right)\\
P\left(\mu_1\left(X\right)-\boldsymbol{g}_{\mu_1\left(X\right)\mid P}\left(P\right)\right)
\end{array}\right)\left(Y-\boldsymbol{g}_{Y\mid P}\left(P\right)\right)\right].
\nonumber
\end{align}
Recall the notation $W = \left(Y,S,X',Z'\right)'$ from Section \ref{sec:model}.
With these definitions and notations, we now propose an orthogonal score function of the form
\begin{align*}
m\left(W;\tilde{\theta},\tilde{\gamma}\right)=\left(\begin{array}{ccccc}m_1\left(W;\tilde{\theta},\tilde{\gamma}\right)'& m_2\left(W;\tilde{\theta},\tilde{\gamma}\right)'& m_3\left(W;\tilde{\theta},\tilde{\gamma}\right)'\end{array}\right)'
\end{align*}
with all the preliminary parameters collected in the concise notation
\begin{align*}
\gamma=
\left(
\frac{f_{P^\ast}\left({\boldsymbol{g}}_{S\mid Z}\right)}{f_{P}\left({\boldsymbol{g}}_{S\mid Z}\right)},
\boldsymbol{g}_{S\mid Z},
\xi_1,
\zeta,
\boldsymbol{g}_{U\mid P},
\Delta_{U\mid P}
\right),
\end{align*}
where
$\zeta\left(z\right)=E\left[\left.\left.\frac{\partial}{\partial p}
\xi_1\left(X,Y,p\right)\right|_{p=\boldsymbol{g}_{S\mid Z}\left(z\right)}
\right\vert Z=z\right]$.
The components of the orthogonal score function are defined by
\begingroup
\allowdisplaybreaks
\begin{align}
m_1\left(W;\tilde{\theta},\tilde{\gamma}\right)
=&
\tilde{\xi}_1\left(X,Y,\tilde{\boldsymbol{g}}_{S\mid Z}\left(Z\right)\right)
-\tilde{\theta}_1
\label{eq:m_1_a}\\
&+
\tilde{\zeta}\left(Z\right)\left(S-\tilde{\boldsymbol{g}}_{S\mid Z}\left(Z\right)\right)
\label{eq:m_1_b}
\\
m_2\left(W;\tilde{\theta},\tilde{\gamma}\right)
=&
\left(\mu_0\left(X\right)',\mu_1\left(X\right)',1\right)'\left(P^\ast\left(\tilde{\boldsymbol{g}}_{S\mid Z}\left(Z\right),Z\right)-\tilde{\boldsymbol{g}}_{S\mid Z}\left(Z\right)\right)-\tilde{\theta}_2
\label{eq:m_2_a}\\
&+
\left(\mu_0\left(X\right)',\mu_1\left(X\right)',1\right)'\left(\partial P^\ast\left(\tilde{\boldsymbol{g}}_{S\mid Z}\left(Z\right),Z\right)-1\right)
\left(S-\tilde{\boldsymbol{g}}_{S\mid Z}\left(Z\right)\right)
\label{eq:m_2_b}
\\
m_3\left(W;\tilde{\theta},\tilde{\gamma}\right)
=&
\tilde{\boldsymbol{g}}_{U\mid P}\left(P^\ast\left(\tilde{\boldsymbol{g}}_{S\mid Z}\left(Z\right),Z\right)\right)-\mathcal{U}\left(W,\tilde{\theta}\right)-\tilde{\theta}_3
\label{eq:m_3_a}\\
&+
\frac{\tilde{f}_{P^\ast}\left(\tilde{\boldsymbol{g}}_{S\mid Z}\left(Z\right)\right)}{\tilde{f}_{P}\left(\tilde{\boldsymbol{g}}_{S\mid Z}\left(Z\right)\right)}\left(\mathcal{U}\left(W,\tilde{\theta}\right)-\tilde{\boldsymbol{g}}_{U\mid P}\left(\tilde{\boldsymbol{g}}_{S\mid Z}\left(Z\right)\right)\right)
\label{eq:m_3_b}\\&+
\left(
\tilde{\Delta}_{U\mid P}\left(P^\ast\left(\tilde{\boldsymbol{g}}_{S\mid Z}\left(Z\right),Z\right)\right)\partial P^\ast\left(\tilde{\boldsymbol{g}}_{S\mid Z}\left(Z\right),Z\right)
-
\frac{\tilde{f}_{P^\ast}\left(\tilde{\boldsymbol{g}}_{S\mid Z}\left(Z\right)\right)}{\tilde{f}_{P}\left(\tilde{\boldsymbol{g}}_{S\mid Z}\left(Z\right)\right)}\tilde{\Delta}_{U\mid P}\left(\tilde{\boldsymbol{g}}_{S\mid Z}\left(Z\right)\right)
\right)\nonumber
\\
&\times
\left(S-\tilde{\boldsymbol{g}}_{S\mid Z}\left(Z\right)\right),
\label{eq:m_3_c}
\end{align}
\endgroup
where $\mathcal{U}\left(W,\theta\right)=Y-\left(1-S\right)\mu_0\left(X\right)'\beta_0-S\mu_1\left(X\right)'\beta_1$ and $\partial P^\ast\left(p,z\right)=\frac{\partial}{\partial p}P^\ast\left(p,z\right)$.
See Appendix \ref{sec:decomposition_of_moment_functions} for an alternative representation of $m\left(W;\theta,\gamma\right)$
This score function $m$ is \textit{orthogonal} in the sense that
\begin{equation}\label{eq:orthogonality}
\frac{\partial}{\partial r}\int m\left(w;\theta,\gamma+r\left(\breve{\gamma}-\gamma\right)\right)F_W\left(dw\right)|_{r=0}=0
\qquad
\text{for all } \breve{\gamma} \in \Gamma.
\end{equation}
See Appendix \ref{sec:eq:orthogonality} for a derivation of this property as well as the definition of $\Gamma$.
To arrive at the orthogonal score, we take advantage of two convenient features of $\theta$ and $\gamma$.
First, all the nuisance parameters $\gamma$ take forms of mean-square projections or densities.
Second, the nuisance parameters $\gamma$ enter the equations for $\theta$ through integrals and derivatives, and therefore the Gateaux derivative operator (and thus the orthogonalization operator) can directly act on $\gamma$.
These two features allow us to obtain the adjustment terms using the result of \citet[][Proposition 4]{newey1994asymptotic}.
The orthogonality \eqref{eq:orthogonality} allows for the score to be insensitive to local perturbations $\breve{\gamma}$ of $\gamma$.
In particular, the first-order effect of the estimation error $\hat\gamma-\gamma$ is zero, and the remaining effects are of a smaller order:
$$
\int m\left(w;\theta,\hat{\gamma}\right)F_W\left(dw\right)=o_p\left(n^{-1/2}\right).
$$
In other words, the effects of the estimation error $\hat{\gamma}_\ell - \gamma$ of preliminary parameters are asymptotically negligible relative to the convergence rate of the empirical mean, which is of order $n^{-1/2}$.
Although we postpone rigorous discussions until Section \ref{sec:large_sample} and proofs of the main theorem therein, we here discuss intuitions behind the orthogonality property.
If we knew the true functions ${\xi}_1$, ${\boldsymbol{g}}_{U\mid P}$, and ${\boldsymbol{g}}_{S\mid Z}$, then we would use the moment conditions
$$
E\left[
\left(
\begin{array}{c}
{\xi}_1\left(X,Y,{\boldsymbol{g}}_{S\mid Z}\left(Z\right)\right)-\tilde{\theta}_1\\
\left(\mu_0\left(X\right)',\mu_1\left(X\right)',1\right)'\left(P^\ast\left({\boldsymbol{g}}_{S\mid Z}\left(Z\right),Z\right)-{\boldsymbol{g}}_{S\mid Z}\left(Z\right)\right)-\tilde{\theta}_2\\
{\boldsymbol{g}}_{U\mid P}\left(P^\ast\left({\boldsymbol{g}}_{S\mid Z}\left(Z\right),Z\right)\right)-\mathcal{U}\left(W,\tilde{\theta}\right)-\tilde{\theta}_3\\
\end{array}
\right)
\right]
=0,
$$
which correspond to line \eqref{eq:m_1_a} for $m_1$, line \eqref{eq:m_2_a} for $m_2$, and line \eqref{eq:m_3_a} for $m_3$.
However, since we estimate ${\xi}_1$, ${\boldsymbol{g}}_{U\mid P}$, and ${\boldsymbol{g}}_{S\mid Z}$,
influence function adjustments need to be added for these estimators.
First, the influence function adjustment for the estimated ${\xi}_1$ is zero.
This is because, as is well known in the literature on partial linear models \citep{robinson1988root}, the moment condition $E[{\xi}_1\left(X,Y,{\boldsymbol{g}}_{S\mid Z}\left(Z\right)\right)-\tilde{\theta}_1]=0$ satisfies the orthogonality property with respect to ${\xi}_1$.
Second, the influence function adjustment to account for an estimation of ${\boldsymbol{g}}_{U\mid P}$ appears in \eqref{eq:m_3_b}.
An estimation of ${\boldsymbol{g}}_{U\mid P}\left(P\right)$ entails an influence function adjustment by $\mathcal{U}\left(W,{\theta}\right)-{\boldsymbol{g}}_{U\mid P}\left({\boldsymbol{g}}_{S\mid Z}\left(Z\right)\right)$.
Since the function ${\boldsymbol{g}}_{U\mid P}$ is evaluated at $P^\ast\left({\boldsymbol{g}}_{S\mid Z}\left(Z\right),Z\right)$ in line \eqref{eq:m_3_a} for $m_3$, this adjustment is scaled by the coefficient ${{f}_{P^\ast}\left({\boldsymbol{g}}_{S\mid Z}\left(Z\right)\right)}/{{f}_{P}\left({\boldsymbol{g}}_{S\mid Z}\left(Z\right)\right)}$ as in \eqref{eq:m_3_b}.
Third, the influence function adjustments for the estimated ${\boldsymbol{g}}_{S\mid Z}$ appear in line \eqref{eq:m_1_b} for $m_1$, line \eqref{eq:m_2_b} for $m_2$, and line \eqref{eq:m_3_c} for $m_3$.
An estimation of ${\boldsymbol{g}}_{S\mid Z}\left(Z\right)$ entails an influence function adjustment by $S-{\boldsymbol{g}}_{S\mid Z}\left(Z\right)$.
Since the function ${\boldsymbol{g}}_{S\mid Z}$ appears inside other functions in lines \eqref{eq:m_1_a}, \eqref{eq:m_2_a}, \eqref{eq:m_3_a} and \eqref{eq:m_3_b}, these adjustments are scaled by their derivatives of the respective functions, which result in the coefficients in front of $S-{\boldsymbol{g}}_{S\mid Z}\left(Z\right)$ in lines \eqref{eq:m_1_b}, \eqref{eq:m_2_b} and \eqref{eq:m_3_c}.
Note that all these adjustments take the form prescribed by \citet[][Proposition 4]{newey1994asymptotic}.
This is because, as mentioned earlier, our moment functions share two convenient properties.
First, $\gamma$ takes forms of mean-square projections or densities.
Second, $\gamma$ enter the equations for $\theta$ so that the Gateaux derivative operator (and thus the orthogonalization operator) can directly act on $\gamma$.
Formal mathematical analyses rationalizing these intuitions are found in Appendix \ref{sec:orthogonality_lemma}.
Given a random sample $\{W_i\}_{i=1}^n$ of size $n$, we now use the above orthogonal score to construct a double debiased estimator with cross fitting or sample splitting.
Let $L > 1$ be a natural number of folds, and randomly partition the sample index set $\{1,...,n\}$ into $L > 1$ subsets $I_1,...,I_L$ of approximately equal size.
For every subsample index $\ell \in \{1,...,L\}$, let
$$\hat{\gamma}_\ell=
\left(
\frac{\hat{f}_{P^\ast}\left(\hat{\boldsymbol{g}}_{S\mid Z}\right)}{\hat{f}_{P}\left(\hat{\boldsymbol{g}}_{S\mid Z}\right)},
\hat{\boldsymbol{g}}_{S\mid Z},
\hat{\xi}_1,
\hat{\zeta},
\hat{\boldsymbol{g}}_{U\mid P},
\hat{\Delta}_{U\mid P}
\right)$$
denote a preliminary parameter estimate obtained by using all observations $i \in \{1,...,n\} \backslash I_\ell$.\footnote{For the components of the preliminary parameter estimate $\hat{\gamma}_\ell$, we omit $\ell$ for the sake of notational simplicity.}
We define our double debiased estimator $\hat\theta$ of $\theta$ as the solution to
\begin{align}\label{eq:double_debiased_theta}
\frac{1}{L} \sum_{\ell=1}^L \frac{1}{|I_\ell|} \sum_{i \in I_\ell} m\left(W_i; \hat\theta , \hat{\gamma}_\ell\right) = 0.
\end{align}
Accordingly, our double debiased estimator $\widehat{PRTE}$ of the PRTE is defined by
\begin{align}\label{eq:double_debiased_prte}
\widehat{PRTE} = \Lambda\left( \hat\theta \right).
\end{align}
\section{Large Sample Theory for Inference}\label{sec:large_sample}
In this section, we investigate the large sample theory for the double debiased estimators, $\hat\theta$ and $\widehat{PRTE}$, defined in \eqref{eq:double_debiased_theta} and \eqref{eq:double_debiased_prte}, respectively.
To this end, we formally present and discuss relevant assumptions.
First, we consider a counterfactual policy $P^\ast$ satisfying the following conditions.
\begin{assumption}\label{assn:pstar}
The support of $P^\ast$ is a subset of the support of $P$.
The cumulative distribution function of $P^\ast$ is absolutely continuous with the Lebesgue measure.
\end{assumption}
Recall that many treatment parameters can be represented in terms of the MTE in similar manners to the PRTE -- see \citet{heckman/vytlacil:2005}.
Assumption \ref{assn:pstar} rules out some of the parameters such as $ATT$, $ATU$, $ATE$, $LATE\left(x\right)$, and $TUT\left(x\right)$.
Second, assume that all the preliminary parameter components of the orthogonal moment function are identified in the following sense.
\begin{assumption}\label{Assn_identification}
$\gamma$ is identifiable, and the following matrix is invertible:
$$
E\left[
\left(\begin{array}{c}
\left(1-P\right)\left(\mu_0\left(X\right)-\boldsymbol{g}_{\mu_0\left(X\right)\mid P}\left(P\right)\right)\\
P\left(\mu_1\left(X\right)-\boldsymbol{g}_{\mu_1\left(X\right)\mid P}\left(P\right)\right)
\end{array}\right)
\left(\begin{array}{c}
\left(1-P\right)\left(\mu_0\left(X\right)-\boldsymbol{g}_{\mu_0\left(X\right)\mid P}\left(P\right)\right)\\
P\left(\mu_1\left(X\right)-\boldsymbol{g}_{\mu_1\left(X\right)\mid P}\left(P\right)\right)
\end{array}\right)'
\right].
$$
\end{assumption}
The invertibility condition serves to identify the parameter vector $\left(\beta_0',\beta_1'\right)'$ in the partial linear model as in \citet{carneiro/lee:2009}.
The requirement for $\gamma$ to be identified essentially boils down to the identification of the density functions and the conditional expectation functions, which is satisfied for a large class of data generating processes, and is usually taken for granted in the nonparametrics literature.
We next make three high-level conditions (Assumptions \ref{Assn_slow_conver}--\ref{Assn:Taylor_Errr}) regarding large sample behaviors of the preliminary parameter estimators.
While we stress that each of these three assumptions is a mild requirement for most common nonparametric estimators, we supplement each of the three non-primitive assumptions by lower-level sufficient conditions in Appendix \ref{sec:lowe_level_sufficient_conditions}.
\begin{assumption}\label{Assn_slow_conver}
For every $\hat{\gamma}_\ell=\hat{\gamma}_1,\ldots,\hat{\gamma}_L$, the following objects are $o_p(1)$.\footnotesize
\begingroup
\allowdisplaybreaks
\begin{align*}
&\int \left(\hat{\zeta}\left(z\right)-{\zeta}\left(z\right)\right)\left(\hat{\boldsymbol{g}}_{S\mid Z}\left(z\right)-{\boldsymbol{g}}_{S\mid Z}\left(z\right)\right)F_W\left(dw\right).
\\
&\int\frac{\hat{f}_{P^\ast}\left(\hat{\boldsymbol{g}}_{S\mid Z}\left(z\right)\right)}{\hat{f}_{P}\left(\hat{\boldsymbol{g}}_{S\mid Z}\left(z\right)\right)}\left(\hat{\Delta}_{U\mid P}\left(\hat{\boldsymbol{g}}_{S\mid Z}\left(z\right)\right)-{\Delta}_{U\mid P}\left({\boldsymbol{g}}_{S\mid Z}\left(z\right)\right)
\left({\boldsymbol{g}_{S\mid Z}}\left(z\right)-\hat{\boldsymbol{g}}_{S\mid Z}\left(z\right)\right)\right)F_W\left(dw\right)
\\
&\int \left(\hat{\Delta}_{U\mid P}\left(P^\ast\left(\hat{\boldsymbol{g}}_{S\mid Z}\left(z\right),z\right)\right)\partial P^\ast\left(\hat{\boldsymbol{g}}_{S\mid Z}\left(z\right),z\right)-{\Delta}_{U\mid P}\left(P^\ast\left({\boldsymbol{g}}_{S\mid Z}\left(z\right),z\right)\right)\partial P^\ast\left({\boldsymbol{g}}_{S\mid Z}\left(z\right),z\right)\right)
\left({\boldsymbol{g}_{S\mid Z}}\left(z\right)-\hat{\boldsymbol{g}}_{S\mid Z}\left(z\right)\right)F_W\left(dw\right).
\end{align*}
\endgroup\normalsize
\end{assumption}
Assumption \ref{Assn_slow_conver} requires that some preliminary parameter estimators to have the $n^{1/4}$ rate of convergence.
This slow rate of convergence suffices because of the orthogonal score we developed and proposed in Section \ref{sec:orthogonal_score} -- see the intuitions and discussions in Section \ref{sec:intuitions_discussions}.
As far as this condition is satisfied, the asymptotic distribution of our parameters of interest, namely $\hat\theta$ and $\widehat{PRTE}$, will not be affected by estimation errors of these preliminary parameter estimates.
Furthermore, this condition can be satisfied by common nonparametric estimators of density and conditional expectation functions.
For completeness, we show lower-level sufficient conditions for Assumption \ref{Assn_slow_conver} in Appendix \ref{sec:lowe_level_sufficient_conditions}.
\begin{assumption}\label{Assn_CV}
The following objects are $o_p(1)$.\footnotesize
\begingroup
\allowdisplaybreaks
\begin{align*}
&\int
\left(\hat{\boldsymbol{g}}_{U\mid P}\left(P^\ast\left(\hat{\boldsymbol{g}}_{S\mid Z}\left(z\right),z\right)\right)-{\boldsymbol{g}}_{U\mid P}\left(P^\ast\left(\hat{\boldsymbol{g}}_{S\mid Z}\left(z\right),z\right)\right)\right)
F_W\left(dw\right)
-\int
\frac{\hat{f}_{P^\ast}\left(\hat{\boldsymbol{g}}_{S\mid Z}\left(z\right)\right)}{\hat{f}_{P}\left(\hat{\boldsymbol{g}}_{S\mid Z}\left(z\right)\right)}
\left(\hat{\boldsymbol{g}}_{U\mid P}\left(\hat{\boldsymbol{g}}_{S\mid Z}\left(z\right)-{\boldsymbol{g}}_{U\mid P}\left(\hat{\boldsymbol{g}}_{S\mid Z}\left(z\right)\right)\right)\right)
F_W\left(dw\right)
\\
&\int \left(\hat{\xi}_1\left(x,y,\hat{\boldsymbol{g}}_{S\mid Z}\left(z\right)\right)-{\xi}_1\left(x,y,\hat{\boldsymbol{g}}_{S\mid Z}\left(z\right)\right)\right)F_W\left(dw\right)
\end{align*}
\endgroup\normalsize
\end{assumption}
We show lower-level sufficient conditions for Assumption \ref{Assn_CV} in Appendix \ref{sec:lowe_level_sufficient_conditions}.
The final piece of assumptions is to control higher-order effects of the propensity score estimation on the orthogonal score.
Note that the orthogonality property \eqref{eq:orthogonality} forces the first-order effects of nuisance parameter estimation to be zero.
Since all the nuisance parameters except for the propensity score function $\boldsymbol{g}_{S|Z}$ appear in our score in a linear manner, the orthogonality property precisely vanishes their estimation effects.
On the other hand, since the propensity score function $\boldsymbol{g}_{S|Z}$ appears in our score nonlinearly in general, we need to make sure that the remaining higher-order effects are small enough.
Specifically, we require the following condition.
\begin{assumption}\label{Assn:Taylor_Errr}
$\int\boldsymbol{R}_k\left(z\right)F_Z\left(dz\right)=o_p\left(n^{-1/2}\right)$ for every $\hat{\gamma}_\ell$ and every $k=1,\ldots,4$, where \footnotesize
\begingroup
\allowdisplaybreaks
\begin{eqnarray*}
\boldsymbol{R}_1\left(z\right)&=&
\int\left({\xi}_1\left(x,y,\hat{\boldsymbol{g}}_{S\mid Z}\left(z\right)\right)-{\xi}_1\left(x,y,{\boldsymbol{g}}_{S\mid Z}\left(z\right)\right)\right)F_{\left(Y,X\right)\mid Z=z}\left(dy,dx\right)-{\zeta}\left(z\right)\left(\hat{\boldsymbol{g}}_{S\mid Z}\left(z\right)-{\boldsymbol{g}}_{S\mid Z}\left(z\right)\right)
\\
\boldsymbol{R}_2\left(z\right)&=&
\boldsymbol{g}_{\left(\mu_0\left(X\right)',\mu_1\left(X\right)',1\right)'\mid Z}\left(z\right)\left(P^\ast\left(\hat{\boldsymbol{g}}_{S\mid Z}\left(z\right),z\right)-P^\ast\left({\boldsymbol{g}_{S\mid Z}}\left(z\right),z\right)-\partial P^\ast\left(\hat{\boldsymbol{g}}_{S\mid Z}\left(z\right),z\right)
\left(\hat{\boldsymbol{g}}_{S\mid Z}\left(z\right)-\boldsymbol{g}_{S\mid Z}\left(z\right)\right)\right)
\\
\boldsymbol{R}_3\left(z\right)&=&
{\boldsymbol{g}}_{U\mid P}\left(P^\ast\left(\hat{\boldsymbol{g}}_{S\mid Z}\left(z\right),z\right)\right)-{\boldsymbol{g}_{U\mid P}}\left(P^\ast\left({\boldsymbol{g}_{S\mid Z}}\left(z\right),z\right)\right)
-
{\Delta}_{U\mid P}\left(P^\ast\left({\boldsymbol{g}}_{S\mid Z}\left(z\right),z\right)\right)\partial P^\ast\left({\boldsymbol{g}}_{S\mid Z}\left(z\right),z\right)
\left(\hat{\boldsymbol{g}}_{S\mid Z}\left(z\right)-{\boldsymbol{g}}_{S\mid Z}\left(z\right)\right)
\\
\boldsymbol{R}_4\left(z\right)&=&
\frac{\hat{f}_{P^\ast}\left(\hat{\boldsymbol{g}}_{S\mid Z}\left(z\right)\right)}{\hat{f}_{P}\left(\hat{\boldsymbol{g}}_{S\mid Z}\left(z\right)\right)}
\left({\boldsymbol{g}_{U\mid P}}\left({\boldsymbol{g}_{S\mid Z}}\left(z\right)\right)-{\boldsymbol{g}}_{U\mid P}\left(\hat{\boldsymbol{g}}_{S\mid Z}\left(z\right)\right)
-{\Delta}_{U\mid P}\left({\boldsymbol{g}}_{S\mid Z}\left(z\right)\right)
\left({\boldsymbol{g}_{S\mid Z}}\left(z\right)-\hat{\boldsymbol{g}}_{S\mid Z}\left(z\right)\right)\right)
\end{eqnarray*}
\endgroup\normalsize
\end{assumption}
Since higher-order effects are usually of smaller magnitudes, Assumption \ref{Assn:Taylor_Errr} is a quite mild condition.
We present lower-level sufficient conditions for Assumption \ref{Assn:Taylor_Errr} in Appendix \ref{sec:lowe_level_sufficient_conditions}.
In addition to these conditions, we also present regularity conditions (Assumptions \ref{Assn_compact}, \ref{Assn_bounded_moments} and \ref{Assn_just_consis}) in Appendix \ref{sec:regularity}.
These additional conditions are relegated to the appendix for the sake of readability, as they are even more standard, even milder, and are even easier to verify with common nonparametric estimators than Assumptions \ref{Assn_slow_conver}--\ref{Assn:Taylor_Errr}.
We obtain the following asymptotic normality result for the double debiased estimator $\hat\theta$ of the intermediate parameter vector.
\begin{theorem}\label{theorem:normal}
If Assumptions \ref{assn:pstar}, \ref{Assn_identification}, \ref{Assn_slow_conver}, \ref{Assn:Taylor_Errr}, \ref{Assn_compact}, \ref{Assn_bounded_moments} and \ref{Assn_just_consis} are satisfied, then
$$
\sqrt{n}\left(\hat\theta-\theta\right)\rightarrow_dN\left(0,\left(\mathcal{M}'\mathcal{M}\right)^{-1}\mathcal{M}'E[m\left(W;{\theta},{\gamma}\right)m\left(W;{\theta},{\gamma}\right)']\mathcal{M}\left(\mathcal{M}'\mathcal{M}\right)^{-1}\right),
$$
where, with $m_{3,2}$ defined in Appendix and \ref{sec:decomposition_of_moment_functions}, the matrix $\mathcal{M}$ takes the form
$$
\mathcal{M}=\left(
\begin{array}{ccc}
I&0&0\\
0&I&0\\
E\left[m_{3,2}\left(W;{\gamma}\right)\right]\frac{\partial}{\partial\theta_1'}\boldsymbol{d}\left(\theta_1\right)&0&I
\end{array}
\right).
$$\end{theorem}
A proof is provided in Appendix \ref{sec:theorem:normal}.
An immediate consequence of this result through the delta method is the following asymptotic normality result for the double debiased estimator $\widehat{PRTE}$ of the PRTE.
\begin{corollary}\label{cor:clt}
If Assumptions \ref{assn:pstar}, \ref{Assn_identification}, \ref{Assn_slow_conver}, \ref{Assn:Taylor_Errr}, \ref{Assn_compact}, \ref{Assn_bounded_moments} and \ref{Assn_just_consis} are satisfied, then
$$
\sqrt{n}\left(\widehat{PRTE}-PRTE\right)\rightarrow_dN\left(0,\lambda\left(\theta\right)\left(\mathcal{M}'\mathcal{M}\right)^{-1}\mathcal{M}'E[m\left(W;{\theta},{\gamma}\right)m\left(W;{\theta},{\gamma}\right)']\mathcal{M}\left(\mathcal{M}'\mathcal{M}\right)^{-1}\lambda\left(\theta\right)'\right).
$$
where $\lambda$ denotes the derivative of $\Lambda$.
\end{corollary}
As we emphasized throughout this paper, this root-$n$ convergence without any influence of preliminary estimation errors $\hat\gamma-\gamma$ is possible as a result of the orthogonality that implies \eqref{eq:orthogonality}.
This method can accommodate a wide array of preliminary estimation techniques including kernel-smoothing, sieve estimation, and shrinkage methods among others.
This flexibility in the choice of preliminary estimation approaches at the current high-level theory is another advantage of using the orthogonal score, unlike conventional methods in absence of orthogonal scores.
A drawback of this root-$n$ asymptotic normality result is that it is not guaranteed to be efficient.
We can estimate every component of the asymptotic variance for $\widehat{PRTE}$ as follows.
First, we can estimate $\lambda\left(\theta\right)$ by $\lambda\left(\hat\theta\right)$ since $\lambda$ is a known function.
Second, estimate $\mathcal{M}$ by
$$
\hat{\mathcal{M}}=\left(
\begin{array}{ccc}
I&0&0\\
0&I&0\\
\frac{1}{n}\sum_{i=1}^nm_{3,2}\left(W_i;\hat{\gamma}\right)\frac{\partial}{\partial\theta_1'}\boldsymbol{d}\left(\hat\theta_1\right)&0&I
\end{array}
\right).
$$
Last, estimate $E[m\left(W;{\theta},{\gamma}\right)m\left(W;{\theta},{\gamma}\right)']$ by
\begin{align*}
\hat\Sigma = \frac{1}{n}\sum_{i=1}^n
\left[\begin{array}{c}
\left(\begin{array}{c}
\hat m_1\left(W_i;\hat{\theta},\hat{\gamma}\right)\\\hat m_2\left(W_i;\hat{\theta},\hat{\gamma}\right)\\\hat m_3\left(W_i;\hat{\theta},\hat{\gamma}\right)
\end{array}\right)
\left(\begin{array}{c}
\hat m_1\left(W_i;\hat{\theta},\hat{\gamma}\right)\\\hat m_2\left(W_i;\hat{\theta},\hat{\gamma}\right)\\\hat m_3\left(W_i;\hat{\theta},\hat{\gamma}\right)
\end{array}\right)'
\end{array}\right].
\end{align*}
The following theorem guarantees the asymptotic validity of the resultant variance estimator.
\begin{theorem}\label{theorem:var_est}
If Assumptions \ref{assn:pstar}, \ref{Assn_identification}, \ref{Assn_slow_conver}, \ref{Assn:Taylor_Errr}, \ref{Assn_compact}, \ref{Assn_bounded_moments} and \ref{Assn_just_consis} are satisfied, then
$$
\lambda\left(\hat\theta\right)\left(\hat{\mathcal{M}}'\hat{\mathcal{M}}\right)^{-1}\hat{\mathcal{M}}'\hat\Sigma\hat{\mathcal{M}}\left(\hat{\mathcal{M}}'\hat{\mathcal{M}}\right)^{-1}\lambda\left(\hat\theta\right)'
$$
is a consistent estimator of the asymptotic variance for $\widehat{PRTE}$.
\end{theorem}
A proof is provided in Appendix \ref{sec:theorem:variance_est}.
\section{Monte Carlo Simulations}\label{sec:simulations}
In this section, we use Monte Carlo simulations to evaluate finite sample performance of the proposed method of double debiased estimation and inference.
We first consider a benchmark design from the literature in Section \ref{sec:benchmark_design}, and demonstrate that even completely nonparametric preliminary estimation can produce desirable finite-sample performance thanks to the orthogonal score.
We second extend the design with higher dimensions of $Z$ in Section \ref{sec:extended_designs}.
\subsection{Nonparametric Preliminary Estimation under a Benchmark Design}\label{sec:benchmark_design}
We generate artificial datasets from the same distribution as in \citet[][supplementary Appendix D]{carneiro/lokshin/umapathi:2017}.
Specifically, we generate
$$
\left(U_0,U_1,U_S\right)=\left(-0.050 \varepsilon_1 + 0.020 \varepsilon_3, 0.012 \varepsilon_1 + 0.010 \varepsilon_2,-1.000 \varepsilon_1\right),
$$
where $\varepsilon_1$ and $\varepsilon_2$ are independent standard normal random variables.
Then, we generate
\begingroup
\allowdisplaybreaks
\begin{align*}
Y_1 =& 0.240 + 0.800 X_1 + 0.400 X_2 + U_1
\\
Y_0 =& 0.020 + 0.500 X_1 + 0.100 X_2 + U_0,
\qquad\text{and}\\
S =& \mathbbm{1}\{ {0.200 + 0.300 Z_1 + 0.100 Z_2} - U_S > 0 \},
\end{align*}
\endgroup
where $X_1 \sim N\left(-2,2^2\right)$, $X_2 \sim N\left(2,2^2\right)$, $Z_1 \sim N\left(-1,3^2\right)$, and $Z_2 \sim N\left(1,3^2\right)$ are mutually independent and are also independent of $\left(U_1,U_0,U_S\right)'$.
The observed outcome is generated in turn as $Y = SY_1 + \left(1-S\right)Y_0$.
We consider two forms of policy changes.
The first takes the form of $P^\ast = P + a\left(1-P\right)$, where $a \in \{0.1,...,0.9\}$.\footnote{Note that $PRTE \rightarrow MPRTE$ \citep[marginal PRTE,][]{carneiro/heckman/vytlacil:2010} as $a \rightarrow 0$ and $PRTE \rightarrow ATU$ as $a \rightarrow 1$.}
The second takes the form of $Z^\ast = Z + a$, where $a \in \{-0.5,...,-0.1,0.1,...,0.5\}$.
See Appendix \ref{sec:details_simulation_setting} for additional details about this simulation setting, including implied analytic expressions for the counterfactual distribution, marginal treatment effects, and the PRTEs.
For each instance of artificial datasets, we estimate the PRTE and its estimated standard error using our proposed method of double debiased estimation and inference with completely nonparametric preliminary estimation.
See Appendix \ref{sec:details_simulation_estimation_p} and \ref{sec:details_simulation_estimation_z} for additional details about concrete estimation and inference procedures.
We experiment with alternative numbers, $L=5$ and $10$, of folds in cross fitting, where the number of observations in each fold is equal.
We report the simulated bias, root mean square error, and coverage frequencies for the nominal probability of 95\%.
Simulated coverage frequencies are computed based on the standard symmetric confidence interval generated by the PRTE estimate and its estimated standard error according to Corollary \ref{cor:clt}.
The number of Monte Carlo iterations is set to 1000 following \citet[][supplementary Appendix D]{carneiro/lokshin/umapathi:2017}.
Since there is no existing limit distribution theory for inference about the PRTE to our best knowledge, we do not have any benchmark in the literature against which to make comparisons of our results.
We therefore present simulation results only for our proposed method.
Tables \ref{tab:simulation_results_p} and \ref{tab:simulation_results_z} summarize simulation results for policy changes of the forms $P^\ast = P + a\left(1-P\right)$ and $Z^\ast = Z + a$, respectively.
The results concerning the mean and bias demonstrate that the proposed double debiased estimator indeed produces small biases (relative to the root mean square error) even in small samples.
In light of these relatively small biases, the results for the root mean square error demonstrate that the estimator converges approximately at the rate of $n^{-1/2}$, consistently with our theory.
The results concerning the coverage frequencies demonstrate that our limit normal distribution results are useful to construct asymptotically valid confidence intervals for the PRTE, except when $a$ is infinitesimal.
The alternative numbers, $L=5$ and $10$, of folds for sample splitting entail quite similar simulation results in terms of all of the bias, root mean square error, and coverage frequencies.
In summary, the proposed estimator enjoys the main properties suggested in this paper even in small samples; namely, it is debiased, converges at the parametric rate, and is asymptotically normal with the proposed asymptotic variance formula, regardless of the number of folds in sample splitting.
\begin{table}
\centering
\scalebox{1.00}{
\begin{tabular}{cccccccc}
\hline\hline
& & & True & \multicolumn{3}{c}{Estimates}& 95\%\\
\cline{5-7}
$a$ & $n$ & $L$ & PRTE & Mean & Bias & RMSE & Coverage\\
\hline
0.1 & 1000 & 5 & 0.243 & 0.324 & 0.081 & 0.889 & 0.828\\
& 2000 & 5 & 0.243 & 0.298 & 0.054 & 0.540 & 0.828\\
\hline
0.2 & 1000 & 5 & 0.237 & 0.309 & 0.073 & 0.516 & 0.863\\
& 2000 & 5 & 0.237 & 0.285 & 0.048 & 0.301 & 0.919\\
\hline
0.3 & 1000 & 5 & 0.230 & 0.320 & 0.090 & 0.420 & 0.901\\
& 2000 & 5 & 0.230 & 0.278 & 0.048 & 0.241 & 0.957\\
\hline
0.4 & 1000 & 5 & 0.225 & 0.298 & 0.073 & 0.376 & 0.912\\
& 2000 & 5 & 0.225 & 0.272 & 0.047 & 0.218 & 0.947\\
\hline
0.5 & 1000 & 5 & 0.219 & 0.287 & 0.068 & 0.347 & 0.922\\
& 2000 & 5 & 0.219 & 0.256 & 0.037 & 0.205 & 0.957\\
\hline
0.6 & 1000 & 5 & 0.213 & 0.273 & 0.060 & 0.323 & 0.928\\
& 2000 & 5 & 0.213 & 0.247 & 0.034 & 0.200 & 0.951\\
\hline
0.7 & 1000 & 5 & 0.207 & 0.259 & 0.052 & 0.312 & 0.923\\
& 2000 & 5 & 0.207 & 0.232 & 0.025 & 0.195 & 0.949\\
\hline
0.8 & 1000 & 5 & 0.201 & 0.245 & 0.044 & 0.306 & 0.908\\
& 2000 & 5 & 0.201 & 0.219 & 0.018 & 0.195 & 0.935\\
\hline
0.9 & 1000 & 5 & 0.194 & 0.223 & 0.028 & 0.311 & 0.872\\
& 2000 & 5 & 0.194 & 0.207 & 0.013 & 0.206 & 0.892\\
\hline
& & & True & \multicolumn{3}{c}{Estimates}& 95\%\\
\cline{5-7}
$a$ & $n$ & $L$ & PRTE & Mean & Bias & RMSE & Coverage\\
\hline
0.1 & 1000 & 10 & 0.243 & 0.321 & 0.078 & 0.876 & 0.829\\
& 2000 & 10 & 0.243 & 0.299 & 0.055 & 0.503 & 0.874\\
\hline
0.2 & 1000 & 10 & 0.237 & 0.295 & 0.058 & 0.513 & 0.878\\
& 2000 & 10 & 0.237 & 0.288 & 0.051 & 0.299 & 0.916\\
\hline
0.3 & 1000 & 10 & 0.230 & 0.308 & 0.077 & 0.405 & 0.904\\
& 2000 & 10 & 0.230 & 0.286 & 0.056 & 0.246 & 0.948\\
\hline
0.4 & 1000 & 10 & 0.225 & 0.293 & 0.068 & 0.362 & 0.917\\
& 2000 & 10 & 0.225 & 0.280 & 0.055 & 0.220 & 0.956\\
\hline
0.5 & 1000 & 10 & 0.219 & 0.285 & 0.066 & 0.340 & 0.929\\
& 2000 & 10 & 0.219 & 0.270 & 0.051 & 0.208 & 0.954\\
\hline
0.6 & 1000 & 10 & 0.213 & 0.271 & 0.058 & 0.321 & 0.929\\
& 2000 & 10 & 0.213 & 0.258 & 0.044 & 0.196 & 0.954\\
\hline
0.7 & 1000 & 10 & 0.207 & 0.259 & 0.052 & 0.307 & 0.924\\
& 2000 & 10 & 0.207 & 0.243 & 0.036 & 0.190 & 0.943\\
\hline
0.8 & 1000 & 10 & 0.201 & 0.240 & 0.039 & 0.308 & 0.908\\
& 2000 & 10 & 0.201 & 0.230 & 0.029 & 0.190 & 0.932\\
\hline
0.9 & 1000 & 10 & 0.194 & 0.218 & 0.023 & 0.307 & 0.894\\
& 2000 & 10 & 0.194 & 0.220 & 0.026 & 0.198 & 0.898\\
\hline\hline
\end{tabular}
}
\caption{Monte Carlo simulation results for the benchmark design with completely nonparametric preliminary estimation under policy changes of the form $P^\ast = P + a \left(1 - P\right)$ for $a \in \{0.1,...,0.9\}$. $n$ denotes the sample size. $L$ is the number of folds in sample splitting. The true PRTE is numerically evaluated. The displayed statistics include the mean, bias, root mean square error, and coverage frequencies with the nominal probability of 95\%.}
\label{tab:simulation_results_p}
\end{table}
\begin{table}
\centering
\scalebox{1.00}{
\begin{tabular}{cccccccc}
\hline\hline
& & & True & \multicolumn{3}{c}{Estimates}& 95\%\\
\cline{5-7}
$a$ & $n$ & $L$ & PRTE & Mean & Bias & RMSE & Coverage\\
\hline
-0.5& 1000 & 5 & 0.223 & 0.187 &-0.035 & 0.166 & 0.946\\
& 2000 & 5 & 0.223 & 0.190 &-0.033 & 0.113 & 0.950\\
\hline
-0.4& 1000 & 5 & 0.222 & 0.186 &-0.036 & 0.167 & 0.945\\
& 2000 & 5 & 0.222 & 0.188 &-0.035 & 0.113 & 0.947\\
\hline
-0.3& 1000 & 5 & 0.222 & 0.184 &-0.038 & 0.169 & 0.942\\
& 2000 & 5 & 0.222 & 0.186 &-0.036 & 0.114 & 0.948\\
\hline
-0.2& 1000 & 5 & 0.221 & 0.183 &-0.038 & 0.171 & 0.933\\
& 2000 & 5 & 0.221 & 0.184 &-0.037 & 0.115 & 0.945\\
\hline
-0.1& 1000 & 5 & 0.221 & 0.182 &-0.039 & 0.177 & 0.927\\
& 2000 & 5 & 0.221 & 0.182 &-0.039 & 0.116 & 0.940\\
\hline
0.1 & 1000 & 5 & 0.220 & 0.202 &-0.018 & 0.178 & 0.924\\
& 2000 & 5 & 0.220 & 0.186 &-0.034 & 0.121 & 0.930\\
\hline
0.2 & 1000 & 5 & 0.219 & 0.206 &-0.014 & 0.177 & 0.923\\
& 2000 & 5 & 0.219 & 0.187 &-0.032 & 0.118 & 0.934\\
\hline
0.3 & 1000 & 5 & 0.219 & 0.207 &-0.012 & 0.176 & 0.930\\
& 2000 & 5 & 0.219 & 0.186 &-0.032 & 0.117 & 0.940\\
\hline
0.4 & 1000 & 5 & 0.218 & 0.207 &-0.012 & 0.176 & 0.927\\
& 2000 & 5 & 0.218 & 0.186 &-0.032 & 0.116 & 0.946\\
\hline
0.5 & 1000 & 5 & 0.218 & 0.206 &-0.012 & 0.174 & 0.934\\
& 2000 & 5 & 0.218 & 0.185 &-0.033 & 0.115 & 0.949\\
\hline\hline
\end{tabular}
}
\caption{Monte Carlo simulation results for the benchmark design with completely nonparametric preliminary estimation under policy changes of the form $Z^\ast = Z + a$ for $a \in \{-0.5,...,-0.1,0.1,...,0.5\}$. $n$ denotes the sample size. $L$ is the number of folds in sample splitting. The true PRTE is numerically evaluated. The displayed statistics include the mean, bias, root mean square error, and coverage frequencies with the nominal probability of 95\%.}
\label{tab:simulation_results_z}
\end{table}
\clearpage
\subsection{Extended Designs with High Dimensions}\label{sec:extended_designs}
In the previous subsection, we demonstrate that even completely nonparametric preliminary estimation yields desirable finite-sample performances.
We next consider extended designs with high dimensions, employ parametric propensity score estimation with absolute shrinkage penalization, and present finite-sample performances of our estimation and inference procedures in this setting.
Suppose that the selection model is specified by
\begingroup
\allowdisplaybreaks
\begin{align*}
S =& \mathbbm{1}\left\{ \frac{3}{10} - \sum_{k=1}^{\text{dim}\left(Z\right)-2} \frac{3}{10^{\left(k+2\right)/2}} - \frac{1}{10^{{\text{dim}\left(Z\right)}/2}} \right.
\\
&\ \
\left. + \frac{3}{10}Z_1 + \sum_{k=1}^{\text{dim}\left(Z\right)-2} \frac{3}{10^{\left(k+2\right)/2}} Z_{k+1} + \frac{1}{10^{\text{dim}\left(Z\right)/2}} Z_{\text{dim}\left(Z\right)} - U_S > 0 \right\},
\end{align*}
\endgroup
where $Z_1 \sim N\left(-1,3^2\right)$ and $Z_2, ..., Z_{\text{dim}\left(Z\right)} \sim N\left(1,3^2\right)$.
This design includes that of the benchmark by \citet[][supplementary Appendix D]{carneiro/lokshin/umapathi:2017} studied in the previous subsection as a special case where $\text{dim}\left(Z\right) = 2$.
Furthermore, this design preserves the same value of PRTE as that in the benchmark by \citet{carneiro/lokshin/umapathi:2017} for each $a$, regardless of the dimension $\text{dim}\left(Z\right) \ge 2$.
In this manner, we devise this extended design to allow for comparisons of simulation results with those presented in Section \ref{sec:benchmark_design}.
We set the high-dimensional setting with $\text{dim}\left(Z\right) = 100$ in this design.
Estimates for the PRTE and their estimated standard errors are obtained using our proposed method of double debiased estimation and inference with preliminary estimation of propensity scores by lasso logit.
See Appendix \ref{sec:details_simulation_parametric_estimation_p} for additional details.
Table \ref{tab:simulation_results_lasso_logit} summarizes simulation results.
Note that we obtain qualitatively similar results to those presented in Table \ref{tab:simulation_results_p}, and hence similar remarks follow.
It is worth noting that the magnitude of the statistics (such as the bias the RMSE) displayed in this table is similar to that in Table \ref{tab:simulation_results_p} despite the difference in $\text{dim}\left(Z\right)$ and the preliminary estimators.
This observation is also consistent with the property of the orthogonal score that it reduces and asymptotically vanishes effects of preliminary estimation errors.
\begin{table}
\centering
\scalebox{0.95}{
\begin{tabular}{cccccccc}
\hline\hline
& & & True & \multicolumn{3}{c}{Estimates}& 95\%\\
\cline{5-7}
$a$ & $n$ & $L$ & PRTE & Mean & Bias & RMSE & Coverage\\
\hline
0.1 & 1000 & 5 & 0.243 & 0.306 & 0.063 & 0.783 & 0.907\\
& 2000 & 5 & 0.243 & 0.294 & 0.051 & 0.457 & 0.919\\
\hline
0.2 & 1000 & 5 & 0.237 & 0.319 & 0.083 & 0.499 & 0.935\\
& 2000 & 5 & 0.237 & 0.289 & 0.052 & 0.306 & 0.942\\
\hline
0.3 & 1000 & 5 & 0.230 & 0.327 & 0.096 & 0.422 & 0.935\\
& 2000 & 5 & 0.230 & 0.301 & 0.070 & 0.256 & 0.958\\
\hline
0.4 & 1000 & 5 & 0.225 & 0.314 & 0.090 & 0.385 & 0.940\\
& 2000 & 5 & 0.225 & 0.281 & 0.057 & 0.242 & 0.949\\
\hline
0.5 & 1000 & 5 & 0.219 & 0.310 & 0.091 & 0.382 & 0.931\\
& 2000 & 5 & 0.219 & 0.277 & 0.059 & 0.226 & 0.945\\
\hline
0.6 & 1000 & 5 & 0.213 & 0.281 & 0.068 & 0.358 & 0.929\\
& 2000 & 5 & 0.213 & 0.259 & 0.045 & 0.222 & 0.934\\
\hline
0.7 & 1000 & 5 & 0.207 & 0.268 & 0.061 & 0.355 & 0.911\\
& 2000 & 5 & 0.207 & 0.247 & 0.040 & 0.208 & 0.935\\
\hline
0.8 & 1000 & 5 & 0.201 & 0.237 & 0.036 & 0.340 & 0.911\\
& 2000 & 5 & 0.201 & 0.223 & 0.022 & 0.212 & 0.921\\
\hline
0.9 & 1000 & 5 & 0.194 & 0.219 & 0.025 & 0.378 & 0.860\\
& 2000 & 5 & 0.194 & 0.215 & 0.021 & 0.226 & 0.869\\
\hline\hline
\end{tabular}
}
\caption{Monte Carlo simulation results for the extended design with $\text{dim}\left(Z\right)=100$ based on lasso logit propensity score estimation under policy changes of the form $P^\ast = P + a\left(1-P\right)$ for $a \in \{0.1,...,0.9\}$. $n$ denotes the sample size. $L$ is the number of folds in sample splitting. The true PRTE is numerically evaluated. The displayed statistics include the mean, bias, root mean square error, and coverage frequencies with the nominal probability of 95\%.}
\label{tab:simulation_results_lasso_logit}
\end{table}
\clearpage
\clearpage
\section{Empirical Illustration}\label{sec:empirical}
In this section, we apply the proposed method to an analysis of the effects of counterfactual policies that expose various fractions of the population to upper secondary schooling in Indonesia.
Following the study by \citet{carneiro/lokshin/umapathi:2017}, we use the data from the third wave of the Indonesia Family Life Survey (IFLS) fielded from June through
November 2000 -- we refer readers to their supplementary appendix for further details of this data set.
We also use the same subsample as that of \citet{carneiro/lokshin/umapathi:2017}, that consists of employed males aged 25--60, who have reported non-missing wage and schooling information.
This subsample consists of 2,608 individuals.
We set our variables following \citet{carneiro/lokshin/umapathi:2017}.
The outcome variable $Y$ denotes the log of hourly wages constructed from self-reported monthly wages and hours worked per week.
The binary treatment variable $S$ indicates attendance of upper secondary school or higher.
Control variables $X$ include age, age squared, an indicator for whether the individual was living in a village at age 12, indicators for the province of residence, an indicator of rural residence, distance from the office of the head of the community of residence to the nearest community health post, and indicators for the level of schooling by each parent.
The excluded instrument is the distance from the office of the head of the community of residence in kilometers to the nearest secondary school.
\citet[][Section 4.1]{carneiro/lokshin/umapathi:2017} provide detailed analysis to support the validity of this instrument.
In this paper, we take advantage of their analysis and discussions, and directly adopt their empirical approach in our estimation and inference framework.
See Appendix \ref{sec:model_application} for details of this model specification.
We consider the counterfactual policies of the form $P^\ast = P + a \left(1-P\right)$ with various levels of $a \in \{0.05,0.10,...,0.45,0.50\}$.
Note that such a counterfactual treatment probability $P^\ast$ arises from the counterfactual policy that exposes fraction $a$ of the population to upper secondary schooling.
The estimation and inference approaches that we take in this analysis are the same as those used for our Monte Carlo simulation studies except that we use the probit propensity score estimation following the benchmark study by \citet{carneiro/lokshin/umapathi:2017} -- see Appendix \ref{sec:step_by_step_procedure_application} for a step-by-step procedure of estimation and inference.
We use $L=5$ as the number of folds in cross fitting.\footnote{We also ran another set of estimation with $L=10$ but the results are almost the same both in terms of estimates and their standard errors, similarly to what we observed in our Monte Carlo simulation studies.}
Finally, to mitigate the finite-sample randomness in estimates due to sample splitting, we use the robust re-randomization method following \citet[][Section 3.4]{chernozhukov2018double}.
Table \ref{tab:empirical} summarizes the results.
For each $a \in \{0.05,0.10,...,0.45,0.50\}$, displayed in this table are the point estimate, standard error, and the 95\% confidence interval.
Observe that the policy that assigns a fraction $a=0.05$ of the population to upper secondary schooling is expected to increase the log of hourly wages by 0.169 (with the standard error of 0.064 and the 95\% confidence interval of $[0.044,0.294]$) per treated individual on average.
As the fraction $a$ of the population assigned to treatment increases, these per-individual average policy relevant treatment effects tend to increase.
Specifically, the policy that assigns a fraction $a=0.50$ of the population to upper secondary schooling is expected to increase the log of hourly wages by 0.225 (with the standard error of 0.076 and the 95\% confidence interval of $[0.076,0.375]$) per treated individual on average.
\begin{table}
\centering
\begin{tabular}{ccccc}
\hline\hline
$a$ & $L$ & Estimate & Std. Err. & 95\% CI\\
\hline
0.05 & 5 & 0.169 & 0.064 & [0.044, 0.294]\\
0.10 & 5 & 0.172 & 0.064 & [0.046, 0.298]\\
0.15 & 5 & 0.176 & 0.064 & [0.049, 0.302]\\
0.20 & 5 & 0.180 & 0.065 & [0.053, 0.308]\\
0.25 & 5 & 0.186 & 0.066 & [0.056, 0.316]\\
0.30 & 5 & 0.192 & 0.068 & [0.059, 0.325]\\
0.35 & 5 & 0.199 & 0.070 & [0.062, 0.336]\\
0.40 & 5 & 0.207 & 0.072 & [0.066, 0.347]\\
0.45 & 5 & 0.216 & 0.074 & [0.071, 0.360]\\
0.50 & 5 & 0.225 & 0.076 & [0.076, 0.375]\\
\hline\hline
\end{tabular}
\caption{Estimates, standard errors, and 95\% confidence intervals of the PRTE for the counterfactual policies of the form $P^\ast = P + a \left(1-P\right)$, $a \in \{0.05,0.10,...,0.45,0.50\}$.}
\label{tab:empirical}
\end{table}
\section{Summary}\label{sec:summary}
Estimation of the PRTE involves estimation of multiple preliminary parameters.
These preliminary parameters include propensity scores, conditional expectation functions of the outcome and covariates given the propensity score, and marginal treatment effects.
These preliminary estimators can affect the asymptotic distribution of the PRTE estimator in complicated and intractable manners.
To solve this issue, we propose an orthogonal score for double debiased estimation of the PRTE.
Our proposed orthogonal score allows for the asymptotic distribution of the PRTE estimator to be obtained without any influence of preliminary parameter estimators as far as they satisfy mild convergence rate conditions.
Simulation results confirm our theoretical properties, and demonstrate that the method of estimation and inference works well even in small samples.
Our empirical application demonstrates that the proposed method indeed works with real data.
To our knowledge, our work is the first to develop asymptotic distribution theories for inference about the PRTE.
We hope that our proposed method contributes to empirical analyses of policy relevant treatment effects.
\clearpage