Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.
Identification and Debiased Learning of Causal Effects with General Instrumental Variables
\if11
\fi
\if01
{
center[center omitted — 122 chars of source]
} \fi
abstractInstrumental variable methods are fundamental to causal inference when treatment assignment is confounded by unobserved variables. In this article, we develop a general nonparametric causal framework for identification and learning with multi-categorical or continuous instrumental variables.
Specifically, the mean potential outcomes and the average treatment effect can be identified via a regular weighting function derived from the proposed framework.
Leveraging semiparametric theory, we derive efficient influence functions and construct two consistent, asymptotically normal estimators via debiased machine learning.
The first estimator uses a prespecified weighting function, while the second estimator selects the optimal weighting function adaptively.
Extensions to longitudinal data, dynamic treatment regimes, and multiplicative instrumental variables are further developed. We demonstrate the proposed method by employing simulation studies and analyzing real data from the Job Training Partnership Act program.
{\it Keywords:}
Causal inference,
Debiased machine learning,
Dynamic treatment regimes,
Instrumental variables,
Longitudinal data.
bibunit\section{Introduction}
\subsection{Background}
Observational studies are commonly employed to estimate treatment effects in biomedical and economic research hernan2020causal. In the presence of unmeasured confounding, instrumental variable (IV) methods have been widely used to identify causal effects Robins1994SNMM,imbens1994identification,angrist1996identification,abadie2003semiparametric. These approaches leverage exogenous variation in treatment induced by instruments that influence treatment assignment but are conditionally independent of latent confounding effects.
Under the monotonicity condition that the instrument does not decrease (or increase) the probability of receiving treatment for any individual, imbens1994identification,angrist1996identification show that one can identify causal effects within the complier subgroup with a binary IV.
Based on this work, heckman1999local,heckman2005structural,kennedy2019robust propose estimators for the local IV effect curve, capturing the treatment effect among individuals who would comply when the instrument exceeds a certain threshold. Recently, tchetgen2024nudge relaxes the monotonicity assumption, establishing the identification of the average treatment effect in the nudge subgroup, which involves mixtures of compliers and defiers.
Unlike approaches relying on the monotonicity condition, wang2018bounded, Hartwig2023 identify the average treatment effect (ATE) by imposing certain no additive interaction assumptions between the IV and the latent confounding in the treatment model; liu2025multiplicativeinstrumentalvariablemodel propose an identifying condition that excludes any multiplicative interaction between IV and latent confounders in the treatment model, thereby enabling the identification of the ATE of the treated group. Building upon a no-interaction assumption, cui2023instrumental,wang2023instrumental,michael2024instrumental propose marginal structural mean models and marginal structural Cox models, and fu2022offline,xu2023instrumental propose IV approaches to confounded off-policy evaluation using a similar identification strategy. Furthermore, ye2023instrumented propose an instrumental difference-in-difference method to identify the ATE based on no-interaction assumptions.
In addition, cui2021semiparametric propose to identify optimal treatment regimes under a no unmeasured common effect modifier assumption, and cui2025learning extend the approach to learn optimal treatment regimes for censored survival data.
In parallel, chernozhukov2005iv,chernozhukov2020instrumentalvariablequantileregression develop IV quantile regression for estimating quantile treatment effects, accommodating non-binary instruments.
chen2018sequential proposes sequential instrumental variables censored quantile regression estimators with a binary endogenous variable and a continuous instrument. More recently, chernozhukov2024estimatingcausaleffectsdiscrete introduce a copula invariance condition to identify treatment effects for the entire population with a binary instrumental variable.
Collectively, these approaches, however, do not offer nonparametric identification strategies for the ATEs and mean potential outcomes in settings where the instrument is non-binary.
Despite the nonparametric identification of ATEs or mean potential outcomes, several classical frameworks have been well developed to handle continuous IVs. The generalized method of moments provides a unifying and flexible parametric framework for estimation under IV settings hansen1982large. It estimates structural parameters of interest by solving the corresponding moment conditions. It encompasses traditional methods like two-stage least squares sargan1958estimation, wooldridge2010econometric as a special case and extends naturally to situations with multiple instruments, heteroskedasticity, and other complexities.
Nonparametric IV methods represent another important and widely adopted framework, particularly suitable for settings where the structural relationships between variables are complex and cannot be adequately captured by parametric models newey2003instrumental,darolles2011nonparametric,severini2012efficiency,AiChen2012,newey2013nonparametric,Chen_2018.
By avoiding restrictive assumptions on functional forms, these methods enable flexible, data-driven estimation of treatment and outcome models. In particular, proximal causal inference methods are closely related to nonparametric IV techniques and have emerged as effective approaches to address unmeasured confounding when identifying ATEs miao2018identifying,tchetgen2024introduction,cui2024semiparametric.
\subsection{Contributions}
We now outline the contents of the paper.
In Section (ref), we present the fundamental assumptions and introduce a novel concept termed regular weighting function, whose existence critically depends on the strength of the association between IV and treatment. Motivated by solving nonparametric IV problems, we further introduce a novel condition, termed additive IV, which provides a sufficient condition for the existence of nonparametric IV solutions. The additive IV condition can be viewed as a natural generalization of the no-interaction assumption. Moreover, we extend the assumption of the no unmeasured common effect modifier to the general IV setting, which also provides a sufficient condition for the existence of nonparametric IV solutions.
In Section (ref), we employ semiparametric theory vanderLaan2003, tsiatis2006semiparametric to derive the efficient influence functions (EIFs) for target causal estimands defined by various regular weighting functions. For identifying the ATE, we characterize the optimal weighting function, which yields an estimator achieving the semiparametric efficiency bound under the setting of homoskedastic latent confounding. Building upon the debiased machine learning framework schick1986asymptotically, Chernozhukov2018, we propose two novel estimators for the ATE.
The first estimator uses a prespecified weighting function, while the second estimator does not require a prespecified weighting function and selects the optimal weighting function adaptively.
In Section (ref), we extend the identification strategy from the point-exposure setting to longitudinal data.
We develop a cross-fitting procedure for estimation and establish that the proposed estimator is asymptotically normal, and its variance estimator is consistent.
Further extensions to dynamic treatment regimes and multiplicative IVs are provided in the Supplementary Material, demonstrating the versatility of our framework for identification and estimation.
Our contributions are fourfold.
First, we introduce a novel concept termed regular weighting function, which provides a new perspective on the IV relevance condition.
Second, we propose a novel identification strategy in the general IV setting by solving a special class of nonparametric IV problems, and we establish sufficient conditions for the existence of such solutions. Third, we derive EIFs and construct two consistent, asymptotically normal estimators via debiased machine learning. We note that the second estimator is non-standard, which does not require a prespecified weighting function and selects the optimal weighting function adaptively.
Fourth, we extend our identification framework to accommodate more complex settings, including longitudinal data structures, dynamic treatment regimes, and multiplicative IV models.
Simulation and empirical studies are conducted to demonstrate the validity and applicability of the proposed methods.
\section{Identification}
\subsection{Preliminaries}
Let \( A \in \mathcal{A} := \{0, \ldots, M\} \) denote a multi-categorical treatment variable, where \( M = 1 \) corresponds to the binary treatment setting.
Let \( Z \in \mathcal{Z}\subseteq\mathbb{R}^{|\mathcal{Z}|} \) denote the IV, which may be multi-categorical or continuous. Let \( U\in\mathcal{U} \subseteq\mathbb{R}^{|\mathcal{U}|}\) represent unmeasured confounders, \( L\in\mathcal{L} \subseteq\mathbb{R}^{|\mathcal{L}|}\) observable confounders, \( Y\in\mathcal{Y} \subseteq\mathbb{R}\) the observed outcome, and \( Y(a) \) the potential outcome at the treatment level \( A = a \).
The observed data consist of \( O = \{Z, A, Y, L\} \in \mathcal{Z}\times \mathcal{A}\times \mathcal{Y}\times\mathcal{L}\).
We introduce four fundamental assumptions in the IV setting.
\begin{assumption}[Consistency]
$Y=Y(A)$.
\end{assumption}
\begin{assumption}[Latent ignorability]
For any $a\in\mathcal{A}$, $Y(a)\perp\!\!\!\perp \{A,Z\}\mid U,L$.
\end{assumption}
\begin{assumption}[IV independence]
$Z\perp\!\!\!\perp U\mid L$.
\end{assumption}
\begin{assumption}[IV relevance]
For any $a\in\mathcal{A}$ and $l\in\mathcal{L}$, $Z\not\perp\!\!\!\perp I\{A=a\}\mid L=l$.
That is, there exist two distinct values $z_0,z_1\in\mathcal{Z}$ such that
$\Pr(A=a\mid Z=z_0,L=l)\neq \Pr(A=a\mid Z=z_1,L=l).$
\end{assumption}
Assumption (ref) posits that, conditional on both the observed covariates and unmeasured confounders, the potential outcome \( Y(a) \) is independent of the treatment and IV.
Assumption (ref) states that \( Z \) is independent of the unmeasured confounders \( U \) given the observed covariates \( L \).
Assumption (ref) requires that \( Z \) has a nontrivial effect on the treatment \( A \), conditional on any level of \( L \). This condition is slightly stronger than \( A \not\perp\!\!\!\perp Z \mid L \), which only requires the existence of some \( l \in \mathcal{L} \), \( a \in \mathcal{A} \), and $z_0, z_1 \in \mathcal{Z}$ such that $\Pr(A=a\mid Z=z_0,L=l)\neq \Pr(A=a\mid Z=z_1,L=l).$
\subsection{Regular weighting function}
Next, we introduce a novel concept of a regular weighting function, which is closely related to the IV relevance and the positivity condition.
\begin{definition}[Regular weighting function]
For each $a\in\mathcal{A}$, a function $\pi(Z,L)$ is a regular weighting function (RWF) for $A = a$ if it is uniformly bounded and there exists a positive constant $\epsilon_0$ such that
$|\mathrm{Cov}\!\{I\{A=a\}, \pi(Z,L) \mid L\} | \geq \epsilon_0$ uniformly for all $L$.
\end{definition}
The existence of an RWF for every \( a \in \mathcal{A} \) implies the IV relevance condition in Assumption (ref). Specifically, requiring the absolute value of the conditional covariance to be uniformly bounded below by a positive constant $\epsilon_0$ rules out scenarios where $\pi(Z,L)$ is irrelevant to $I\{A=a\}$, which would undermine the stability and validity of the identification of causal effects.
To guarantee the existence of an RWF, the following strong IV relevance assumption is required.
\begin{assumption}[Strong IV relevance]
There exists a positive constant $\epsilon_0$ such that for any $l\in\mathcal{L}$ and $a\in\mathcal{A}$,
$\mathrm{Var}\!\{\Pr(A=a\mid Z,L) \mid L=l\} \geq \epsilon_0.$
\end{assumption}
\begin{remark}
Here we directly impose a stronger version of the IV relevance condition in Assumption (ref), which is needed to derive the asymptotic properties of our proposed estimator. However, for identification or semiparametric analysis, only a weaker form is required, namely that for all \(a \in \mathcal{A}\),
$\mathrm{Var}\{\Pr(A=a \mid Z,L) \mid L=l\} \neq 0$.
A similar distinction also arises in our definition of RWFs.
\end{remark}
Assumption (ref) implies that the IV has a non-negligible effect on $A$, entailing the IV relevance condition in Assumption (ref). In fact, $\mathrm{Var}\!\{\Pr(A=a\mid Z,L) \mid L\}$ equals
\begin{align*}
\Pr(A=a\mid L)\times
\left\{\mathbb{E}[\Pr(A=a\mid Z,L)|A=a,L]-\mathbb{E}[\Pr(A=a\mid Z,L)|L]\right\}.
\end{align*}
Thus, Assumption (ref) also guarantees that for any $l \in \mathcal{L}$, $\Pr(A=a \mid L=l) \neq 0$, corresponding to the classical positivity condition hernan2020causal.
The following proposition further emphasizes the role of Assumption (ref) when identifying an RWF for $A=a$.
\begin{proposition}[Existence of RWFs]
There exists an RWF $\pi(Z,L)$ for $A$ if and only if Assumption (ref) holds. If there exists an RWF $\pi(Z,L)$ for $A=a$, then $\Pr(A=a\mid Z,L)$ must be an RWF for $A=a$.
\end{proposition}
In particular, if there exists an RWF $\pi(Z,L)$ for $A=a$, then there actually exists a whole family of RWFs, which can be generated by multiplying $\pi(Z,L)$ by any function $f(L)$ that is uniformly bounded above and away from zero.
In fact, according to Proposition (ref), Assumption (ref) implies that the function $\Pr(A=a \mid Z,L)$ itself is an RWF, thereby guaranteeing the existence of at least one RWF.
\subsection{Identification of the mean potential outcome}
In this subsection, we propose a strategy to identify the potential outcome mean \(\mathbb{E}[Y(a)]\) by formulating and solving a class of nonparametric IV problems. Throughout, we define \( A^{(a)} := I\{A = a\} \) for each \( a \in \mathcal{A} \) for convenience of notation. Our first primary goal is, for each \( a \in \mathcal{A} \), to identify a function \( f_a^o(A^{(a)}, L) \) that satisfies the conditional moment restriction
\begin{equation}
\mathbb{E}[A^{(a)} Y \mid Z, L] = \mathbb{E}[f_a^o(A^{(a)}, L) \mid Z, L].
\end{equation}
This conditional moment equation is common in the nonparametric IV literature newey2003instrumental.
The following theorem establishes the uniqueness of the solution to Equation (ref) and provides an explicit representation in terms of any RWF \(\pi(Z, L)\).
\begin{theorem}[Uniqueness and closed form solution]
Under Assumption (ref), for each $a\in\mathcal{A}$, if there is a solution \( f_a^o(A^{(a)}, L) \) to Equation (ref), it is unique. Moreover, for any RWF \( \pi(Z, L) \) for $A=a$, the solution \( f_a^o(A^{(a)}, L) \) satisfies
\begin{equation}
\begin{aligned}
f_a^o(0,L) &= -\dfrac{\mathrm{Cov}\!\{A^{(a)}Y, \pi(Z,L) \mid L\}}{\mathrm{Cov}\!\{A^{(a)}, \pi(Z,L) \mid L\}} \mathbb{E}[A^{(a)}\mid Z,L] + \mathbb{E}[A^{(a)}Y\mid Z,L],\\
f_a^o(1,L) &= \dfrac{\mathrm{Cov}\!\{A^{(a)}Y, \pi(Z,L) \mid L\}}{\mathrm{Cov}\!\{A^{(a)}, \pi(Z,L) \mid L\}} \{1-\mathbb{E}[A^{(a)}\mid Z,L]\} + \mathbb{E}[A^{(a)}Y\mid Z,L].
\end{aligned}
\end{equation}
\end{theorem}
However, Theorem (ref) only ensures uniqueness without guaranteeing the existence of a solution to Equation (ref). Indeed, Equation (ref) is over-identified when $Z$ is non-binary, meaning that a solution may not exist in general.
In the literature on proximal causal inference, the existence of a solution to the bridge equation is typically guaranteed under completeness conditions tchetgen2024introduction.
Likewise, within the IV framework, the additive IV condition plays a central role in ensuring the existence of a solution to the nonparametric IV problem in Equation (ref).
\begin{definition}[Additive IV]
For each $a \in \mathcal{A}$, $Z$ is an additive IV for $A = a$ if there exist functions $b(U,L)$ and $c(Z,L)$ such that
$\Pr(A = a \mid Z, U, L) = b(U, L) + c(Z, L).$
Furthermore, $Z$ is an additive IV for $A$ if it is an additive IV for $A=a$ for all $a\in\mathcal{A}$.
\end{definition}
According to wang2018bounded, the definition of additive IV is motivated by the no-interaction condition between $Z$ and $U$ in the treatment model.
tchetgen2021genius,sun2022selective,ye2024geniusmawiirobustmendelianrandomization also adopt analogous no-interaction conditions to identify the ATE with an invalid IV.
The following proposition gives an alternative definition for additive IV.
\begin{proposition}
For each $a \in \mathcal{A}$, $Z$ is an additive IV for $A = a$ if and only if for any $\pi(Z,L)$,
$\mathrm{Cov}\{I\{A=a\},\pi(Z,L)\mid U,L\}=\mathrm{Cov}\{I\{A=a\},\pi(Z,L)\mid L\}.$
\end{proposition}
Specifically, when \( Z \) is binary, the condition in Definition (ref) holds if and only if
$$\Pr(A=a \mid Z=1, U, L)- \Pr(A=a \mid Z=0, U, L) \perp\!\!\!\perp U \mid L,$$
implying that the differential effect of the instrument on treatment $A$ is conditionally independent of the unmeasured confounders \( U \), given the observed covariates \( L \).
Next, we establish the identification strategy of the mean potential outcome under an additive IV.
\begin{theorem}[Identification of the mean potential outcome]
Under Assumptions (ref)--(ref), for each \( a\in\mathcal{A} \), if $Z$ is an additive IV for $A=a$, then there exists a unique solution \( f_a^o(A^{(a)}, L) \) to Equation (ref), given by
\begin{equation*}
\begin{aligned}
f_a^o(0,L) &= \mathrm{Cov}\!\{\mathbb{E}[Y(a)\mid U,L], \Pr(A=a\mid Z,U,L)\mid L\},\\
f_a^o(1,L) &= \mathrm{Cov}\!\{\mathbb{E}[Y(a)\mid U,L], \Pr(A=a\mid Z,U,L)\mid L\} + \mathbb{E}[Y(a)\mid L].
\end{aligned}
\end{equation*}
Furthermore, if $\pi(Z,L)$ is an RWF for $A=a$,
\begin{align}
\mathbb{E}[Y(a)] = \mathbb{E}\left[\dfrac{\mathrm{Cov}\!\{A^{(a)}Y, \pi(Z,L) \mid L\}}{\mathrm{Cov}\!\{A^{(a)}, \pi(Z,L) \mid L\}}\right].
\end{align}
\end{theorem}
Notably, our identification strategy holds even when the treatment space $\mathcal{A}$ is multi-categorical. Intuitively, this is because we transform the multi-categorical treatment $A$ into a binary variable $A^{(a)}$ for each $a \in \mathcal{A}$.
Moreover, if there are no latent confounders $U$, then \( f_a^o(0, L) \equiv 0 \) by Theorem (ref). In other words, a significant deviation of \( f_a^o(0, L) \) from zero indicates the existence of latent confounding. This observation provides a novel criterion for detecting unmeasured confounder $U$, which is beyond the scope of our article.
Furthermore, when \( Z \) is a non-binary IV, Theorem (ref) provides a practical approach to assess whether \( Z \) qualifies as an additive IV for \( A = a \), which constitutes a relatively strong condition. Specifically, under Assumptions (ref)--(ref), if \( Z \) is indeed an additive IV for \( A = a \), the choice of the RWF \( \pi(Z, L) \) does not affect the expression value on the right-hand side of Equation (ref).
Therefore, two distinct RWFs \( \pi_1(Z, L) \) and \( \pi_2(Z, L) \) can be selected and the corresponding values can be compared. A discrepancy between these two values would indicate a violation of the additive IV condition.
\subsection{Identification of the average treatment effect}
In practical applications with binary treatment \( A \), researchers are also interested
in identifying the ATE. The following theorem summarizes the identification results for the ATE under a weaker assumption than additive IV.
\begin{theorem}[Identification of the ATE]
Under Assumptions (ref)--(ref), and
\begin{align}
\mathrm{Cov}\!\{Y(1)-Y(0),\mathbb{E}[A\mid Z=z,U,L]-\mathbb{E}[A\mid U,L]\mid L\}=0,
for any z\in\mathcal{Z},
\end{align}
the nonparametric IV problem
\begin{equation}
\mathbb{E}[Y\mid Z,L]=\mathbb{E}[f^o(A,L)\mid Z,L]
\end{equation}
has a unique solution \( f^o(A, L) \) as
\begin{align*}
f^o(0,L) &= \mathrm{Cov}\!\{\mathbb{E}[Y(1)-Y(0)\mid U,L], \Pr(A=1\mid U,L)\mid L\} + \mathbb{E}[Y(0)\mid L],\\
f^o(1,L) &= \mathrm{Cov}\!\{\mathbb{E}[Y(1)-Y(0)\mid U,L], \Pr(A=1\mid U,L)\mid L\} + \mathbb{E}[Y(1)\mid L].
\end{align*}
Furthermore, for any RWF \( \pi(Z, L) \) for $A$, the ATE is identified as
\begin{align}
\mathbb{E}[Y(1)-Y(0)] = \psi_{\pi}^o := \mathbb{E}\left[\dfrac{\mathrm{Cov}\!\{Y, \pi(Z,L) \mid L\}}{\mathrm{Cov}\!\{A, \pi(Z,L) \mid L\}}\right].
\end{align}
\end{theorem}
\begin{remark}
In fact, under Assumptions (ref)--(ref), if $\pi(Z,L)$ is an RWF for $A$,
\begin{align*}
\psi_{\pi}^o= \mathbb{E}\left[\mathbb{E}[Y(1)-Y(0)\mid U,L]\dfrac{\mathrm{Cov}\!\{A, \pi(Z,L) \mid U,L\}}{\mathrm{Cov}\!\{A, \pi(Z,L) \mid L\}}\right].
\end{align*}
This result suggests that, even when $Z$ does not satisfy Equation (ref), the estimand $\psi_{\pi}^o$ can still be interpreted as a weighted average of the conditional ATE $\mathbb{E}[Y(1)-Y(0)\mid U,L]$, as long as $\mathrm{Cov}\!\{A, \pi(Z,L) \mid U,L\}$ and $\mathrm{Cov}\!\{A, \pi(Z,L) \mid L\}$ have the same sign.
This representation is analogous to the assumption-lean inference of vansteelandt2022assumption,vansteelandt2024assumption, in which the estimands remain meaningful and interpretable even under model misspecification.
\end{remark}
\begin{remark}
For a fixed RWF $\pi(Z,L)$, Equation (ref) holds when
\begin{align}
\mathrm{Cov}\!\{Y(1)-Y(0),\,\mathrm{Cov}\!\{A, \pi(Z,L) \mid U,L\}\mid L\}=0.
\end{align}
One can check that Equation (ref) holds
if and only if Equation (ref) is satisfied for all $\pi(Z,L)$.
Both of them can be viewed as extensions of the “no unmeasured common effect modifier” assumption proposed by cui2021semiparametric.
\end{remark}
\begin{remark}
When \( L = \emptyset \) and $Z$ are univariate, the parameter $\psi_{\pi}^o$ in Equation (ref) reduces to a form resembling the two-stage least squares estimator if we set \(\pi(Z,L) = Z\). Moreover, several classical works interpret $\psi_{\pi}^o$ as $\tau_0$, the solution to the conditional moment equation
\(\mathbb{E}[Y - \tau_0 A \mid Z] = 0\)
hansen1982large,newey1994large, where \(\tau_0\) represents the effect of a marginal change in the endogenous variable \( A \) on the outcome. Consequently, our identification strategy can be viewed as a nonparametric generalization of the two-stage least squares and generalized method of moments, which additionally accounts for confounding effects through \( L \).
\end{remark}
\begin{remark}
For identifying the average dose-response curve of a continuous treatment, we provide a similar identification result in the Supplementary Material.
\end{remark}
\section{Semiparametric theory for learning the average treatment effect}
\subsection{Prespecified weighting function scenario}
Recall that under Assumptions (ref)--(ref), if Equation (ref) holds, the choice of the RWF \( \pi(Z, L) \) does not affect the value of \( \psi_{\pi}^o \) in Equation (ref), as shown in Theorem (ref).
However, the semiparametric efficiency bound of $\psi_{\pi}^o$ still depends on the choice of \( \pi(Z, L) \). Consequently, it is necessary to derive the EIFs corresponding to all possible RWFs \( \pi(Z, L) \).
To this end, we define the following nuisance functions:
\begin{align*}
&\delta^o(L) := \mathbb{E}[A \mid L],
&& \eta^o(L) := \mathbb{E}[Y \mid L],\\
&\kappa_{\pi}^o(L) := \mathrm{Cov}\!\{A, \pi(Z,L) \mid L\},
&& \zeta_{\pi}^o(L) := \mathbb{E}[Y \pi(Z,L) \mid L],\\
&\rho_{\pi}^o(L) := \mathbb{E}[\pi(Z,L) \mid L],
&& \gamma_{\pi}^o(L) := \dfrac{\mathrm{Cov}\!\{Y, \pi(Z,L) \mid L\}}{\mathrm{Cov}\!\{A, \pi(Z,L) \mid L\}}.
\end{align*}
For notational convenience, we collect them into a unified nuisance vector:
\begin{align}
\alpha_{\pi}^o(L) := [\delta^o(L), \kappa_{\pi}^o(L), \rho_{\pi}^o(L), \eta^o(L), \gamma_{\pi}^o(L)].
\end{align}
Note that \(\zeta_{\pi}^o(L)\) serves only as an intermediate nuisance function, which will be used in our proof and does not appear in the unified vector \(\alpha_{\pi}^o(L)\).
Next, we derive the EIF of \(\psi_{\pi}^o\) for any choice of RWF \(\pi(Z,L)\).
\begin{theorem}
Under Assumptions (ref)--(ref),
for any RWF $\pi(Z,L)$ for $A$, the EIF for $\psi_{\pi}^o$ in Equation (ref) is $\varphi_{\pi}(O;\psi_{\pi}^o,\alpha_{\pi}^o)$, where
$\varphi_{\pi}(O;\psi_{\pi},\alpha_{\pi})$ equals
\begin{align*}
\dfrac{ \left\{\pi(Z,L)-\rho_{\pi}(L)\right\}}{\kappa_{\pi}(L)} \{Y-\eta(L)\} - \psi_{\pi}
+ \left(1 - \dfrac{\{\pi(Z,L)-\rho_{\pi}(L)\}\{A-\delta(L)\}}{\kappa_{\pi}(L)} \right) \gamma_{\pi}(L).
\end{align*}
\end{theorem}
\begin{remark}
We carry out our semiparametric analysis in the fully nonparametric model,
and the unique influence function corresponds to the EIF.
\end{remark}
This theorem underpins the construction of efficient estimators for \( \psi_{\pi}^o \). In practice, the true nuisance vector \( \alpha_{\pi}^o \) is unknown and must be estimated from the data.
The following proposition quantifies the bias introduced by substituting an estimated \( \alpha_{\pi} \) for the true vector. This mixed bias property is well-documented in the existing literature rotnitzky2021characterization.
\begin{proposition}[Mixed bias property]
Under Assumptions (ref)--(ref),
for any RWF $\pi(Z,L)$ and any fixed nuisance function $\alpha_{\pi}$, $\varphi_{\pi}(O;\psi_{\pi}^o,\alpha_{\pi})$ satisfies that
\begin{align*}
&\mathbb{E}[\varphi_{\pi}(O;\psi_{\pi}^o,\alpha_{\pi})]
=\mathbb{E}\left[\dfrac{1}{\kappa_{\pi}(L)}\left\{\begin{array}{l}
\{\kappa_{\pi}(L)-\kappa_{\pi}^o(L)\}\{\gamma_{\pi}(L)-\gamma_{\pi}^o(L)\}\\
+\{\rho_{\pi}(L)-\rho_{\pi}^o(L)\}\{\eta(L)-\eta^o(L)\}\\
+\{\rho_{\pi}(L)-\rho_{\pi}^o(L)\}\{\delta(L)-\delta^o(L)\}\gamma_{\pi}(L)
\end{array}\right\}\right].
\end{align*}
\end{proposition}
Next, we aim to select a function \( \pi(Z, L) \) from the class of RWFs that minimizes the asymptotic variance of the estimator for \( \psi_{\pi}^o \). Specifically, we seek the optimal \( \pi(Z, L) \) that minimizes the second moment of the EIF.
Minimizing this quantity yields the most statistically efficient estimator among all estimators based on different RWFs.
Intuitively, Proposition (ref) suggests that
\(\pi_*^o(Z,L) := \Pr(A = 1 \mid Z, L)\) is a natural candidate for the optimal weighting function.
The following proposition characterizes the \(\pi(Z,L)\) that achieves this minimum variance bound.
\begin{proposition}[Optimal RWF]
Under the conditions of Theorem (ref), suppose that the solution $f^o(A,L)$ to the nonparametric IV problem in Equation (ref) satisfies
\begin{equation}
\mathbb{E}\left[\{Y-f^o(A,L)\}^2\middle| Z,L\right]\perp\!\!\!\perp Z\mid L.
\end{equation}
Then the quantity $\mathbb{E}[\varphi_{\pi}(O;\psi_{\pi}^o,\alpha_{\pi}^o)^2]$ reaches its lower bound when $\pi(Z,L)=\Pr(A=1\mid Z,L)$.
\end{proposition}
The condition in Equation (ref) implies that, after conditioning on $L$, the instrument $Z$ does not provide additional information about the residual variation in $Y$. wiemann2023optimal leverage a similar assumption to derive the “optimal instruments” in the special case where $L=\emptyset$.
\subsection{Adaptive weighting function scenario}
As established in Proposition (ref), choosing
\(\pi_*^o(Z,L) := \Pr(A = 1 \mid Z, L)\) achieves the optimal efficiency bound for estimating \(\psi_{\pi}^o\). Moreover, Proposition (ref) implies that if \(\pi_*^o(Z,L)\) does not constitute a valid regular weighting function for \(A\), then no valid alternative exists.
These two observations motivate our strategy of adaptively estimating the optimal weighting function.
Consequently, we take \(\pi_*^o(Z,L)\) as the regular weighting function for \(A\). Since \(\pi_*^o(Z,L)\) is unknown in practice, we treat it as a nuisance function and estimate it in a data-adaptive manner.
Concretely, we define the following nuisance functions
\begin{align*}
&\delta^o(L):=\mathbb{E}[A\mid L],
&&\eta^o(L):=\mathbb{E}[Y\mid L],\\
&\kappa^o(L):=\mathrm{Cov}\!\{A,\pi_*^o(Z,L)\mid L\},
&&\zeta^o(L):=\mathbb{E}[Y\pi_*^o(Z,L)\mid L],\\
&\xi^o(Z,L):=\mathbb{E}[Y\mid Z,L],
&&\gamma^o(L):=\dfrac{\mathrm{Cov}\!\{Y, \pi_*^o(Z,L) \mid L\}}{\mathrm{Cov}\!\{A, \pi_*^o(Z,L) \mid L\}}.
\end{align*}
We unify these nuisance functions into a nuisance vector:
\begin{equation}
\beta^o(Z,L):=[\pi_*^o(Z,L),\delta^o(L),\kappa^o(L),
\xi^o(Z,L),\eta^o(L),\gamma^o(L)].
\end{equation}
Note that $\zeta^o(L)$ is only an intermediate nuisance function and is not included in $\beta^o(Z,L)$.
According to Theorem (ref), the ATE can be identified as
\begin{equation}
\psi_{*}^o:=\mathbb{E}\left[\gamma^o(L)\right]
=\mathbb{E}\left[\dfrac{\pi_*^o(Z,L)-\delta^o(L)}{\kappa^o(L)}Y\right].
\end{equation}
We now proceed to derive the EIF for $\psi_{*}^o$.
\begin{theorem}
Under Assumptions (ref)-(ref) and (ref),
the EIF for $\psi_{*}^o$ is given by $\varphi(O;\psi_{*}^o,\beta^o)$, where
\begin{align*}
\varphi(O;\psi_{*},\beta):=&\dfrac{\pi_*(Z,L)-\delta(L)}{\kappa(L)}Y-\psi_{*}
+\dfrac{\gamma(L)}{\kappa(L)}\left\{
\kappa(L) + (A-\pi_*(Z,L))^2 - (A-\delta(L))^2
\right\}\\
&+\dfrac{1}{\kappa(L)}
\left\{\xi(Z,L)(A-\pi_*(Z,L))-\eta(L)(A-\delta(L))\right\}.
\end{align*}
\end{theorem}
\begin{remark}
Similarly to Theorem (ref),
this theorem still holds even if $Z$ does not satisfy Equation (ref) for $A$, because the definition of $\psi_{*}^o$ does not rely on the additive IV condition.
\end{remark}
Next, the EIF characterization in Theorem (ref) forms the foundation for analyzing the robustness of the proposed estimator, which is demonstrated in the next proposition.
\begin{proposition}[Mixed bias property]
Under Assumptions (ref)-(ref) and (ref),
for any fixed nuisance vector $\beta$ in Equation (ref), $\varphi(O;\psi_{*}^o,\beta)$ satisfies
\begin{align*}
\mathbb{E}[\varphi(O;\psi_{*}^o,\beta)]=
\mathbb{E}\left[\dfrac{1}{\kappa(L)}\left\{\begin{array}{l}
\{\gamma(L)-\gamma^o(L)\}\{\kappa(L) - \kappa^o(L)\}\\
- \gamma(L)(\delta^o(L)-\delta(L))^2\\
+\gamma(L)(\pi_*^o(Z,L)-\pi_*(Z,L))^2\\
-(\xi(Z,L)-\xi^o(Z,L))(\pi_*(Z,L)-\pi_*^o(Z,L))\\
+(\eta(L)-\eta^o(L))(\delta(L)-\delta^o(L))
\end{array}\right\}\right].
\end{align*}
\end{proposition}
\subsection{Cross-fitting procedure}
In this subsection, we adopt the cross-fitting procedure Chernozhukov2018 to construct debiased estimators for $\psi_{\pi}^o$ and $\psi_{*}^o$ in Equations (ref) and (ref). Without loss of generality, assume that the sample size \(n\) is evenly divisible by the number of folds \(K\).
We randomly partition the sample into \(K\) folds of equal size. Let \(I_k\) denote the set of indices belonging to the \(k\)-th fold, and let \(I_{-k}\) denote its complement. Denote by \(|I_k|\) the size of the fold \(I_k\). For any random variable \(O\), define the empirical average over fold \(k\) as
\(
\mathbb{E}_{nk}[O] := \sum_{i \in I_k} O_i/|I_k|.
\)
We further define the $L_2$ norm of the nuisance vector $\alpha_{\pi}(L)$ from Equation (ref) as
\[
\|\alpha_{\pi}(L)\|_2^2 := \|\delta(L)\|_2^2 + \|\eta(L)\|_2^2 + \|\kappa_{\pi}(L)\|_2^2 + \|\rho_{\pi}(L)\|_2^2 + \|\gamma_{\pi}(L)\|_2^2,
\]
and the $L_2$ norm of the nuisance vector $\beta(Z,L)$ from Equation (ref) as
\[
\|\beta(Z,L)\|_2^2 := \|\pi_*(Z,L)\|_2^2 + \|\xi(Z,L)\|_2^2 + \|\delta(L)\|_2^2 + \|\eta(L)\|_2^2 + \|\kappa(L)\|_2^2 + \|\gamma(L)\|_2^2.
\]
Next, for any fixed fold \( I_k\), the nuisance estimators \( \hat\alpha_{\pi}^{(n,k)} \) are trained using only the observations in \( I_{-k} \) with any suitable machine learning methods. By construction, \( \hat\alpha_{\pi}^{(n,k)} \) is independent of the samples in $I_k$. We derive the estimator $\hat\psi_{\pi}^{(n)}$ as the solution to the equation
\begin{align}
\sum_{k=1}^K\mathbb{E}_{nk}\left[\varphi_{\pi}(O;\hat\psi_{\pi}^{(n)},\hat\alpha_{\pi}^{(n,k)})\right]=0.
\end{align}
Next, we establish the consistency and asymptotic normality of \(\hat\psi_{\pi}^{(n)}\) defined in Equation (ref).
\begin{theorem}[Asymptotic normality of $\hat\psi_{\pi}^{(n)}$]
Under Assumptions (ref)--(ref), suppose that
\(\pi(Z,L)\) is an RWF for $A$. Assume further that, for any \(k=1,\ldots,K\),
\(\mathbb{E}[\|\hat \alpha_{\pi}^{(n,k)}-\alpha_{\pi}^o\|_2^2]=o(1)\), and that
\begin{align*}
\left\{\begin{array}{l}
\|\hat\kappa_{\pi}^{(n,k)}-\kappa_{\pi}^o\|_2\times \|\hat\gamma_{\pi}^{(n,k)}-\gamma_{\pi}^o\|_2\\
+\|\hat\rho_{\pi}^{(n,k)}-\rho_{\pi}^o\|_2\times
\|\hat\delta^{(n,k)}-\delta^o\|_2\\
+\|\hat\rho_{\pi}^{(n,k)}-\rho_{\pi}^o\|_2\times
\|\hat\eta^{(n,k)}-\eta^o\|_2\\
\end{array}\right\}
=o_p(n^{-1/2}).
\end{align*}
Then $\sqrt{n}\left(\hat\psi_{\pi}^{(n)}-\psi_{\pi}^o\right)/\sigma_{\pi}^o$ converges in distribution to $\mathcal{N}(0,1)$, where the asymptotic variance is defined as $(\sigma_{\pi}^o)^2:=\mathbb{E}[\varphi_{\pi}(O;\psi_{\pi}^o,\alpha_{\pi}^o)^2]$. In addition, if we define
$$(\hat\sigma_{\pi}^{(n)})^2:=\sum_{k=1}^K\mathbb{E}_{nk}[\varphi_{\pi}(O;\hat\psi_{\pi}^{(n)},\hat\alpha_{\pi}^{(n,k)})^2]/K,$$
then $(\hat\sigma_{\pi}^{(n)})^2$ converges in probability to $(\sigma_{\pi}^o)^2$.
\end{theorem}
Analogously, for any fixed fold \(k = 1, \ldots, K\), nuisance estimators \(\hat\beta^{(n,k)}\) are trained using only the observations in \(I_{-k}\) with any suitable machine learning method. The estimator \(\hat\psi_{*}^{(n)}\) is then defined as the solution to
\begin{align}
\sum_{k=1}^K \mathbb{E}_{nk}\!\left[\varphi\!\left(O;\hat\psi_{*}^{(n)},\hat\beta^{(n,k)}\right)\right] = 0.
\end{align}
Next, we establish the consistency and asymptotic normality of \(\hat\psi_{*}^{(n)}\) defined in Equation (ref).
\begin{theorem}[Asymptotic normality of $\hat\psi_{*}^{(n)}$]
Under Assumptions (ref)-(ref) and (ref), suppose that Equation (ref) holds. Assume that for any $k=1,\ldots,K$, $\mathbb{E}[\|\hat \beta^{(n,k)}-\beta^o\|_2^2]=o(1)$, and that
\begin{align*}
\left\{
\begin{array}{l}
\|\hat\gamma^{(n,k)}-\gamma^o\|_2 \times \|\hat\kappa^{(n,k)}-\kappa^o\|_2
+\|\hat\delta^{(n,k)}-\delta^o\|_2^2+\|\hat\pi_*^{(n,k)}-\pi_*^o\|_2^2\\
+\|\hat\xi^{(n,k)}-\xi^o\|_2\times \|\hat\pi_*^{(n,k)}-\pi_*^o\|_2
+\|\hat\eta^{(n,k)}-\eta^o\|_2\times \|\hat\delta^{(n,k)}-\delta^o\|_2
\end{array}\right\}
=o_p(n^{-1/2})
\end{align*}
Then $\sqrt{n}\left(\hat\psi_{*}^{(n)}-\psi_{*}^o\right)/\sigma_{*}^o$ converges in distribution to $\mathcal{N}(0,1)$, where the asymptotic variance is defined as $(\sigma_{*}^o)^2:=\mathbb{E}[\varphi(O;\psi_{*}^o,\beta^o)^2].$ In addition, if we define
$$(\hat\sigma_{*}^{(n)})^2:=\sum_{k=1}^K\mathbb{E}_{nk}[\varphi(O;\hat\psi_{*}^{(n)},\hat\beta^{(n,k)})^2]/K,$$
then $(\hat\sigma_{*}^{(n)})^2$ converges in probability to $(\sigma_{*}^o)^2$.
\end{theorem}
\section{Extension to longitudinal data}
\subsection{Identification strategy}
In this section, we extend the identification strategy introduced in Section (ref) to a longitudinal setting with a sequence of IVs observed at each time point. Specifically, consider a longitudinal study with measurements collected at discrete time points \(T+1\), indexed by \(t = 0, 1, \ldots, T\), where \(T\) is a fixed nonnegative integer. The special case \(T = 0\) reduces to the panel data framework discussed in Section (ref).
For notation, let \(\overline{a}_t := [a_0, a_1, \ldots, a_t]\), \(\underline{a}_t := [a_t, a_{t+1}, \ldots, a_T]\), \(a_t^s := [a_t, a_{t+1}, \ldots, a_s]\), and \(\overline{a} := [a_0, a_1, \ldots, a_T]\). By convention, we set $\underline{a}_{T+1}$ and $a_t^{t-1}$ as empty for any \(t\), and note that \(\overline{a} = \underline{a}_0 = \overline{a}_T\).
At each time point \(t\), let \(L_t\in\mathcal{L}_t\) denote the vector of observed confounders, \(U_t\in\mathcal{U}_t\) the vector of unobserved confounders, \(Z_t\in\mathcal{Z}_t\) the IV (which may be multi-categorical or continuous), and \(A_t\in\mathcal{A}_t\) the discrete treatment assignment. The observed data are given by
\(
O := [\overline{Z}_T, \overline{A}_T, \overline{L}_T, Y],
\)
where \(Y\) is the outcome of interest, observed only at time \(T+1\).
At each time point \(t=0,\ldots, T\), define the history \(H_t := [\overline{Z}_{t-1}, \overline{A}_{t-1}, \overline{L}_t]\in\mathcal{H}_t\), and let \(H_{T+1} := O\) denote all the observed data. Notably, the histories satisfy the recursive relation
\(
H_t = [H_{t-1}, Z_{t-1}, A_{t-1}, L_t].
\)
Let \(Y(\overline{a})\) denote the potential outcome in the treatment history \(\overline{A} = \overline{a}\). Next, we extend the preceding assumptions to a longitudinal setting with IVs.
\begin{assumption}[Consistency]
\(Y=Y(\overline{A})\).
\end{assumption}
\begin{assumption}[Latent ignorability]
For any fixed $t$, $\{Z_t,A_t\}\perp\!\!\!\perp Y(\underline{a}_t) \mid H_t,\overline{U}_t$.
\end{assumption}
\begin{assumption}[IV independence]
For any fixed $t$, $Z_t\perp\!\!\!\perp \overline{U}_t\mid H_t$.
\end{assumption}
\begin{assumption}[Strong IV relevance]
For each time $t=0,\ldots,T$, there exists a positive constant $\epsilon_0$ such that for any $h_t\in\mathcal{H}_t$ and $a_t\in\mathcal{A}_t$,
$\mathrm{Var}\!\{\Pr(A_t=a_t\mid Z_t,H_t) \mid H_t=h_t\} \geq \epsilon_0.$
\end{assumption}
Assumptions (ref)--(ref) can be regarded as a longitudinal extension of Assumptions (ref)-(ref) and (ref).
For illustration, Figure (ref) displays a sequential directed acyclic graph (DAG) for the IV setting with \(T = 2\) under the one-step Markov property, where Assumptions (ref), (ref) hold.
\begin{figure}
\resizebox{0.4\textwidth}{!}{
\begin{tikzpicture}[
->,
>=Stealth,
node distance=1.5cm,
thick,
every node/.style={
draw=black,
rounded corners=3pt,
minimum size=0.5cm,
inner sep=2pt,
font=,
fill=white
},
every edge/.style={
draw=black,
->,
shorten >=1pt,
shorten <=1pt
}
]
\node (Z0) at (0,-1.5) {$Z_0$};
\node (Z1) at (2.25,-1.5) {$Z_1$};
\node (Z2) at (4.5,-1.5) {$Z_2$};
\node (A0) at (0,0) {$A_0$};
\node (A1) at (2.25,0) {$A_1$};
\node (A2) at (4.5,0) {$A_2$};
\node (U0)[fill=gray!20] at (0,1.5) {$U_0$};
\node (U1)[fill=gray!20] at (2.25,1.5) {$U_1$};
\node (U2)[fill=gray!20] at (4.5,1.5) {$U_2$};
\node (L0) at (0,-3) {$L_0$};
\node (L1) at (2.25,-3) {$L_1$};
\node (L2) at (4.5,-3) {$L_2$};
\node (Y) at (6,0) {$Y$};
\draw (Z0) -> (A0);
\draw (Z0) -> (Z1);
\draw (A0) -> (A1);
\draw (A0) -> (Z1);
\draw (A0) -> (U1);
\draw (A0) -> (L1);
\draw (U0) -> (A0);
\draw (U0) -> (U1);
\draw (U0) -> (L1);
\draw (Z1) -> (A1);
\draw (Z1) -> (Z2);
\draw (A1) -> (A2);
\draw (A1) -> (Z2);
\draw (A1) -> (U2);
\draw (A1) -> (L2);
\draw (U1) -> (A1);
\draw (U1) -> (U2);
\draw (U1) -> (L2);
\draw (Z2) -> (A2);
\draw (U2) -> (A2);
\draw (L0) to [bend left=30] (A0);
\draw (L0) -> (Z0);
\draw (L0) -> (L1);
\draw (L0) to [bend left=30] (U0);
\draw (L1) to [bend left=30] (A1);
\draw (L1) -> (Z1);
\draw (L1) -> (L2);
\draw (L1) to [bend left=30] (U1);
\draw (L2) to [bend left=30] (A2);
\draw (L2) -> (Z2);
\draw (L2) -> (Y);
\draw (L2) to [bend left=30] (U2);
\draw (U2) -> (Y);
\draw (A2) -> (Y);
\end{tikzpicture}}
\caption{
Sequential DAG when $T=2$.
Gray nodes indicate unobserved variables, while white nodes indicate observed ones.
}
\end{figure}
Next, we generalize the definitions of RWF and additive IV from the point-exposure setting to accommodate longitudinal data.
\begin{definition}[Longitudinal RWF]
A function \(\pi_t(Z_t, H_t)\) is said to be an RWF for \(A_t\) if it is uniformly bounded, and for each \(a_t \in \mathcal{A}_t\), there exists a constant \(\epsilon_0 > 0\) such that
\[
\left| \mathrm{Cov}\!\left\{ I\{A_t = a_t\}, \pi_t(Z_t, H_t) \mid H_t =h_t\right\} \right|
\geq \epsilon_0,
\quad \text{for all } h_t\in\mathcal{H}_t.
\]
\end{definition}
\begin{definition}[Longitudinal additive IV]
For any fixed \(t\), we say that \(Z_t\) is an additive IV for \(A_t\) if there exist functions \(b_{t,a_t}(\overline{U}_t,H_t)\) and \(c_{t,a_t}(Z_t,H_t)\) such that, for all \(a_t \in \mathcal{A}_t\),
\[
\Pr(A_t = a_t \mid Z_t, \overline{U}_t, H_t)
= b_{t,a_t}(\overline{U}_t, H_t) + c_{t,a_t}(Z_t, H_t).
\]
\end{definition}
These two definitions are required for identifying the longitudinal mean potential outcomes. For a fixed sequence of RWFs \(\pi_t(Z_t,H_t)\), define \(\gamma_{T+1,\underline{a}_{T+1}}^o(H_{T+1}) := Y.\)
Then, for \(t = T, \ldots, 0\), recursively define the nuisance function
\[
\gamma_{t,\underline{a}_t}^o(H_t)
:= \frac{\mathrm{Cov}\!\left\{ I\{A_t = a_t\} \, \gamma_{t+1,\underline{a}_{t+1}}^o(H_{t+1}), \pi_t(Z_t, H_t) \mid H_t \right\}}
{\mathrm{Cov}\!\left\{ I\{A_t = a_t\}, \pi_t(Z_t, H_t) \mid H_t \right\}}.
\]
Intuitively, \(\gamma_{t,\underline{a}_t}^o(H_t)\) can be interpreted as an estimator of the conditional mean potential outcomes \(\mathbb{E}[Y(\underline{a}_t) \mid H_t]\).
Building on this intuition, we derive the following identification result for longitudinal additive IVs.
\begin{theorem}[Longitudinal additive IV identification]
Under Assumptions (ref)--(ref), let $s$, $r$ be two integers with \(0 \le s \le T+1\), \(r \ge 0\), and \(s + r \le T+1\). Suppose that for each \(t=0,\ldots,T\), \(Z_t\) serves as an additive IV for \(A_t\), and that \(\pi_t(Z_t, H_t)\) is an RWF for \(A_t\). Then, the mean potential outcomes \(\mathbb{E}\bigl[Y(\underline{a}_{s})\bigr]\) can be expressed as
\begin{equation}
\begin{aligned}
\mathbb{E}\left[
\prod_{t=s}^{T-r}
\dfrac{\left(\pi_t(Z_t,H_t)-\mathbb{E}[\pi_t(Z_t,H_t)\mid H_t]\right)I\{A_t=a_t\}}{\mathrm{Cov}\!\{I\{A_t=a_t\},\pi_t(Z_t,H_t)\mid H_t\}}
\gamma_{T+1-r,a_{T+1-r}}^o(H_{T+1-r})
\right].
\end{aligned}
\end{equation}
\end{theorem}
In particular, consider the special case of Theorem (ref) with \(s = 0\) and \(r = 0\), which corresponds to identifying the mean of potential outcomes from the initial time point to the final time \(T\). Without loss of generality, if we set \(\pi(Z_t,H_t) = Z_t\), the identification formula simplifies to
\begin{align}
\mathbb{E}[Y(\overline{a})] = \psi_{\overline{a}}^o :=
\mathbb{E}\Biggl[
\prod_{t=0}^{T}
\frac{(Z_t - \mathbb{E}[Z_t \mid H_t]) \, I\{A_t = a_t\}}{\mathrm{Cov}\!\{ I\{A_t = a_t\}, Z_t \mid H_t \}}
\times Y
\Biggr].
\end{align}
This expression is directly analogous to the inverse probability weighting approach in settings without unmeasured confounding.
Alternatively, setting \(s = 0\) and \(r = T+1\), the identification formula reduces to
\(
\psi_{\overline{a}}^o = \mathbb{E}\bigl[\gamma_{0,\overline{a}}^o(H_0)\bigr],
\)
which corresponds to the outcome regression approach or the g-formula in longitudinal causal inference hernan2020causal.
\subsection{Semiparametric theory}
In this subsection, without loss of generality, we focus on the case where $Z_t$ is univariate and \(\pi_t(Z_t, H_t)=Z_t\) is an RWF for \(A_t\). We derive the EIFs for \(\psi_{\overline{a}}^o\) in Equation (ref) when \(\pi_t(Z_t, H_t) = Z_t\); for a general RWF \(\pi_t(Z_t, H_t)\), we can define \(Z_t^{\pi} := \pi_t(Z_t, H_t)\) and replace \(Z_t\) with \(Z_t^{\pi}\) in the subsequent analysis.
For notational convenience, define \(\gamma_{T+1,\underline{a}_{T+1}}^o(H_{T+1}) := Y\) and \(A_t^{(a_t)} := I\{A_t = a_t\}\).
Next, for \(t = T, \ldots, 0\), recursively define nuisance functions:
\begin{align*}
&\kappa_{t,a_t}^o(H_t) := \mathrm{Cov}\!\{A_t^{(a_t)}, Z_t \mid H_t\},\\
&\delta_{t,a_t}^o(H_t) := \mathbb{E}[A_t^{(a_t)} \mid H_t],
&& \eta_{t,a_t}^o(H_t) := \mathbb{E}\bigl[A_t^{(a_t)} \gamma_{t+1,a_{t+1}}^o(H_{t+1}) \mid H_t\bigr],\\
&\rho_t^o(H_t) := \mathbb{E}[Z_t \mid H_t],
&& \gamma_{t,\underline{a}_t}^o(H_t) := \frac{\mathrm{Cov}\!\{ A_t^{(a_t)} \gamma_{t+1,\underline{a}_{t+1}}^o(H_{t+1}), Z_t \mid H_t \}}
{\mathrm{Cov}\!\{ A_t^{(a_t)}, Z_t \mid H_t \}}.
\end{align*}
We unify these nuisance functions into one nuisance vector:
\begin{align}
\alpha_{\overline{a}}^o := \{\alpha_{t,\underline{a}_t}^o\}_{t=0}^T, \qquad
\alpha_{t,\underline{a}_t}^o := [
\delta_{t,a_t}^o,
\kappa_{t,a_t}^o,
\rho_t^o,
\eta_{t,\underline{a}_t}^o,
\gamma_{t,\underline{a}_t}^o].
\end{align}
We now proceed to derive the EIF for \(\psi_{\overline{a}}^o\) in Equation (ref) in the following theorem.
\begin{theorem}
Under Assumptions (ref)--(ref), suppose that for each \(t=0,\ldots,T\),
\(\pi_t(Z_t, H_t)=Z_t\) is an RWF for \(A_t=a_t\). Then, the EIF for $\psi_{\overline{a}}^o$ is $\varphi_{\overline{a}}(O;\psi_{\overline{a}}^o,\alpha_{\overline{a}}^o)$, where
\begin{align*}
&\varphi_{\overline{a}}(O;\psi_{\overline{a}},\alpha_{\overline{a}}):=
\prod_{t=0}^{T}
\dfrac{\left\{Z_t-\rho_t(H_t)\right\}A_t^{(a_t)}}
{\kappa_{t,a_t}(H_t)}
Y-\psi_{\overline{a}}
+\sum_{t=0}^T\left(\prod_{s=0}^{t-1}\dfrac{\left\{Z_s-\rho_s(H_s)\right\}A_s^{(a_s)}}
{\kappa_{s,a_s}(H_s)}\right)
\&\times \left\{
\left(1-\dfrac{\{Z_t-\rho_t(H_t)\}\{A_t^{(a_t)}-\delta_{t,a_t}(H_t)\}}{\kappa_{t,a_t}(H_t)} \right)\gamma_{t,\underline{a}_t}(H_t) -
\dfrac{(Z_t-\rho_t(H_t))\eta_{t,\underline{a}_t}(H_t)}{\kappa_{t,a_t}(H_t)}\right\}.
\end{align*}
\end{theorem}
Notably, when $t=0$, the EIF in Theorem (ref)
has a structure similar to that derived in Theorem (ref).
The next proposition derives the mixed bias property for the EIF in Theorem (ref).
\begin{proposition}[Mixed bias property]
Under the conditions of Theorem (ref), for any fixed nuisance vector $\alpha_{\overline{a}}$ in Equation (ref), $\mathbb{E}[\varphi_{\overline{a}}(O;\psi_{\overline{a}}^o,\alpha_{\overline{a}})]$ equals
\begin{align*}
&\mathbb{E}\left[\begin{array}{l}
\sum_{t=0}^T\left(\prod_{s=0}^{t-1}\dfrac{\left\{Z_s-\rho_s(H_s)\right\}A_s^{(a_s)}}
{\kappa_{s,a_s}(H_s)}\right)
\times \dfrac{1}{\kappa_{t,a_t}(H_t)}\\
\times\left(\begin{array}{l}
\left\{\kappa_{t,a_t}(H_t)- \kappa_{t,a_t}^o(H_t)\right\} \left\{\gamma_{t,\underline{a}_t}(H_t)-\gamma_{t,\underline{a}_t}^o(H_t)\right\}\\
+\left\{\rho_t(H_t) - \rho_t^o(H_t)\right\}\left\{\eta_{t,\underline{a}_t}(H_t)-\eta_{t,\underline{a}_t}^o(H_t)\right\}\\
+\left\{\rho_t(H_t) - \rho_t^o(H_t)\right\}\left\{\delta_{t,a_t}(H_t)-\delta_{t,a_t}^o(H_t)\right\} \gamma_{t,\underline{a}_t}(H_t)
\end{array}\right)
\end{array}\right].
\end{align*}
\end{proposition}
\subsection{Cross-fitting procedure}
We develop a cross-fitting procedure for estimating $\psi_{\overline{a}}^o$. This estimator can be intuitively understood as resulting from a backward fitting strategy. Let
\(\hat{\mathbb{E}}^{(n,k)}[\phi(O) \mid H_t]\) denote an estimate of the conditional expectation \(\mathbb{E}[\phi(O)\mid H_t]\), and
\(\widehat{\text{Cov}}^{(n,k)}\{\phi_1(O), \phi_2(O)\mid H_t\}\) denote an estimate of the conditional covariance
\(\text{Cov}\{\phi_1(O),\phi_2(O)\mid H_t\}\). Both estimates are obtained using only the observations in \(I_{-k}\) and fitted using an appropriate machine learning method. The Algorithm (ref) summarizes the backward cross-fitting procedure.
\begin{algorithm}[htbp]
\caption{Backward Cross-Fitting Procedure for Longitudinal Data}
Randomly divide the samples evenly into \(K\) folds \(\{I_k\}_{k=1}^K\)\;
\For{$k = 1$ \KwTo $K$}{
Initialize \(t \gets T\) and set \(\hat\Psi_{T+1,\underline{a}_{T+1}}^{(n,k)} := Y\)\;
\While{$t \ge 0$}{
Fit the nuisance functions \(\hat \alpha_{t,\overline{a}}^{(n,k)}(H_t)\) using samples in \(I_{-k}\):
\begin{align*}
&\hat\eta_{t,\underline{a}_t}^{(n,k)}(H_t) := \hat{\mathbb{E}}^{(n,k)}[A_t^{(a_t)} \hat\Psi_{t+1,\underline{a}_{t+1}}^{(n,k)} \mid H_t], &
\hat \delta_{t,a_t}^{(n,k)}(H_t) := \hat{\mathbb{E}}^{(n,k)}[A_t^{(a_t)} \mid H_t],\\
&\hat\gamma_{t,\underline{a}_t}^{(n,k)}(H_t) :=
\frac{\widehat{\mathrm{Cov}}^{(n,k)}\{ A_t^{(a_t)} \hat\Psi_{t+1,\underline{a}_{t+1}}^{(n,k)}, Z_t \mid H_t\}}
{\widehat{\mathrm{Cov}}^{(n,k)}\{ A_t^{(a_t)}, Z_t \mid H_t \}}, &
\hat \rho_t^{(n,k)}(H_t) := \hat{\mathbb{E}}^{(n,k)}[Z_t \mid H_t],\\
&\hat\kappa_{t,a_t}^{(n,k)}(H_t) := \widehat{\mathrm{Cov}}^{(n,k)}\{A_t^{(a_t)},Z_t \mid H_t\}.
\end{align*}
Update for all samples:
\begin{align*}
\hat\Psi_{t,\underline{a}_{t}}^{(n,k)} :=&
\frac{1}{\hat\kappa_{t,a_t}^{(n,k)}(H_t)}
\left(Z_t - \hat\rho_t^{(n,k)}(H_t)\right)
\left(A_t^{(a_t)} \hat\Psi_{t+1,\underline{a}_{t+1}}^{(n,k)} - \hat\eta_{t,\underline{a}_t}^{(n,k)}(H_t)\right)
+ \&\left(1 -
\frac{(Z_t-\hat\rho_t^{(n,k)}(H_t))(A_t^{(a_t)} - \hat\delta_{t,a_t}^{(n,k)}(H_t))}{\hat\kappa_{t,a_t}^{(n,k)}(H_t)}\right)
\hat\gamma_{t,\underline{a}_t}^{(n,k)}(H_t).
\end{align*}
Decrement \(t \gets t-1\)\;
}
}
Output the estimator and its variance:
\begin{equation}
\hat\psi_{\overline{a}}^{(n)} := \frac{1}{K} \sum_{k=1}^K \mathbb{E}_{nk}[\hat\Psi_{0,\overline{a}}^{(n,k)}],\quad
(\hat\sigma_{\overline{a}}^{(n)})^2 := \frac{1}{K} \sum_{k=1}^K \mathbb{E}_{nk} \Big[(\hat\Psi_{0,\overline{a}}^{(n,k)} - \mathbb{E}_{nk}[\hat\Psi_{0,\overline{a}}^{(n,k)}])^2 \Big].
\end{equation}
\end{algorithm}
\begin{figure}[htbp]
\resizebox{0.8\textwidth}{!}{
\begin{tikzpicture}[
arrow/.style={-{Stealth[length=2mm, width=1.2mm]}, thick, gray!70},
every node/.style={minimum width=1.6cm, minimum height=0.9cm, align=center, font=},
state/.style={
draw=gray!70,
rounded corners=3pt,
top color=white,
bottom color=blue!5,
very thick,
drop shadow
},
est/.style={
circle,
draw=orange!60!black,
top color=white,
bottom color=orange!20,
thick,
drop shadow
}
]
\node[state] (T3) at (0, 0) {$\hat\Psi_{3,\underline{a}_3}^{(n,k)},H_3$};
\node[state] (T2) at (3, 0) {$\hat\Psi_{2,\underline{a}_2}^{(n,k)},H_2$};
\node[state] (T1) at (2*3, 0) {$\hat\Psi_{1,\underline{a}_1}^{(n,k)},H_1$};
\node[state] (T0) at (3*3, 0) {$\hat\Psi_{0,\underline{a}_0}^{(n,k)},H_0$};
\node[est] (alpha2) at (3, -1.8) {$\hat\alpha_{2,\overline{a}}^{(n,k)}$};
\node[est] (alpha1) at (2*3, -1.8) {$\hat\alpha_{1,\overline{a}}^{(n,k)}$};
\node[est] (alpha0) at (3*3, -1.8) {$\hat\alpha_{0,\overline{a}}^{(n,k)}$};
\node[state] (T3') at (0, -2*1.8) {$\hat\Psi_{3,\underline{a}_3}^{(n,k)},H_3$};
\node[state] (T2') at (3, -2*1.8) {$\hat\Psi_{2,\underline{a}_2}^{(n,k)},H_2$};
\node[state] (T1') at (2*3, -2*1.8) {$\hat\Psi_{1,\underline{a}_1}^{(n,k)},H_1$};
\node[state] (T0') at (3*3, -2*1.8) {$\hat\Psi_{0,\underline{a}_0}^{(n,k)},H_0$};
\node[draw=none] at (0, 0.9) {\scriptsize \textbf{$t=3$}};
\node[draw=none] at (3, 0.9) {\scriptsize \textbf{$t=2$}};
\node[draw=none] at (2*3, 0.9) {\scriptsize \textbf{$t=1$}};
\node[draw=none] at (3*3, 0.9) {\scriptsize \textbf{$t=0$}};
\node[draw=none, anchor=east, text=blue!50!black] at (-1, 0) {\textbf{ Evaluation set} $I_k$};
\node[draw=none, anchor=east, text=orange!60!black] at (-1, -2*1.8) {\textbf{ Training set} $I_{-k}$};
\draw[arrow] (T3) -- (T2);
\draw[arrow] (T2) -- (T1);
\draw[arrow] (T1) -- (T0);
\draw[arrow] (T3') -- (T2');
\draw[arrow] (T2') -- (T1');
\draw[arrow] (T1') -- (T0');
\draw[arrow, orange!70!black] (T3') -- (alpha2);
\draw[arrow, orange!70!black] (T2') -- (alpha1);
\draw[arrow, orange!70!black] (T1') -- (alpha0);
\draw[arrow, dashed, blue!60!black] (alpha2) -- (T2');
\draw[arrow, dashed, blue!60!black] (alpha2) -- (T2);
\draw[arrow, dashed, blue!60!black] (alpha1) -- (T1');
\draw[arrow, dashed, blue!60!black] (alpha1) -- (T1);
\draw[arrow, dashed, blue!60!black] (alpha0) -- (T0');
\draw[arrow, dashed, blue!60!black] (alpha0) -- (T0);
\end{tikzpicture}
}
\caption{Illustration of the backward cross-fitting procedure in the $k$-th fold for $T=2$.}
\end{figure}
Figure (ref) provides a graphical illustration of
Algorithm (ref) for \(T=2\). The upper row corresponds to the
evaluation set \(I_k\), while the bottom row represents the training set
\(I_{-k}\). The middle row displays the nuisance estimators
\(\hat\alpha_{t,\overline{a}}^{(n,k)}\), which are fitted using the training data and
adopted to form the final predictions. Importantly,
\(\hat\alpha_{t,\overline{a}}^{(n,k)}\) is constructed exclusively from samples in
\(I_{-k}\), thereby ensuring its independence from the evaluation set \(I_k\).
Moreover, in Algorithm (ref), we introduce the random variable $\hat\Psi_{t,\underline{a}_{t}}^{(n,k)}$. By induction, one can verify that for any $t=0,\ldots,T$,
\begin{align*}
&\hat\Psi_{t,\underline{a}_{t}}^{(n,k)}=\prod_{s=t}^{T} \dfrac{\left\{Z_s-\hat\rho_s^{(n,k)}(H_s)\right\}A_s^{(a_s)}} {\hat\kappa_{s,a_s}^{(n,k)}(H_s)} Y +\sum_{s=t}^T\left(\prod_{r=t}^{s-1}\dfrac{\left\{Z_r-\hat\rho_r^{(n,k)}(H_r)\right\}A_r^{(a_r)}} {\hat\kappa_{r,a_r}^{(n,k)}(H_r)}\right) \&\times \left\{\begin{array}{c}
\left(1-\dfrac{\{Z_s-\hat\rho_s^{(n,k)}(H_s)\}\{A_s^{(a_s)}-\hat\delta_{s,a_s}^{(n,k)}(H_s)\}}{\hat\kappa_{s,a_s}^{(n,k)}(H_s)} \right)\hat\gamma_{s,\underline{a}_s}^{(n,k)}(H_s) \\- \dfrac{(Z_s-\hat\rho_s^{(n,k)}(H_s))\times\hat\eta_{s,\underline{a}_s}^{(n,k)}(H_s)}{\hat\kappa_{s,a_s}^{(n,k)}(H_s)}
\end{array}\right\}.
\end{align*}
Intuitively, \(\hat\Psi_{t,\underline{a}_t}^{(n,k)}\) serves as an augmented estimator of the
true conditional mean potential outcome
\(\gamma_{t,\underline{a}_t}^o(H_t)\). In particular, if all nuisance functions are
correctly specified, then
$\mathbb{E}[\hat\Psi_{t,\underline{a}_t}^{(n,k)} \mid H_t]
= \gamma_{t,\underline{a}_t}^o(H_t).$
In addition, one can verify that
\(
\varphi_{\overline{a}}(O;\psi_{\overline{a}},\hat\alpha_{\overline{a}}^{(n,k)}) =
\hat\Psi_{0,\overline{a}}^{(n,k)} - \psi_{\overline{a}}.
\)
This representation motivates the construction of the estimator $\hat\psi_{\overline{a}}^{(n)}$ in Equation (ref).
An estimator of nuisance function $\hat{\alpha}_{t,\overline{a}}^{(n,k)}$ is referred to as a DR-Learner(or IF-Learner), and the theoretical properties of the nuisance vector $\hat{\alpha}_{t,\overline{a}}^{(n,k)}$ have been well-established in curth2021estimatingstructuraltargetfunctions,nie2020quasioracleestimationheterogeneoustreatment,foster2023orthogonalstatisticallearning, morzywolek2024weightedorthogonallearnersheterogeneous.
Next, we establish that, under the convergence rates required
for nuisance estimators, $\hat{\psi}_{\overline{a}}^{(n)}$ is asymptotically normal,
and its variance estimator $\hat{\sigma}_{\overline{a}}^{(n)}$ is consistent.
\begin{theorem}[Asymptotic normality of $\hat\psi_{\overline{a}}^{(n)}$]
Under Assumptions (ref)--(ref), suppose that for each \(t = 0, \ldots, T\), \(Z_t\) serves as an additive IV for \(A_t\), and that \(\pi_t(Z_t, H_t) = Z_t\) is an RWF for \(A_t\).
Further, for each \(t = 0, \ldots, T\) and \(k = 1, \ldots, K\), suppose that the following rate condition holds for the nuisance functions $\hat\alpha_{t,\overline{a}}^{(n,k)}$ defined in Algorithm (ref):
\begin{align*}
\left\{\begin{array}{l}s
\|\hat\kappa_{t,a_t}^{(n,k)}- \kappa_{t,a_t}^{o}\|_2
\times\|\hat\gamma_{t,\underline{a}_t}^{(n,k)}-\gamma_{t,\underline{a}_t}^o\|_2\\
+\|\hat\rho_t^{(n,k)} - \rho_t^o\|_2\times
\|\hat\eta_{t,\underline{a}_t}^{(n,k)}-\eta_{t,\underline{a}_t}^o\|_2\\
+\|\hat\rho_t^{(n,k)} - \rho_t^o\|_2\times
\|\hat\delta_{t,a_t}^{(n,k)}- \delta_{t,a_t}^{o}\|_2
\end{array}\right\}=o_p(n^{-1/2}).
\end{align*}
Furthermore, assume that for any fixed $k$ and time $t$,
\(\mathbb{E}[\|\hat\alpha_{t,\overline{a}}^{(n,k)}- \alpha_{t,\overline{a}}^o\|_2^2]=o(1)\).
Then $\sqrt{n}\{\hat\psi_{\overline{a}}^{(n)}-\psi_{\overline{a}}^o\}/\sigma_{\overline{a}}^o$ converges in distribution to $\mathcal{N}(0,1)$, where $\hat\psi_{\overline{a}}^{(n)}$ is defined in Equation (ref), and $(\sigma_{\overline{a}}^o)^2:=\mathbb{E}[\varphi_{\overline{a}}(O;\psi_{\overline{a}}^o,\alpha_{\overline{a}}^o)^2]$.
In addition, $\hat \sigma_{\overline{a}}^{(n)}$ converges in probability to $\sigma_{\overline{a}}^o$.
\end{theorem}
Notably, the RWF \(\pi_t(Z_t,H_t)\) can be can adaptively selected as demonstrated in Section (ref) to obtain an efficient estimator of \(\psi_{\overline{a}}^o\). A detailed discussion of this adaptive selection is provided in the Supplementary Material.
Intuitively, a reliable candidate for adaptive RWF is the conditional probability function \(\pi_t^o(Z_t,H_t) = \Pr(A_t = a_t \mid Z_t,H_t)\), which directly characterizes the treatment assignment mechanism.
Finally, although we focus on the evaluation of static treatment rules, the proposed methods can be extended to facilitate the evaluation of dynamic treatment rules, which is also demonstrated in the Supplementary Material. This extension would allow the incorporation of time-varying treatment strategies and enable more comprehensive assessments of treatment effects over time, enhancing the applicability of our approach to dynamic decision-making processes.
\section{Simulations}
\subsection{Point exposure scenario}
In this section, we conduct simulation studies with a binary treatment and a continuous IV to illustrate the asymptotic results established in Section (ref).
We generate $U$ and $L$ from the uniform distribution $U(-1,1)$ independently. The continuous IV is constructed as $Z = L + \sin(3L) + 2\epsilon_Z$, where $\epsilon_Z$ is an exogenous error term drawn from the standard normal distribution. Let “Ber” denote the binomial distribution. The binary treatment $A$ is generated under the following two designs:
\begin{enumerate}[label=(A\arabic*),leftmargin=1cm]
• $A \sim \mathrm{Ber}\big(1,0.7\Phi(-2Z + 2L) + 0.3\Phi(3U - L)\big)$, where $\Phi$ denotes the cumulative distribution function of the standard normal distribution.
• $A \sim \mathrm{Ber}\big(1,\{1 + \exp(-(Z - L + U))\}^{-1}\big)$.
\end{enumerate}
It is straightforward to verify that $Z$ is an additive IV for $A$ in (ref), while in (ref) the additive IV conditions are violated.
Next, we independently generate $\epsilon_Y \sim N(0,1)$. The outcome $Y$ is then generated according to the following two models:
\begin{enumerate}[label=(Y\arabic*),leftmargin=1cm]
• $Y = Y(A) = 2U - 2L + 4AL + \epsilon_Y$.
• $Y = Y(A) = (1 - A)\{3\cos(2U) - 3\cos(2L)\} + A\{3\sin(2U) + 2L\} + \epsilon_Y$.
\end{enumerate}
The condition $Y(1) - Y(0)\perp\!\!\!\perp U \mid L$ is satisfied in (ref) but violated in (ref). It is straightforward to verify that the ATEs in (ref) and (ref) are both zero.
We set the sample size $n=5000$ and use $K=2$ folds for cross-fitting. Nuisance functions are estimated via spline methods implemented in the \texttt{R} package \texttt{mgcv} wood2011fast.
The results are summarized in Table (ref).
We report the finite sample behavior of the ATE estimators
\(\hat\psi_{\pi}^{(n)}\) and \(\hat\psi_{*}^{(n)}\), defined in
Equations (ref) and
(ref), respectively.
The corresponding results are also presented for \(\mathbb{E}[Y(1)]\) (treated) and \(\mathbb{E}[Y(0)]\) (control).
Under the model misspecification settings (ref), (ref), the ATE estimator tends to exhibit larger bias due to the violation of Equation (ref). We also find that the adaptive weighting method generally yields a smaller variance, suggesting superior efficiency.
\begin{table}
\caption{Simulation results under four DGPs. In each setting, the estimation is repeated 1000 times to calculate the average bias (Bias), empirical standard deviation (SD), average estimated standard error (SE), and the empirical 95% coverage rate (CR).}
\begin{tabular}{cc| ccc |ccc}
\toprule
\multirow{2}{*}{DGP}&\multirow{2}{*}{Metric}&
\multicolumn{3}{c}{Adaptive weight}&\multicolumn{3}{|c}{Prespecified ($\pi(Z,L)=Z$)}\\
&& Treated & Control & ATE & Treated & Control & ATE\\
\midrule
\multirow{4}{*}{(ref), (ref)}
&Bias &.0011&.0011&.0019& .0018 &.0032&.0028 \\
&SE &.0543&.0545&.0807& .0626 &.0626&.0920 \\
&SD &.0577&.0581&.0844& .0612 &.0648&.0922 \\
&CR &.927&.938&.939& .947 &.935&.950 \\
[5pt]
\multirow{4}{*}{(ref), (ref)}
&Bias &.0023&.0078&.0064& .0043 &.0044&.0012 \\
&SE &.0864&.0611&.1060& .1006 &.0691&.1220 \\
&SD &.0906&.0637&.1092& .0987 &.0703&.1200 \\
&CR &.926&.946&.941& .951 &.951&.952 \\
[5pt]
\multirow{4}{*}{(ref), (ref)}
&Bias &.0006&.0020&.0013& .0006 &.0013&.0003 \\
&SE &.0554&.0561&.0820& .0565 &.0569&.0833 \\
&SD &.0547&.0575&.0830& .0563 &.0557&.0828 \\
&CR &.954&.934&.938& .956 &.956&.944 \\
[5pt]
\multirow{4}{*}{(ref), (ref)}
&Bias &.0000&.0318&.0315& .0027 &.0269& .0301 \\
&SE &.0886&.0615&.1079& .0908 &.0621& .1100 \\
&SD &.0891&.0619&.1089& .0915 &.0618& .1110 \\
&CR &.948&.912&.936& .955 &.924& .941 \\
\bottomrule
\end{tabular}
\end{table}
\subsection{Longitudinal scenario}
To simplify the setting, we fix the time horizon at \( T = 1 \). For each time point \( t = 0,\,1 \), we independently generate noise terms \( \epsilon_{U_t} \), \( \epsilon_{L_t} \), \( \epsilon_{Z_t} \), and \( \epsilon_Y \) from standard normal distributions. Based on these, the variables are simulated according to the following DGP:
\begin{align*}
&L_0 \sim 1.5\epsilon_{L_0},\quad U_0 \mid H_0 \sim 1.5\epsilon_{U_0};\\
&Z_0 \mid H_0,U_0 \sim 0.3 L_0 + \sin(1.5 L_0)+ 2\epsilon_{Z_0};\\
&A_0 \mid H_0,U_0,Z_0 \sim \text{Ber}(1,\,0.7\Phi(-2Z_0+0.6L_0) + 0.3\Phi(3U_0-L_0));\\
&L_1 \mid H_0,U_0,Z_0,A_0 \sim (A_0-0.5)+0.5L_0+0.3U_0+0.5 \epsilon_{L_1};\\
&U_1 \mid H_1,U_0 \sim (A_0-0.5)+0.5U_0+0.3L_1+0.5 \epsilon_{U_1};\\
&Z_1 \mid H_1,\overline{U}_1 \sim 0.5L_1-0.5(A_0-0.5)-0.3Z_0+2\epsilon_{Z_1};\\
&A_1 \mid H_1,\overline{U}_1,Z_1 \sim \text{Ber}(1,\,0.7\Phi(-2Z_1+L_1) + 0.3\Phi(3U_1-L_1));\\
&Y \mid H_1,\overline{U}_1,Z_1,A_1 \sim (A_1-0.5)+2L_1+U_1+0.5\epsilon_{Y}.
\end{align*}
Notably, $Z_0$ and $Z_1$ serve as additive IVs for $A_0$ and $A_1$, respectively.
Data are generated using the \texttt{R} package \texttt{simcausal} sofrygin2017simcausal. The sample sizes are set to $2000$ and $5000$, with cross-fitting performed using $K=2$ folds. Nuisance functions are estimated via spline methods implemented in the R package \texttt{mgcv}.
Because the true treatment effects are not analytically available under the constructed DGPs, we estimate them using 100{,}000 samples by simulating potential outcomes under modified data-generating processes, where pairs of treatment assignments, $(A_0,A_1)$, are set deterministically to \((0,0)\), \((1,1)\), \((0,1)\), \((1,0)\), \((A_0,1),\) and \((A_0,0)\) (We use the notation \((A_0,0)\) to denote a dynamic treatment regime that follows the natural treatment rule for \(A_0\) while fixing \(A_1 = 0\)).
Table (ref) summarizes the simulation results adopting Algorithm (ref) under different intervention strategies and sample sizes. Across all scenarios, the estimators show small bias, and the estimated standard errors closely match the empirical standard deviations, indicating accurate variance estimation. As expected, increasing the sample size from 2000 to 5000 reduces variability and tightens confidence intervals. The empirical coverage rates are generally close to the nominal 95% level, although a modest decline is observed in certain intervention settings. Overall, the results demonstrate that the proposed method performs well and provides reliable inference in most cases.
\begin{table}
\caption{Simulation results under three DGPs. In each setting, the estimation is repeated 500 times to calculate the average bias (Bias), average estimated standard error (SE), empirical standard deviation (SD), and the empirical 95% coverage rate (CR).}
\begin{tabular}{ccc|cccccc}
\toprule
&\multirow{2}{*}{Size}&
\multirow{2}{*}{Metric}&\multicolumn{6}{c}{Intervention}\\
&&& $(0,0)$ & $(0,1)$ & $(1,0)$ & $(1,1)$ & $(A_0,0)$ & $(A_0,1)$ \\
\midrule
& \multirow{4}{*}{2000}
& Bias &.0234&.0071&.0147&.0189&.0016&.0033\\
& & SE &.6090&.5267&.5266&.5960&.1401&.1403\\
& & SD &.6102&.5382&.5203&.6343&.1395&.1460\\
& & CR &.955&.959&.960&.938&.938&.929\\
[5pt]
& \multirow{4}{*}{5000}
& Bias &.0134&.0078&.0057&.0155&.0057&.0037\\
& & SE &.2143&.1882&.1962&.2115&.0658&.0651\\
& & SD &.2211&.1949&.1870&.2004&.0607&.0612\\
& & CR &.940&.930&.942&.956&.956&.950\\
\bottomrule
\end{tabular}
\end{table}
\section{Empirical illustration}
We perform a policy analysis using the dataset provided in the Supplementary Material of han2024optimal.
Schooling and post-school training are two central interventions that influence labor market outcomes such as earnings and employment ashenfelter2010handbook.
To enable this analysis, han2024optimal merged data from the Job Training Partnership Act (JTPA) Title II with additional sources on high school (HS) diploma, thereby constructing a dataset suitable for evaluating the effects of HS diplomas and subsidized job training as sequential treatments.
The final sample comprises 9,223 individuals.
We now describe the key features of this dataset.
Let $A_0$ denote whether an individual obtains an HS diploma, and let $A_1$ indicate participation in the job training program.
We define $L_0$ as sex (a baseline covariate) and $L_1$ as pre-program earnings, which act as two time-varying confounders.
The initial treatment $A_0$ affects subsequent pre-program earnings $L_1$, and the assignment to the job training program $A_1$ may depend on $L_1$.
The instruments are $Z_0$, the number of high schools per square mile, and $Z_1$, a random assignment to job training program.
Our target outcome $Y$ indicates whether an individual's terminal earnings exceed the empirical median.
We consider the dynamic treatment regime (DTR) $[g_0, g_1] \in \{0,1,\text{x}\} \times \{0,1,d^+,d^-\}$.
For DTR $g_0$, the value `0' assigns $A_0=0$, `1' assigns $A_0=1$, and `x' assigns the natural selection rule (i.e., the observed assignment). For the treatment rule $g_1$, the value `0' assigns $A_1=0$, `1' assigns $A_1=1$, `$d^+$' assigns $A_1=1$ only when $L_1$ is below the 80% quantile, and `$d^-$' assigns $A_1=1$ only when $L_1$ is above the 80% quantile.
The target is to estimate $\mathbb{E}[Y(g_0(H_0),g_1(H_1))]$.
We use the spline methods in the R package \texttt{mgcv} to estimate the nuisance parameters and calculate the bootstrapped estimated mean and standard deviation based on 1,000 replications.
\begin{table}
\caption{
Estimations for different types of DTRs are repeated 1,000 times using the bootstrap method.
}
\begin{tabular}{c| cccc| cccc| cccc}
\toprule
DTRs
& \texttt{00} & \texttt{01} & \texttt{0$d^+$} & \texttt{0$d^-$}
& \texttt{10} & \texttt{11} & \texttt{1$d^+$} & \texttt{1$d^-$}
& \texttt{x0} & \texttt{x1} & \texttt{x$d^+$} & \texttt{x$d^-$} \\
\midrule
EST &.292 &.321 &.343 &.264 &.653 &.668 &.732 &.589 &.485 &.550 &.548 &.487 \\
SE &.124 &.089 &.094 &.120 &.155 &.122 &.132 &.147 &.017 &.012 &.013 &.016 \\
SD &.108 &.074 &.081 &.103 &.129 &.089 &.102 &.126 &.017 &.011 &.012 &.015 \\
\bottomrule
\end{tabular}
\end{table}
The results are demonstrated in Table (ref).
The table presents estimated mean of mean potential terminal income (EST), mean estimated standard errors (SE), and empirical standard deviations (SD) for various DTRs.
Specifically, the DTRs are denoted by combinations of values in the set $\{0,1,\text{x}\}$ for the first treatment $g_0$ and $\{0,1,d^+,d^-\}$ for the second treatment $g_1$.
The treatment values and their respective effects are summarized for 12 distinct combinations.
We observe that the SE and SD for \{\texttt{x0}, \texttt{x1}, \texttt{x$d^+$}, \texttt{x$d^-$}\} are smaller compared to the other DTRs.
This is likely due to the relatively weak correlation between \( Z_0 \) and \( A_0 \).
Furthermore, the SE and SD values for the same DTR are fairly consistent with each other, demonstrating the validity of the variance estimator in our algorithm.
Considering the EST, the DTRs $\{\texttt{10}, \texttt{11}, \texttt{1d}^+, \texttt{1d}^-\}$ generally yield higher estimates than their counterparts $\{\texttt{00}, \texttt{01}, \texttt{0d}^+, \texttt{0d}^-\}$. The natural treatment rules $\{\texttt{x0}, \texttt{x1}, \texttt{xd}^+, \texttt{xd}^-\}$ produce estimates that lie between these two groups. This pattern suggests that obtaining an HS diploma has a positive effect on income.
Meanwhile, the estimated average terminal income is higher for DTRs $\{\texttt{0d}^+, \texttt{1d}^+, \texttt{xd}^+\}$, which assign the program only to low-earning individuals, than for DTRs $\{\texttt{00}, \texttt{10}, \texttt{x0}\}$ or $\{\texttt{01}, \texttt{11}, \texttt{x1}\}$. In contrast, DTRs $\{\texttt{0d}^-, \texttt{1d}^-, \texttt{xd}^-\}$, which assign the program only to high-earning individuals, yield lower estimates. This pattern aligns with the findings of han2024optimal, suggesting that restricting the job-training program to low-earners positively affects terminal income.
\section{Supplementary Material}
The supplementary material contains the additional results on dynamic treatment regimes, multiplicative IVs, connection to handling continuous instruments via discretization, identification with continuous treatments, and proofs.
The codes of the simulation results and real data analysis are publicly available at \url{https://github.com/chensy123-sys/Additive-IV}.
\putbib[bibfile_main]
bibunit