# Empirical Gateaux Derivatives for Causal Inference

Michael Jordan<EMAIL_ADDRESS>University of California, Berkeley
Yixin Wang<EMAIL_ADDRESS>University of Michigan  Angela Zhou
111Alphabetical order<EMAIL_ADDRESS>University of Southern
California, Marshall School of Business

###### Abstract

We study a constructive algorithm that approximates Gateaux derivatives for
statistical functionals by finite-differencing, with a focus on causal
inference functionals. We consider the case where probability distributions
are not known a priori but also need to be estimated from data. These
estimated distributions lead to empirical Gateaux derivatives, and we study
the relationships between empirical, numerical, and analytical Gateaux
derivatives. Starting with a case study of the interventional mean (average
potential outcome), we instantiate the exact relationship between finite-
differences and the analytical Gateaux derivative. We then derive requirements
on the rates of numerical approximation in perturbation and smoothing that
preserve the statistical benefits of one-step adjustments, such as rate-
double-robustness. We then study more complicated functionals such as dynamic
treatment regimes and the linear-programming formulation for policy
optimization in infinite-horizon Markov decision processes. The newfound
ability to approximate bias adjustments in the presence of arbitrary
constraints illustrates the usefulness of constructive approaches for Gateaux
derivatives. We also find that the statistical structure of the functional
(rate-double robustness) can permit less conservative rates of finite-
difference approximation. This property, however, can be specific to
particular functionals, e.g. it occurs for the average potential outcome
(hence average treatment effect) but not the infinite-horizon MDP policy
value.

##  1 Introduction

Inferential targets in causal machine learning are often statistical
functionals of the data distribution, for example an average treatment effect
or an average policy value for infinite-horizon offline reinforcement
learning. Estimation of these statistical functionals can be vulnerable to the
first-stage bias introduced by the estimation of nuisance functions such as
outcome regressions or transition probabilities. However, by leveraging the
causal structure, it is possible to derive bias-adjusted estimators that
weaken the need for accurate estimation of these nuisance functions. The
celebrated doubly robust estimator is one such bias adjustment, and there are
general frameworks for deriving such estimators, including semiparametric
efficiency, double robustness, Neyman orthogonality, and “debiased/double”
machine learning, as summarized in Bickel et al. (1993); Robins et al. (1994);
Newey (1994a); Chernozhukov et al. (2016, 2018); Hines et al. (2022); Kennedy
(2016, 2022). These general frameworks can also be applied to more complicated
functionals in longitudinal causal inference and offline reinforcement
learning (also known as offline policy evaluation/learning) (Chakraborty and
Moodie, 2013; Van der Laan and Robins, 2003; Jiang and Li, 2016; Thomas and
Brunskill, 2016; Bibaut et al., 2019; Kallus and Uehara, 2020, 2022; Tang et
al., 2019; Kallus and Zhou, 2020, 2022).

A drawback of these general frameworks is that significant analytic effort is
often necessary to obtain concrete estimators that are appropriate for
particular situations. Indeed, a practitioner may be interested in novel
variants of an estimand, or be working within a constrained class of
functionals that are appropriate for a particular class of problems. For
example, in constrained Markov Decision Processes (MDPs) (Altman, 1999), or
optimization-based estimators more broadly, a practitioner can easily avail
themselves of custom convex-optimization solvers that are computationally
efficient and easy to deploy to add additional constraints. Unfortunately,
such choices may require case-specific re-analysis to establish the desired
statistical properties, and indeed deriving the actual estimator that yields
bias adjustment may require significant analytical effort. It is therefore
important to develop a suite of constructive, numerical or algorithmic methods
that yield desired forms of bias adjustment. Such a suite would be
complementary to analytic derivations.

To this end, we develop off-the-shelf procedures for estimating the
statistical functionals used for bias adjustment in causal inference. We do so
via numerical approximation of the Gateaux derivatives that underlie the
general analytic frameworks, building on prior work in this area. We focus on
interventional effects. Newey (1994b) suggested numerical differentiation to
estimate Gateaux derivatives of functionals unavailable in closed form but
approximated by computational procedures, such as the solution to a system of
differential equations; this approach is later used in Hausman and Newey
(1995). Recent pedagogical surveys of Hines et al. (2022); Kennedy (2022) also
emphasize the usefulness of Gateaux derivatives for deriving influence
functions, which is valid for discrete data and generalizes for continuous
distributions via a smoothing argument discussed in Ichimura and Newey (2015);
Carone et al. (2018). In this spirit, we build on the line of work of Newey
(1994b); Frangakis et al. (2015); Carone et al. (2018); Ichimura and Newey
(2015), which proposes a constructive procedure by numerical approximation of
a Gateaux derivative (in the sense of Hampel (1974); Huber (2004)) by finite
differences.

Numerical differentiation can be a useful paradigm because evaluating the one-
step estimator reduces to multiple evaluations of the black-box statistical
functional. In other areas such as optimization, numerical differentiation
underpins common subroutines that enable analysts to use optimization
algorithms without specialized training in a particular modeling language or
paradigm.222For example, many data scientists use optimization through the
scipy.optimize or r.optim packages: if they do not provide a gradient function
callback, in some configurations, it is by default approximated by numerical
differentiation.  Although more advanced paradigms can compute exact
gradients, they can require specialized training to use or implement.

In addition, the numerical differentiation viewpoint connects the conceptual
interpretation of influence functions as qualitative sensitivity analysis
(i.e. sensitivity of a functional to perturbations in the distribution) to the
use of influence functions in causal inference. Typical presentations of the
latter framework are not concretely connected to the former notion, where
influence functions are well-known in robust statistics. For example, this
former notion was applied in the machine learning setting to probe machine
learning models (Koh and Liang, 2017). We provide a concrete connection from
the viewpoint of exactly quantifying the finite difference approximation
error, which can be interesting independent of computerized estimation for
algorithmic, pedagogical, or conceptual reasons.

In this paper, we build on the previously mentioned line of work on the use of
numerical derivatives to approximate influence functions, focusing on the
implications for statistical estimation. In more detail, we provide exact
characterizations of empirical Gateaux derivatives, in contrast to numerical
derivatives computed based on oracle knowledge of probability distributions.
This work is focused on exact characterizations that can inform statistical
and computational trade-offs. Although we expect that some of these
characterizations may be apparent to experts, we also hope that this level of
concreteness can be helpful to non-experts. We use these characterizations to
study statistical and computational trade-offs: how statistical properties of
functionals in causal inference can weaken computational approximation
requirements, as well as more complicated functionals. These trade-offs could
eventually be practically relevant because the finer-scaled discretization
required for lower approximation error can incur fundamental issues with
floating-point computation that prohibit further gains in accuracy (Fornberg,
1981; Sauer, 2012).

Our contributions are as follows. We begin with a case study of the average
potential outcome. The average treatment effect is the difference of the
average potential outcomes, and our analysis applies also to the average
treatment effect. We discuss the treated mean for simplicity; by symmetry the
same results hold for the control mean. In this setting, we characterize the
exact relationship between finite-differences and the analytical Gateaux
derivative. This characterization helps concretize the proposal of empirical
Gateaux derivatives by allowing us to study the rates of numerical
approximation that can preserve statistical performance. Finally, with these
concretizations in hand, we turn to more complicated motivating examples
including a finite-horizon dynamic treatment regime and infinite-horizon
reinforcement learning. We illustrate with the case of policy learning in
finite-state, finite-action infinite-horizon Markov decision processes how
this constructive approach to Gateaux derivatives can directly extend also to
the case of custom constraints.

##  2 Problem Setup: From Numerical to Empirical Gateaux derivatives

We begin by introducing our first example, that of the mean potential outcome,
and the canonical doubly-robust estimator. We then briefly introduce the key
objects in the more general framework of orthogonalizing or debiasing
statistical functionals by Gateaux derivative adjustments via the one-step
estimator, as well as numerical approximation of these via perturbed black-box
evaluations of the statistical functional via plug-in estimation. After
introducing these frameworks, we clarify what we mean by empirical rather than
numerical Gateaux derivatives and describe our research questions.

##### Setup.

Denote $O\sim P$ as the observations following distribution $P$ belonging to a
statistical model $\mathcal{M}$. We will focus on targeting the estimation of
statistical functionals $\Psi(P)$, where $\Psi:\mathcal{M}\to\mathbb{R}^{q}$
is pathwise differentiable. Throughout, we take as given the assumptions that
grant causal identifiability (i.e. overlap/positivity,
unconfoundedness/ignorability); our concrete examples are all of so-called
backdoor adjustment under these assumptions.

##### The mean potential outcome and augmented inverse propensity weighting.

We very briefly overview the celebrated doubly robust estimate of the mean
potential outcome. The average treatment effect is simply the difference of
these doubly robust estimators for the treated and control means; for
simplicity, we discuss the treated mean only.

######  Example 1 (Mean potential outcome).

The observation $O=(X,A,Y)$ includes covariate $X$, treatment $A\in\\{0,1\\}$,
and outcome $Y\in\mathbb{R}.$ The statistical functional of the mean potential
outcome is:

$\textstyle\Psi(P)=\mathbb{E}[Y(1)]=\mathbb{E}[\mathbb{E}[Y\mid A=1,X]].$

Plug-in estimation by the conditional expectation is termed regression
adjustment. The functional is also identified by inverse propensity weighting,
$\textstyle\Psi(P)=\mathbb{E}[Y(1)]=\mathbb{E}\left[Y\frac{\mathbb{I}[A=1]}{e(X)}\right]$,
where $e(x)=P(A=1\mid X=x)$ is called the propensity score. (We omit explicit
discussion of assumptions such as consistency and ignorability, under which
these identification results hold).

Given two identifying functionals for the mean potential outcome, a natural
question is: can these identification approaches be combined in some way to
improve estimation? The doubly-robust estimator exactly does so (Robins et
al., 1994).

$\textstyle\mathbb{E}[Y(1)]=\mathbb{E}\left[\frac{\mathbb{I}\left[A=1\right]}{p(A=1\mid
X)}(Y-\mathbb{E}[Y\mid A=1,X])+\mathbb{E}[Y\mid A=1,X]\right]$

Appealing properties (mixed-bias and rate-double-robustness) can easily be
verified from its functional form. Although the form of AIPW could be deduced
from the point of view of combining different identification approaches, more
broadly, it is a canonical example of one-step adjustment by influence
functions.

##### Influence function and Gateaux derivative.

A functional $\Psi$ is Gateaux differentiable at P if
$\Psi^{\prime}(H;P)\triangleq\frac{d\Psi(P+\epsilon
H)}{d\epsilon}|_{\epsilon=0}$ exists and is linear and continuous in $H$ for
all $H$.333Another statement of the definition is that $\Psi$ is Gateaux
differentiable at $P$ if there exists $\phi_{P}(H)$ such that
$\phi(\theta+\epsilon H)-\Psi(P)=\epsilon\phi_{P}(H)+o(t),$ as $t\to 0$ (Van
der Vaart, 2000, p.296). This statement coincides with the previous when the
functional maps to the real line.

Hampel (1974) and Huber (2004) define the influence function $\phi(O;P)$ as a
Gateaux derivative with respect to perturbations $H=\delta_{o}-P$ toward a
degenerate Dirac measure $\delta_{o}$ satisfying the identity that $\int
o\delta_{o^{\prime}}(o)=o^{\prime}$,

$\phi(o;P)=\frac{\mathrm{d}\Psi\left(P+\epsilon(\delta_{o}-P)\right)}{\mathrm{d}\epsilon}\Big{|}_{\epsilon=0}.$

##### Related work: Numerical Gateaux derivatives.

The influence function definition via perturbation is not quite enough for
semiparametric inference with continuous distributions, but Carone et al.
(2018) and Ichimura and Newey (2015) show that it suffices to restore absolute
continuity in semiparametric models by taking an additional limit as a
smoothed perturbation converges to the point mass. That is, it suffices to
replace the Dirac delta measure $\delta_{o}$ with a smoothed distribution
$\tilde{\delta}_{{o}^{\prime}}^{\lambda}(o)$, a smoothed delta function
perturbed in the direction of a generic observation ${o}^{\prime}$, evaluated
at $o$, dominated by $P$, such that $\int
g(u)\mathrm{d}\tilde{\delta}_{{o}^{\prime}}^{\lambda}(u)=1.$ A common choice
for $\tilde{\delta}_{o^{\prime}}^{\lambda}$ uses kernel functions
$\tilde{\delta}_{o^{\prime}}^{\lambda}(o)=K_{\lambda}(o-o^{\prime})$ , where
$K_{\lambda}$ is an appropriately bandwidth-normalized kernel function and
notation is overloaded to also designate product kernels. The perturbed
distribution in the direction of a generic observation $o_{i}$ is

$P_{\epsilon,\lambda}^{o_{i}}=P_{\epsilon,\lambda}^{{i}}=(1-\epsilon)P+\epsilon\tilde{\delta}_{o_{{i}}}^{\lambda}.$

Frangakis et al. (2015); Carone et al. (2018) both discuss numerical
approximation by finite-differences: for a fixed $\epsilon,\lambda$ with
$\epsilon\ll\lambda$, evaluate the numerical derivative
$\hat{\phi}(o;P)=\epsilon^{-1}({\Psi(P_{\epsilon,\lambda})-\Psi(P))}{}.$
Further, Hines et al. (2022); Kennedy (2022) highlight influence function
derivations based on the analytic Gateaux derivative.

Carone et al. (2018, Thm. 1) establishes that the Gateaux derivative is the
limit as the paths and smoothing parameter, $\epsilon$ and $\lambda$, tend to
$0$. They generally assume all distributions are known in the analysis and
code. Carone et al. (2018, Thm. 5) provided generic approximation rates in
$(\epsilon,\lambda)$ for multi-point finite difference scheme using results
from the finite difference literature; these establish sufficient conditions
for approximation error that ensure validity of the numerical approximation,
but can be conservative. For a two-point approximation with a uniform kernel,
the rate from Carone et al. (2018, Thm. 5) is
$\left(\epsilon\lambda^{-d}\right)+\lambda^{-2}$. Our analysis later in this
paper weakens these approximation error requirements.

##### Empirical Gateaux derivatives: Plug-in estimation of the functional.

In the context of statistical estimation, the underlying probability
distributions also need to be estimated. Therefore, to ensure validity of the
numerical derivative approach for a general class of functionals, we consider
plug-in estimation of the statistical functional by conducting probability
density estimation of the collection of joint probability estimates
$\tilde{p}(o)$ that comprise $\tilde{P}$. For example, a plug-in
representation of the mean potential outcome is

$\Psi(\tilde{P})=\int\int
y\frac{\tilde{p}(y,A=1,x)}{\tilde{p}(A=1,x)}\tilde{p}(x)\mathrm{d}y\mathrm{d}x.$

Hence, we distinguish between the empirical Gateaux derivative computed by
numerical approximation with estimated probability densities (e.g.
$\tilde{P}(Y,A=1,X)$), the numerical Gateaux derivative obtained by finite-
differencing but with the true probability densities (e.g. $P(Y,A=1,X)$), and
the analytic Gateaux derivative.

$\displaystyle{{\tilde{\phi}}}(O_{i})$
$\displaystyle=\epsilon^{-1}\left({\Psi(\tilde{P}_{\epsilon}^{i})-\Psi(\tilde{P})}{}\right)$
empirical derivative at smoothed and estimated distributions,
$\displaystyle\hat{\phi}(O_{i})$
$\displaystyle=\epsilon^{-1}\left({\Psi(P_{\epsilon}^{i})-\Psi(P)}{}\right)$
numerical derivative at smoothed and true distributions,
$\displaystyle\phi(O_{i})$
$\displaystyle=\frac{\mathrm{d}}{\mathrm{d}\epsilon}\Psi(P_{\epsilon}^{i})\Big{|}_{\epsilon=0}$
analytical Gateaux derivative.

The one-step estimator is an immediate way to use computationally estimated
Gateaux derivatives for statistical estimators (Pfanzagl, 1990), justified by
the asymptotically linear representation (see Appendix A).

$\Psi_{n}=\Psi(\tilde{P})+\frac{1}{n}\sum_{{i}=1}^{n}{{\tilde{\phi}}}(O_{i}),\qquad\text{
where
}{{\tilde{\phi}}}_{i}(O_{i})=\frac{1}{\epsilon}(\Psi(\tilde{P}^{i}_{\epsilon})-\Psi(\tilde{P}))$
(1)

Although for simplicity we consider first-order finite difference methods in
the main text, multipoint methods for numerical differentiation can also be
used.

##### Problem statement.

About the empirical, numerical, and analytical derivatives, we answer:

Q1: How exactly does the empirical Gateaux derivative approximate the
analytical derivative?

Q2: What are the required rates of numerical approximation in perturbation and
smoothing $(\epsilon,\lambda)$ that preserve the beneficial statistical
properties of the constructed estimator?

In Section 4, we first address these questions in the simple case of our
running example: estimation of the mean potential outcome, where we can finely
compare the finite-difference estimator and the classical doubly robust
estimator to characterize statistical properties. In Section 5, we study more
complicated functionals to highlight more general applicability.

######  Example 2 (Example 1, continued).

The influence function (Robins et al., 1994; Hahn, 1998) for a mean potential
outcome is:

${\phi(O)=\frac{\mathbb{I}\left[A=1\right]}{p(A=1\mid X)}(Y-\mathbb{E}[Y\mid
A=1,X])+\mathbb{E}[Y\mid A=1,X]-\Psi(P)}.$

Hence the canonical doubly-robust estimator (augmented inverse propensity
weighting, or AIPW) can be viewed as a one-step adjustment.

We conclude our introductory remarks by noting some limitations of our
approach. We omit discussion of semiparametric efficiency since we focus on
computing influence functions without the discussion of model restrictions and
tangent spaces. Although Carone et al. (2018); Ichimura and Newey (2015) focus
on semiparametric models, we focus on Gateaux derivatives under a
nonparametric model; Carone et al. (2018) suggest one can handle the
semiparametric case by evaluating a Gateaux derivative and projecting onto the
model space directly. Hence our analysis can be viewed as considering the case
of nonparametric models; or analyzing the first step of a procedure for
semiparametric models.

Throughout, we assume the functional is pathwise differentiable. This
assumption does warrant caution, because this is a technical condition that
poses a fundamental barrier for efforts for “automation”; hence empirical
Gateaux derivatives are intended to augment but not supplant or replace
analytical expertise. However, under the (stronger) assumption of Frechet
rather than Gateaux differentiability, previous works give sufficient
conditions that imply the validity of these approaches. Ichimura and Newey
(2015, sec. 4) provides a generic sufficient condition for asymptotic
linearity assuming additionally that the Frechet derivative is Lipschitz,
while Kandasamy et al. (2015) show the remainder term is second order; both
under convergence rate conditions on the density inputs. See Section D.2 for
more discussion, including a class of problems where general guarantees are
possible, i.e. strongly convex stochastic programs.

##  3 Related Work

In this section, we discuss the most closely related directions but do not
attempt to provide a survey or overview of influence functions/semiparametric
statistics; see Bickel et al. (1993); Tsiatis (2006); Fisher and Kennedy
(2021); Kennedy (2022) for a broader overview or the appendix for more
detailed discussions. In Section D.1, for additional context, we also include
additional work on numerical derivatives in optimization and statistical
machine learning more broadly.

##### Numerical approximation of Gateaux derivatives.

Recent pedagogical reviews (Fisher and Kennedy, 2021; Hines et al., 2022;
Kennedy, 2022) emphasize the Gateaux derivative and its use as a heuristic
tool for deriving influence functions. While our initial characterization can
be viewed as a finite-difference approach to the former, we also consider more
complicated optimization functionals, with a focus on how dual variables
characterize the Gateaux derivatives therein. Carone et al. (2018) derive
high-level sufficient conditions for approximation error via general results
for finite-difference approximations. Relative to those general sufficient
conditions, we focus on specific improvements in approximation error in
different functionals. For the specific case of the mean potential
outcome/average treatment effect, Frangakis et al. (2015) use additional
structural knowledge. Namely, they interchange integration and differentiation
via analytical insight to obtain functional evaluations in terms of the
conditional expectation. For more complex functionals, making these
determinations is similar to analytically deriving the influence function,
i.e. less “computerized”. In general, it may be difficult to deduce the form
of the nuisance functions without using knowledge of the estimand or working
with graphical structure or algorithms. The approach in Appendix B of Bravo et
al. (2020) considers the specific case of debiasing moment conditions with
respect to first stages but also implicitly leverages this change in order of
integration.

##### Other work on computerizing causal inference with different approaches.

Several recent very exciting works in causal inference build towards
“automatic” or “computerized” semiparametric inference via different
approaches (Kandasamy et al., 2015; Bhattacharya et al., 2020; Chernozhukov et
al., 2021, 2020, 2022; Farrell et al., 2021; Jung et al., 2021b). To summarize
the difference here, we focus on concretizing an approach based on numerical
derivatives, which requires only the zeroth-order computation oracle of
perturbed functional evaluations, without additional nuisance estimation. This
different approach could be particularly helpful for incorporating debiased
causal inference in the design and analysis of algorithms. But of course,
conversely, it may not lead to as improved estimation performance as other
approaches that modify the mode of estimation or estimate additional nuisance
functions.

Kandasamy et al. (2015) outlines a general approach for Gateaux derivatives.
Chernozhukov et al. (2021) builds on variational characterizations of Ichimura
and Newey (2022) via Riesz representation for generalized regression
residuals. Chernozhukov et al. (2020, 2022) and other works further develop
“automatic” estimation based on this approach. Our references throughout to
Ichimura and Newey (2015) are primarily for the smoothing arguments required
in a semiparametric setting; meanwhile our last application to optimization is
not a generalized regression residual. Farrell et al. (2021) develops a
framework for modeling heterogeneity with deep networks, similar to a
functional coefficient model that admits an influence function computable by
automatic differentiation (or numerical differentiation of Hessians/Jacobians
of a given parametrization). Bhattacharya et al. (2020) focus on graphical
algorithms and characterizations of the efficient influence function for
identified latent variable models, for ATE-like functionals, given a causal
graph. Alaa and Van Der Schaar (2019) develop approaches for validating causal
machine learning methods with influence functions.

##### Influnce functions or Gateaux derivatives in other settings.

More broadly in machine learning, influence functions have been used for
qualitative sensitivity analysis or bias correction (Koh and Liang, 2017;
Giordano et al., 2019a, b; Wilson et al., 2020); further connections could be
interesting. However, the use of influence functions in the causal setting is
different because one computes influence functions of the statistical
functional with respect to perturbations in the distribution, rather than
perturbing the moment function that identifies the estimator (i.e
investigating changes in the neural network model under datapoint
perturbations to the loss minimization). Recently, in data-driven
optimization, finite-difference approaches for Gateaux derivatives have been
used to debias in-sample bias of the stochastic optimization functional in
general, see e.g. Gupta et al. (2021). Guo et al. (2022) builds on the finite
difference method of Ito et al. (2018) for debiasing optimization with causal
outcome coefficients but without exact characterization of the adjustment.

##  4 Case study: Empirical Gateaux Derivative of the Mean Potential Outcome

To illustrate the numerical approach to Gateaux derivatives, we first work
through the running example of the mean potential outcome and concretize
precisely how the numerical differentiation incurs error relative to the
analytical derivative, i.e. the doubly robust or augmented IPW (AIPW)
estimator.

###  4.1 Relating the empirical and analytical Gateaux derivatives

##### Setup.

We target the estimation of a mean potential outcome, Example 1. For the
smoothed perturbation, we specialize to kernel functions $K_{\lambda}(u)$ from
kernel density estimation satisfying that $\int_{-\infty}^{\infty}K(u)du=1$
(normalization to a probability density) and $K(-u)=K(u),\forall u$
(symmetry).444Examples include the Gaussian kernel with
$K(u)={({2\pi})^{-\frac{1}{2}}}e^{-u^{2}/2},$ or uniform with
$K(u)=\nicefrac{{1}}{{2}}\mathbb{I}\left[\left|u\right|\leq 1\right].$ We
generally consider
$K_{\lambda}(u)=\lambda^{-d}K(u/\lambda),K(x)=K\left(x_{1},\ldots,x_{d}\right)$.
The mean potential outcome is identified as follows,

$\textstyle\Psi(P)=\mathbb{E}[Y(1)]\textstyle=\mathbb{E}[\mathbb{E}[Y\mid
A=1,X]]=\int\int yp(y\mid A=1,x)\mathrm{d}y\mathrm{d}x=\int\int
y\frac{p(y,A=1,x)}{p(A=1,x)}p(x)\mathrm{d}y\mathrm{d}x.$ (2)

We obtain $\Psi(\tilde{P}_{\epsilon})$ by plugging in the $\tilde{p}$
probability density estimates. Note this is distinct from the plug-in
estimator in empirical process theory, which plugs in the empirical CDF; we
plug-in statistical estimates. See Algorithm 1 for a summary of the procedure.

##### Comparison of empirical Gateaux derivative and augmented inverse
propensity weighting.

We have argued that, in the context of statistical estimation, approximating
Gateaux derivatives requires plug-in estimation of the statistical functional
with estimated probabilities. Towards resolving Q1, we observe that the
perturbed probability densities induce the functional form of the nuisance
functions in the analytical derivative. It is in contrast to the analytical
derivative, which describes the estimand in terms of, e.g. conditional
expectations, and allows the analyst to choose the functional form for
estimating the nuisance function.

We will denote the $\tilde{P}$-induced conditional expectation as

${\mathbb{E}_{\tilde{P}}[Y\mid A=1,X=x]=\int y{\tilde{p}(y\mid
A=1,x)}\mathrm{d}y}.$

For example, the density-induced conditional expectation of kernel density
estimates is exactly Nadaraya-Watson regression.555For example,
$\textstyle\hat{\mu}^{NW}(x)=\int\frac{y\sum_{i=1}^{n}K_{h}\left(x-X_{i}\right)K_{h}\left(y-Y_{i}\right)}{\sum_{j=1}^{n}K_{h}\left(x-X_{j}\right)}\mathrm{d}y=\frac{\sum_{i=1}^{n}K_{h}\left(x-X_{i}\right)\int
yK_{h}\left(y-Y_{i}\right)\mathrm{d}y}{\sum_{j=1}^{n}K_{h}\left(x-X_{j}\right)}=\frac{\sum_{i=1}^{n}K_{h}\left(x-X_{i}\right)y_{i}}{\sum_{j=1}^{n}K_{h}\left(x-X_{j}\right)}.$
Given this, we study empirical derivatives with kernel density estimation. We
study KDEs because they provide a classic textbook example of density
estimates that induce conditional expectation estimators; our analysis is
modular in that other density estimation approaches can be used and their
rates of convergence would simply be substituted in Lemma 1. In Appendix D we
discuss other choices of density estimates: any “modern” generative modeling
approach/density estimation approach can be used for plug-in evaluation of the
statistical functional. Of course, all such nonparametric estimation
approaches suffer from the curse of dimensionality in general so additional
structure is required in general to satisfy the product-rate conditions.

Our characterization will show exactly how plug-in evaluation with finite
differences is nearly equivalent to evaluating the doubly-robust estimator
with a Nadaraya-Watson regressor as the conditional outcome regression for
$\mathbb{E}[Y\mid A=1,X]$ and a kernel density estimate of the propensity
score; up to an additive bias that arises due to the smoothed perturbation
distribution, and what we will term a smoothed nuisance evaluation. We next
discuss what this smoothed nuisance function is, and how we study its
approximation error relative to estimation error of the induced nuisance.

##### Integrating the smoothed perturbation induces smoothed nuisance
functions.

Integrating with respect to the smoothed perturbation distribution can be
interpreted as a smoothed evaluation of the resulting nuisance functions. We
define the $\tilde{\delta}_{x_{0}}^{\lambda}$-smoothed conditional expectation
which smooths the evaluation point, i.e. smooths evaluation around $x_{0}$
rather than precisely at $x_{0}$.

######  Definition 1 ($\tilde{\delta}_{x_{0}}^{\lambda}$-smoothed conditional
expectation).

$\textstyle\tilde{\mathbb{E}}_{P}[Y\mid A=1,X=x_{0}]=\int\mathbb{E}_{P}[Y\mid
A=1,X=u]\;\tilde{\delta}_{x_{0}}^{\lambda}(u)\;\mathrm{d}u$ (3)

Similar ideas appeared in the literature on non-regular inference or the
localization of a global functional (Newey, 1994b). For example, (van der Laan
et al., 2018) proposes smoothing a non-regular target estimand, such as a
dose-response curve or density evaluated at a point $x_{0}$, by conducting a
smoothed evaluation locally around $x_{0}$. In that work, a function evaluated
at $x_{0},$ $g(x_{0})$, is approximated with a kernel smooth over $x$ in a
neighborhood of $x_{0}$ with bandwidth $h$; so ${\Psi_{\lambda}(x_{0})=\int
K_{\lambda}(u-x_{0})g(u)\mathrm{d}u}.$ Jung et al. (2021a) uses a similar
smoothing/localization device but for a different problem in general. The
following lemma, similar to (van der Laan et al., 2018, Lemma 25.1),
summarizes convergence rate implications of these corresponding smoothed
nuisances; these are direct consequences of typical analyses of kernel density
estimation.

######  Lemma 1.

Let $\tilde{\mu}_{\epsilon}(X)=\tilde{E}_{\tilde{P}_{\epsilon}}[Y\mid A=1,X]$
and $\tilde{e}_{\epsilon}(X)=\tilde{p}_{\epsilon}(A=1,X)/\tilde{p}(X)$ denote
the nuisances induced by plug-in estimation of probability distributions
$\tilde{p}$. For any bounded function $g(x)$ in a Hölder class of degree
$\beta$,
$\textstyle(\int\tilde{\delta}_{x}^{\lambda}(u)g(u)\mathop{}\\!\mathrm{d}u-g(x))^{2}\mathop{}\\!\mathrm{d}x=O(\lambda^{\beta}).$
Assume that $\tilde{\delta}^{\lambda}$ is a product kernel in a Hölder class
of degree $\beta$, and
$\textstyle\int|{u}|^{\beta}|{\tilde{\delta}_{u}^{\lambda}}|\mathop{}\\!\mathrm{d}u<\infty$
and $\int u^{s}\tilde{\delta}_{u}^{\lambda}\mathop{}\\!\mathrm{d}u=0$ for
$s\leq\beta.$ Assume $Y$ is bounded. Then the perturbed nuisances satisfy the
following rates:

$\displaystyle\lVert\tilde{\mu}_{\epsilon}(X)-\mu(X)\rVert=\lVert\tilde{\mu}(X)-\mu(X)\rVert+O(\epsilon\lambda^{-d/2}),$
$\displaystyle\lVert\tilde{e}_{\epsilon}(X)-e(X)\rVert=\lVert\tilde{e}(X)-e(X)\rVert+O(\epsilon\lambda^{-d/2}).$

The perturbed nuisances are asymptotically consistent, so long as $\tilde{p}$
induces asymptotically consistent nuisances and we fix $\epsilon$ as a
sequence vanishing in $n$ at appropriate rate to counteract the growth in
$\lambda$ due to “roughness”, i.e. the variability of
$\tilde{\delta}^{\lambda}$. The difference between the function evaluation and
its smoothed nuisance arises from the bias analysis of kernel density
estimators and the dimension-dependence of ensuring absolute continuity, but
it is a non-stochastic, entirely deterministic argument. While the analysis
improves mildly in the dimension dependence on $\lambda$ compared to the
sufficient general analysis of Carone et al. (2018, Thm. 5), this still
highlights the generally unfavorable dependence of the numerical approach on
the dimension due to smoothing.

Algorithm 1 Empirical Gateaux derivatives

1:Inputs: $\tilde{P}$ probability density estimates, data $\\{O_{i}\\}_{1:n},$
perturbation and smoothing $(\epsilon,\lambda)$ parameters, functional $\Psi$

2:for $i=1,\dots n$ do

3:
$\tilde{P}_{\epsilon}^{i}=(1-\epsilon)\tilde{P}+\epsilon\tilde{\delta}_{o_{i}}^{\lambda}$

4:
${{\tilde{\phi}}}(O_{i})\leftarrow\epsilon^{-1}(\Psi(\tilde{P}_{\epsilon}^{i})-\Psi(\tilde{P}))$

5:end for

6:Output: one-step adjusted estimate from empirical Gateaux derivatives,
$\Psi(\tilde{P})+\frac{1}{n}\sum_{i}{{\tilde{\phi}}}(o_{i})$

Our first proposition establishes the exact approximation error induced by
finite-differencing.

######  Proposition 1 (Gateaux derivative of probability density
representation).

Consider perturbations in the direction of $o_{i}=(x_{i},a_{i},y_{i})$. Let
$\tilde{p}(x),\tilde{p}(A=1,x),\tilde{p}(y,A=1,x)$ denote kernel density
estimates. Let $\Psi(P)$ be as in Equation 2.

$\displaystyle\frac{\Psi(\tilde{P}_{\epsilon_{i}})-\Psi(\tilde{P})}{\epsilon}=\int\frac{\tilde{p}(x)\mathbb{I}\left[a_{i}=1\right]\tilde{\delta}_{x_{i}}^{\lambda}(x)}{\tilde{p}_{\epsilon}(A=1,x)}\left\\{\left(\int
y\;\tilde{\delta}_{y_{i}}^{\lambda}(y)\mathrm{d}y\right)-\mathbb{E}_{\tilde{P}}[Y\mid
A=1,X=x]\right\\}\;\mathrm{d}x$
$\displaystyle\qquad+\left(\tilde{\mathbb{E}}_{\tilde{P}_{\epsilon}}[Y\mid
A=1,X=x_{{i}}]-\Psi(\tilde{P})\right).$

We compare our expansion of Proposition 1 to the canonical AIPW estimator in
Example 2. We specialized to the case of a uniform kernel for
$\tilde{\delta}_{y_{i}}^{\lambda}(y)$ so that $\int
y\;\tilde{\delta}_{y_{i}}^{\lambda}(y)\mathrm{d}y=y_{i}$. We next observe that
the convolution with respect to the smoothing kernel introduces additive bias
relative to evaluating terms in the expression at the observation $o_{i}$.

######  Corollary 1.

When the perturbation observation is the observation datapoint
$O_{i}=(X_{i},A_{i},Y_{i})$,

$\displaystyle{{\tilde{\phi}}}(O_{i})$
$\displaystyle=\frac{\mathbb{I}\left[A_{i}=1\right]}{\tilde{p}_{\epsilon}(A=1\mid
X_{i})}\left(Y_{i}-{\mathbb{E}_{\tilde{P}}[Y\mid
A=1,X_{i}]}\right)+\left({\mathbb{E}}_{\tilde{P}_{\epsilon}}[Y\mid
A=1,X_{{i}}]-\Psi(\tilde{P})\right)+O(\lambda^{\beta}).$

Lemma 1 allows us to further simplify and conclude that the decomposition of
Proposition 1 is close to the canonical AIPW estimator up to the
$O(\lambda^{\beta})$ bias induced by smoothed evaluation.

###  4.2 Implications for estimation

Next, we use our characterization of Proposition 1, Lemma 1, and Corollary 1
to study the statistical properties of our computational/numerical
approximation. Corollary 1 arises from deterministic equivalences and allows
us to study the relationship to the typical “oracle” AIPW estimator; albeit
with induced $\epsilon$-perturbed nuisances.

An appealing property of bias-adjusted treatment effect estimation is rate-
double robustness: because we incur the product of convergence rates of the
nuisances, we may enjoy parametric $n^{-\frac{1}{2}}$-rate convergence of the
target functional while nuisances converge at slower rates, for example both
at $n^{-\frac{1}{4}};$ see Chernozhukov et al. (2018); Rotnitzky et al. (2021)
for more detail. We use the rate-double robustness property of the target and
Lemma 1 to infer the required rate conditions on $(\epsilon,\lambda)$ that
would retain the beneficial statistical properties of the canonical AIPW
estimator. The next result combines our previous characterizations with the
standard analysis of AIPW.

######  Assumption 1 (Regularity conditions).

Assume the following regularity conditions hold:

  1. (i)

$Y$ is bounded.

  2. (ii)

$p(y\mid A=1,x),p(A=1\mid x)$ are in Hölder classes of minimum degree $\beta$.

  3. (iii)

The estimates $p(y,a,x),p(a,x),p(x)$ belong to a Donsker function class.

  4. (iv)

Assume $\tilde{\mu},\tilde{e}$ satisfy the product-rate condition:

$\lVert\tilde{\mu}(X)-\mu(X)\rVert\times\lVert\tilde{e}(X)-e(X)\rVert=O_{p}(n^{-\frac{1}{2}}).$

Assume $\tilde{\mu},\tilde{e}$ are RMSE-consistent with $r_{\mu},r_{e}$,
respectively, so that $r_{\mu}+r_{e}\geq\frac{1}{2}$.

######  Theorem 1 (Rate double-robustness).

Consider the one step-estimator with the empirical Gateaux derivative. Under
1, when $\epsilon\lambda^{-d/2}=o(n^{-\max(r_{\mu},r_{e})}),$ and
$\lambda^{\beta}=o(n^{-\frac{1}{2}})$:

$\left(\Psi(\tilde{P})+\frac{1}{n}\sum_{{i}=1}^{n}{{\tilde{\phi}}}(O_{i})\right)-\Psi(P)=O_{p}(n^{-\frac{1}{2}}).$

We now interpret the qualitative estimation implications of the analysis.
Theorem 1 states the rate conditions required of numerical approximation in
$\epsilon,\lambda$ to preserve the statistical property of rate double
robustness. Assuming that the nuisances induced by the probability density
estimations satisfy a product-rate condition that the product of their RMSEs
is faster than the parametric $O(n^{-\frac{1}{2}})$ rate, the form of
Proposition 1 suggests that the rate of RMSE convergence of perturbed
nuisances is required to be faster than the rates of the unperturbed nuisance
functions, $\max(r_{\mu},r_{e})$, which can be a slower rate than the
$\epsilon=o(n^{-\frac{1}{2}})$ implied by the generic analysis of finite
differences. For direct comparison, the rate of Carone et al. (2018, Thm. 5)
suggests the rate requirement of
$\left(\epsilon\lambda^{-d_{1}}\right)^{s}+\lambda^{-2}=o(n^{-\frac{1}{2}}).$

In general, the dimension dependence on $\lambda$ in Lemma 1 will prevent
setting a slower rate for $\epsilon$ for nontrivial dimensions. Note that
without the smoothing requirement, the potential improvement in the rate of
$\epsilon$ could be on the order of generic rate improvements implied by a
central difference scheme. That is, we would have attained the approximation
benefits of generically approximation-improved finite difference schemes via
the specialized statistical structure of the target adjustments.

Our analysis did make use of our characterization in Corollary 1, i.e. a
deterministic approximation error from AIPW with induced nuisance functions.

##  5 Applications: Mean Potential Outcome and Dynamic Settings

In terms of empirics, we first show an plot of numerical approximation error
ranging over $(\epsilon,\lambda)$ in the simple case of the mean potential
outcome as validation. We apply the tools we derived in the previous section
to analyze more complicated settings. We study longitudinal data (dynamic
treatment regime) and the policy linear program for an infinite-horizon Markov
decision with linear constraints, with attention to how different functionals
have different interactions between approximation error and statistical
properties.

(a)

(b)

Figure 1: (a) Epsilon-lambda plot for mean potential outcome: MAE in
approximation of one-step adjustment. (b) Error in estimation of
$\mathbb{E}[Y(1)]$ over $n$.

###  5.1 Mean potential outcome

We first conduct a validation of our decomposition and theoretical
characterization. Our data generating process includes a piecewise linear
outcome model (such that kernel regression is misspecified). We conduct plug-
in evaluation of the mean potential outcome with kernel density estimates. In
Figure 1(a) we include an $(\epsilon,\lambda)$-plot for the simple case of
AIPW (see Carone et al. (2018) for more discussion and examples). We consider
a one-dimensional case with uniformly distributed $X$, piecewise-linear $Y$,
and smooth propensity scores that are logistic in $\sin(X)$. We use $n=500$
and fix the bandwidth $h=0.05.$ We use colors to denote magnitude of the mean
absolute error, included in text on the heatmap. Without loss of generality,
we study the estimation of a mean under missingness, $\mathbb{E}[Y(1)]$.
Figure 1(b) illustrates the estimation error of various strategies with the
comparable kernel-based estimates (DM is regression adjustment): even in this
simple setting AIPW offers benefits and sample efficiency relative to DM. We
include further numerical experiments in the Appendix.

###  5.2 Dynamic treatment regimes

We study the estimation properties of the empirical Gateaux derivative for
multi-stage DTR from the probability density representation. In $T$-stage
dynamic treatment regimes, the causal quantity of interest is the mean
potential outcome $\mathbb{E}\left[Y({\bar{a}})\right]$, where
$\bar{a}=(a_{0},\ldots,a_{T})$ is the (deterministic) treatment strategy.
Assuming $Y({\bar{a}})$ is sequentially ignorable given treatment and
covariate history $(\bar{A}_{t},\bar{X}_{t})$ at each time $t$, this causal
quantity can be identified by the $g$-formula,

$\textstyle\mathbb{E}\left[Y({\bar{a}})\right]=\int\mathbb{E}\left[Y\,|\,\bar{A}=\bar{a},\bar{X}=\bar{x}\right]\prod_{t=1}^{T}\tilde{p}(x_{t}\,|\,\bar{a}_{t-1},\bar{x}_{t-1})\mathop{}\\!\mathrm{d}\bar{x},$

where $\bar{A}=(A_{0},\ldots,A_{t}),\bar{X}=(X_{0},\ldots,X_{t})$. To derive
its influence function, we take the empirical Gateaux derivative by
considering the perturbation in the direction of
$o_{i}=(\bar{x}_{i},\bar{a}_{i},y_{i})$, where
$\bar{x}_{i}=(x_{i0},\ldots,x_{it})$ and $\bar{a}_{i}=(a_{i0},\ldots,a_{it})$,

$\displaystyle\Psi(\tilde{P}^{i}_{\epsilon})$
$\displaystyle\textstyle=\int\left(\int
y\frac{\tilde{p}_{\epsilon}(y,\bar{A}=\bar{a},\bar{X}=\bar{x})}{\tilde{p}_{\epsilon}(\bar{A}=\bar{a},\bar{X}=\bar{x})}\mathop{}\\!\mathrm{d}y\right)\prod_{t=1}^{T}\frac{\tilde{p}_{\epsilon}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}_{\epsilon}(\bar{a}_{t-1},\bar{x}_{t-1}))}\tilde{p}_{\epsilon}(x_{0})\mathop{}\\!\mathrm{d}\bar{x}.$
(4)

Below we characterize how this empirical Gateaux derivative at the smoothed
distribution differs from the one at the (unsmoothed) estimated distribution.
The following result is analogous to Propositions 1 and 1 but is extended to
the dynamic treatment regimes.

######  Proposition 2.

$\displaystyle\epsilon^{-1}({\Psi(\tilde{P}^{i}_{\epsilon})-\Psi(\tilde{P})})$
$\displaystyle=\left(\mathbb{E}_{\tilde{p}(\bar{x}_{1:T})}\left[\mathbb{E}_{\tilde{P}_{\epsilon}}\left[Y\,|\,\bar{A},\bar{X}\right]\mid\bar{x}_{0}\right]-\Psi(\tilde{P})\right)+\frac{\mathbb{I}\left[\bar{A}=\bar{a}_{i}\right]}{\tilde{p}_{\epsilon}(\bar{a}\,|\,\bar{x})}{}\left[Y_{i}-\mathbb{E}_{\tilde{p}(\bar{x})}\left[\mathbb{E}_{\tilde{P}}\left[Y\,|\,\bar{A},\bar{X}\right]\right]\right]$
$\displaystyle\quad+\sum_{s=1}^{T}\frac{\mathbb{I}\left[\bar{A}=\bar{a}_{i,T-s}\right]}{\tilde{p}_{\epsilon}(\bar{a}_{T-s}\,|\,\bar{x}_{T-s})}\left\\{\mathbb{E}_{\tilde{P}}\left[\mathbb{E}_{\tilde{P}_{\epsilon}}\left[Y\,|\,\bar{A},\bar{X}\right]\mid\bar{a}_{i,T-s},\bar{x}_{i,T-s+1}\right]\right.$
(5)
$\displaystyle\quad\quad\quad\quad\left.-\mathbb{E}_{\tilde{P}}\left[\mathbb{E}_{\tilde{P}_{\epsilon}}\left[Y\,|\,\bar{A},\bar{X}\right]\mid\bar{a}_{i,T-s},\bar{x}_{i,T-s}\right]\right\\}+O(\epsilon^{2})+O(\lambda^{\beta})$

The expansion verifies that the requirements for numerical approximation are
similar in $\epsilon$ as in the case of a single-timestep mean potential
outcome, despite the additional complexity of the dynamic treatment regime
function as a nested expectation. Again, smoothing incurs overall unfavorable
dependence in the dimension, as the extension of Lemma 1 to
$\mathbb{E}_{\tilde{P}_{\epsilon}}[Y\mid\bar{A},\bar{X}]$ also incurs
unfavorable dependence on the dimension. The history-dependent nuisances imply
that the dimension additionally grows in time horizon.

###  5.3 Perturbations for (constrained) infinite-horizon off-policy
evaluation

We show how this framework can be applied to derive generic bias corrections
for infinite-horizon off-policy optimization in offline reinforcement
learning, via the canonical linear programming characterization of (Puterman,
2014; De Farias and Van Roy, 2003). Our application is relevant to a recent
line of work on off-policy evaluation and learning (OPE/L) (Jiang and Li,
2016; Thomas and Brunskill, 2016; Bibaut et al., 2019; Kallus and Uehara,
2022; Tang et al., 2019; Kallus and Zhou, 2020, 2022).

##### Setup.

We recall the classical linear programming formulation of finding an optimal
policy in a tabular (finite-state, finite-action) Markov decision process, due
to (Puterman, 2014; De Farias and Van Roy, 2003). Recent work has revisited
this formulation with interest in developing primal-dual algorithms (Wang,
2017; Lee and He, 2018, 2018; Abbasi-Yadkori et al., 2019; Serrano and Neu,
2020). The dynamics follow an infinite-horizon Markov decision process with
discount factor $\gamma$. The offline dataset comprises of trajectories of
state, action, next state observations:
$\\{(s^{0}_{i},a^{0}_{i},\dots,s_{i}^{t},\dots,s_{i}^{T})\\}_{1:N}$. We derive
our bias adjustment under a statistical model where the stationary
distribution factorizes as ${p(s,a,s^{\prime})=d(s,a)P(s\mid s,a)}.$ From the
joint state-action-state occupancy distribution $p(s,a,s^{\prime})$, we
estimate the marginalization $d(s,a)$, the stationary state-action occupancy
distribution, which we use to estimate the transition probability
$P(s^{\prime}\mid a,s)=\frac{p(s,a,s^{\prime})}{d(s,a)}.$ Let $\mu_{0}(s)$
denote the initial state distribution, estimated from offline data and assumed
to be equivalent in the original confounded and evaluation setting. Because we
focus on the tabular case with discrete distributions, we do not require
smoothing.

##### The optimal policy linear program.

The optimal value function solves the following linear program in a finite-
state and finite-action setting; the objective value at optimality is the
policy value. Here, $P_{a}\in\mathbb{R}^{|S|\times|S|}$ is transition matrix
at action $a$, with $P_{a}(s,s^{\prime})=P(s^{\prime}\mid s,a)$, $\mu_{0}(s)$
is the marginal occupancy distribution of the initial state. We assume, as is
fairly common in theoretical analyses of reinforcement learning, that the
reward function is a known function of state and action, $r(s,a)$. Let
${r}_{a}=\\{r(s,a)\\}_{s\in\mathcal{S}}$ be the vector of reward values per
state, for fixed action $a$.

$\textstyle\Psi_{D}(P)=\underset{V}{\operatorname{min}}\\{(1-\gamma){\mu_{0}}^{\top}{V}\colon\left(I-\gamma
P_{a}\right){V}-{r}_{a}\geq 0,\quad\forall a\in\mathcal{A}\\}$ (6)

The dual formulation is well-known to parametrize the stationary occupancy
probabilities of the optimal policy (at optimality):

$\textstyle\Psi_{P}(P)=\underset{\mu}{\operatorname{max}}\left\\{\sum_{a\in\mathcal{A}}\mu_{a}^{\top}{r}_{{a}}\colon\sum_{a\in\mathcal{A}}(I-\gamma
P_{a}^{\top})\mu_{a}=(1-\gamma){\mu_{0}},\quad\mu_{a}\geq 0,\quad\forall
a\in\mathcal{A}\right\\}$ (7)

We verify the bias correction from the Gateaux derivative
$\frac{\mathrm{d}}{\mathrm{d}\epsilon}\Psi_{D}(P_{\epsilon})\Big{|}_{\epsilon=0}$.
Note that evaluating this functional at the perturbed distribution corresponds
to evaluating the linear program with perturbed left-hand side coefficients,
i.e. perturbing the constraints of the coefficient matrix as well as the
objective (marginalizing over a perturbed initial state distribution. We
evaluate the Gateaux derivative in the direction of
$(\tilde{s},\tilde{a},\tilde{s}^{\prime}).$ We assume nondegeneracy (unique
optimal solution); this can be relaxed by using more sophisticated
perturbation analysis.

######  Assumption 2.

The optimal solution is unique; there is no degeneracy.

######  Proposition 3.

Perturb in the direction of a generic observation
$(\tilde{s},\tilde{a},\tilde{s}^{\prime}).$

$\frac{\mathrm{d}}{\mathrm{d}\epsilon}\Psi_{D}(P_{\epsilon})\Big{|}_{\epsilon=0}=(1-\gamma)V^{*}(\tilde{s})-\frac{\mu^{*}(\tilde{s},\tilde{a})}{d(\tilde{s},\tilde{a})}\left(r(\tilde{s},\tilde{a})+\gamma
V^{*}(\tilde{s}^{\prime})-V^{*}(\tilde{s})\right)-\Psi_{D}(P)$ (8)

##### Sketch of argument.

The statement of Proposition 3 is not surprising per se, since we can
immediately verify the mean-zero sample version of the Bellman residual, and
so one might plausibly derive this form from double robustness. Rather, we
provide an alternative argument that applies perturbation analysis of
optimization programs with the specific perturbations from the finite
difference. We build on the sensitivity analysis characterization of Freund
(1985) which studies perturbations of linear programs with respect to
constraint matrix perturbations. For $\epsilon$ small enough, the active basis
remains the same. When the perturbation matrix is $P+\epsilon G$ for some
matrix $G$, Freund (1985) notes that the derivative of the optimal value is
$d^{\prime}(\epsilon)=-\overline{\pi}G\overline{x}$ (where
$(\overline{\pi},\overline{x})$ are the dual and primal optimal solutions,
respectively) can be obtained from the derivative of a matrix inverse, and
obtains higher order expansions. We generalize this argument for the specific
form of the perturbed transition probabilities.

We next compute the approximation error that arises from evaluating finite
differences, so long as $\epsilon$ is small enough to maintain the same active
basis, which is empirically verifiable.

######  Proposition 4 (Error analysis from finite-differencing).

Perturb in the direction of a generic observation
$(\tilde{s},\tilde{a},\tilde{s}^{\prime}).$ Then:

$\epsilon^{-1}(\Psi(P_{\epsilon})-\Psi(P))=(1-\gamma)V^{*}_{\epsilon}(\tilde{s})-\frac{\mu^{*}(\tilde{s},\tilde{a})}{d_{\epsilon}(\tilde{s},\tilde{a})}\left(r(\tilde{s},\tilde{a})+\gamma
V(\tilde{s}^{\prime})-V_{\epsilon}(\tilde{s})\right)-\Psi_{D}(P)+O(\epsilon)$
(9)

In this setting, unlike the earlier case of the treatment mean in Corollary 1
and Theorem 1, even in a favorable discrete-state case without smoothing
requirements, it does not appear that rate-double robustness would admit
weaker numerical requirements on $\epsilon$. Ultimately, this arises because
finite-difference approximation in this case involves perturbed matrix
inversion, rather than composition of sum and product rules of differentiation
as in the case of the previous functional. In the analysis, Equation 47 shows
that $V_{\epsilon}(s)-V(s)=O(\epsilon).$ so that Proposition 4 holds with the
unperturbed nuisances and the same order of approximation error.

##### Interpreting the finite-difference approach.

We take a moment to translate the benefits of a finite-difference approach to
estimating a one-step adjustment in the context of reinforcement learning. In
reinforcement learning, a common paradigm is that of “model-based”
reinforcement learning, i.e. learning a dynamics model (estimating transition
distribution) and plugging it into a planning problem or the linear
programming formulation. A benefit of finite-differences is that one can
approximate doubly-robust adjustments by additional black-box function
evaluations that they would already be computing, without estimating new
nuisance functions.

##### Generalization to constraints.

We argue the main benefit of this approach is that by computing these finite
differences, we can obtain generic bias corrections via a more broadly
applicable argument that does not require interpretation of the closed-form
solution to an optimization problem. For example, consider a relevant subclass
of constrained Markov decision processes comprising of additional constraints
on the policy variables: the linear programming formulation is particularly
appealing because linear constraints can be directly added (Altman, 1999).
Another benefit of the finite-differences approach is that one-step estimation
can be obtained by evaluating the functional at feasible points. In our case,
this corresponds to optimizing with respect to feasible probability
distributions.

######  Example 3 (Constrained policy optimization).

$\textstyle\Psi_{P}(P)=\operatorname{max}\left\\{\sum_{a\in\mathcal{A}}\mu_{a}^{\top}{r}_{{a}}\colon\sum_{a\in\mathcal{A}}(I-\gamma
P_{a}^{\top})\mu_{a}=(1-\gamma){\mu_{0}},\;\;\mu_{a}\geq 0,\;\;\forall
a\in\mathcal{A},\;\mu\in\mathcal{P}\right\\}$

The additional constraint $\mu\in\mathcal{P}$ reflects linearly representable
constraints on state-occupancy under the optimal policy. For example, relevant
constraints could require the occupancy in certain risky states to be bounded.
Paternain et al. (2019) develops a framework with a “caution function”, a
convex function on the state-action occupancy measures; this analysis applies
to linear caution functions.

To compare to related work, we note that double robustness has been studied
(for policy evaluation) in Tang et al. (2019). Our derivation via perturbation
analysis is similar to Tang et al. (2019, Thm. 4.1)’s illustration via the
Lagrangian. This is unsurprising since the Lagrangian provides the dual
characterization that also underlies perturbation analysis. Relative to the
doubly robust estimators derived elsewhere, this influence function may be
higher variance.

##  6 Conclusion

##### Future directions.

An important direction for future work remains empirical comparisons of the
performance of estimators based on the numerical differentiation approach. (In
our experiences, we found some variation in results based on various
configurations of numerical integration). The smoothing perturbation, although
it restores absolute continuity, introduces technical barriers in the implied
approximation rates with regards to dimension-dependence. Other arguments that
reduce the dependence on smoothing would be of interest. It would be
interesting to generalize some comparisons of finite vs. analytical
differences further, such as to Z- and M-estimators. Our analysis also relied
on a priori knowledge of beneficial statistical properties: it would be
interesting to verify these adaptively.

##### Concluding remarks.

Our analysis suggests that the current major technical limitations include the
dimension-dependence due to smoothing. Regarding Q2, our analysis of different
functionals surfaced fine-grained differences in opportunities for statistical
properties to weaken requirements of computational approximation: Theorem 1
and Proposition 4 have different, estimand-dependent implications.

##### Acknowledgements.

Angela Zhou gratefully acknowledges support from the Foundations of Data
Science Institute and the Simons Institute’s Program on Causality. Part of
this work was done while the author was visiting the Simons Institute for the
Theory of Computing. Yixin Wang acknowledges National Science Foundation Grant
NSF-CHE-2231174.

## References

  * Abbasi-Yadkori et al. [2019] Y. Abbasi-Yadkori, P. L. Bartlett, X. Chen, and A. Malek.  Large-scale markov decision problems via the linear programming dual.  _arXiv preprint arXiv:1901.01992_ , 2019.
  * Alaa and Van Der Schaar [2019] A. Alaa and M. Van Der Schaar.  Validating causal inference models via influence functions.  In _International Conference on Machine Learning_ , pages 191–201. PMLR, 2019.
  * Altman [1999] E. Altman.  _Constrained Markov decision processes: stochastic modeling_.  Routledge, 1999.
  * Bhattacharya et al. [2020] R. Bhattacharya, R. Nabi, and I. Shpitser.  Semiparametric inference for causal effects in graphical models with hidden variables.  _arXiv preprint arXiv:2003.12659_ , 2020.
  * Bibaut et al. [2019] A. F. Bibaut, I. Malenica, N. Vlassis, and M. J. van der Laan.  More efficient off-policy evaluation through regularized targeted learning.  _arXiv preprint arXiv:1912.06292_ , 2019.
  * Bickel and Ritov [2003] P. J. Bickel and Y. Ritov.  Nonparametric estimators which can be” plugged-in”.  _The Annals of Statistics_ , 31(4):1033–1053, 2003.
  * Bickel et al. [1993] P. J. Bickel, C. A. Klaassen, P. J. Bickel, Y. Ritov, J. Klaassen, J. A. Wellner, and Y. Ritov.  _Efficient and adaptive estimation for semiparametric models_ , volume 4.  Johns Hopkins University Press Baltimore, 1993.
  * Bonnans and Shapiro [1998] J. F. Bonnans and A. Shapiro.  Optimization problems with perturbations: A guided tour.  _SIAM review_ , 40(2):228–264, 1998.
  * Bravo et al. [2020] F. Bravo, J. C. Escanciano, and I. Van Keilegom.  Two-step semiparametric empirical likelihood inference.  _The Annals of Statistics_ , 48(1):1–26, 2020\.
  * Carone et al. [2018] M. Carone, A. R. Luedtke, and M. J. van der Laan.  Toward computerized efficient estimation in infinite-dimensional models.  _Journal of the American Statistical Association_ , 2018.
  * Chakraborty and Moodie [2013] B. Chakraborty and E. Moodie.  Statistical methods for dynamic treatment regimes.  _Springer-Verlag. doi_ , 10:978–1, 2013.
  * Chernozhukov et al. [2016] V. Chernozhukov, J. C. Escanciano, H. Ichimura, W. K. Newey, and J. M. Robins.  Locally robust semiparametric estimation.  _arXiv preprint arXiv:1608.00033_ , 2016.
  * Chernozhukov et al. [2018] V. Chernozhukov, D. Chetverikov, M. Demirer, E. Duflo, C. Hansen, W. Newey, and J. Robins.  Double/debiased machine learning for treatment and structural parameters, 2018.
  * Chernozhukov et al. [2020] V. Chernozhukov, W. Newey, R. Singh, and V. Syrgkanis.  Adversarial estimation of riesz representers.  _arXiv preprint arXiv:2101.00009_ , 2020.
  * Chernozhukov et al. [2021] V. Chernozhukov, W. K. Newey, V. Quintas-Martinez, and V. Syrgkanis.  Automatic debiased machine learning via neural nets for generalized linear regression.  _arXiv preprint arXiv:2104.14737_ , 2021.
  * Chernozhukov et al. [2022] V. Chernozhukov, W. Newey, V. M. Quintas-Martinez, and V. Syrgkanis.  Riesznet and forestriesz: Automatic debiased machine learning with neural nets and random forests.  In _International Conference on Machine Learning_ , pages 3901–3914. PMLR, 2022.
  * De Farias and Van Roy [2003] D. P. De Farias and B. Van Roy.  The linear programming approach to approximate dynamic programming.  _Operations research_ , 51(6):850–865, 2003.
  * Duchi and Ruan [2021] J. C. Duchi and F. Ruan.  Asymptotic optimality in stochastic optimization.  _The Annals of Statistics_ , 49(1):21–48, 2021\.
  * Duchi et al. [2015] J. C. Duchi, M. I. Jordan, M. J. Wainwright, and A. Wibisono.  Optimal rates for zero-order convex optimization: The power of two function evaluations.  _IEEE Transactions on Information Theory_ , 61(5):2788–2806, 2015.
  * Farrell et al. [2021] M. H. Farrell, T. Liang, and S. Misra.  Deep learning for individual heterogeneity: an automatic inference framework.  2021\.
  * Fisher and Kennedy [2021] A. Fisher and E. H. Kennedy.  Visually communicating and teaching intuition for influence functions.  _The American Statistician_ , 75(2):162–172, 2021\.
  * Flaxman et al. [2004] A. D. Flaxman, A. T. Kalai, and H. B. McMahan.  Online convex optimization in the bandit setting: gradient descent without a gradient.  _arXiv preprint cs/0408007_ , 2004.
  * Fornberg [1981] B. Fornberg.  Numerical differentiation of analytic functions.  _ACM Transactions on Mathematical Software (TOMS)_ , 7(4):512–526, 1981.
  * Frangakis et al. [2015] C. E. Frangakis, T. Qian, Z. Wu, and I. Diaz.  Deductive derivation and turing-computerization of semiparametric efficient estimation.  _Biometrics_ , 71(4):867–874, 2015.
  * Freund [1985] R. M. Freund.  Postoptimal analysis of a linear program under simultaneous changes in matrix coefficients.  In _Mathematical Programming Essays in Honor of George B. Dantzig Part I_ , pages 1–13. Springer, 1985.
  * Giordano et al. [2019a] R. Giordano, M. I. Jordan, and T. Broderick.  A higher-order swiss army infinitesimal jackknife.  _arXiv preprint arXiv:1907.12116_ , 2019a.
  * Giordano et al. [2019b] R. Giordano, W. Stephenson, R. Liu, M. Jordan, and T. Broderick.  A swiss army infinitesimal jackknife.  In _The 22nd International Conference on Artificial Intelligence and Statistics_ , pages 1139–1147. PMLR, 2019b.
  * Grant and Boyd [2008] M. C. Grant and S. P. Boyd.  Graph implementations for nonsmooth convex programs.  In _Recent advances in learning and control_ , pages 95–110. Springer, 2008.
  * Guo et al. [2022] W. Guo, M. I. Jordan, and A. Zhou.  Off-policy evaluation with policy-dependent optimization response.  _arXiv preprint arXiv:2202.12958_ , 2022.
  * Gupta et al. [2021] V. Gupta, M. Huang, and P. Rusmevichientong.  Debiasing in-sample policy performance for small-data, large-scale optimization.  _arXiv preprint arXiv:2107.12438_ , 2021.
  * Hahn [1998] J. Hahn.  On the role of the propensity score in efficient semiparametric estimation of average treatment effects.  _Econometrica_ , pages 315–331, 1998.
  * Hampel [1974] F. R. Hampel.  The influence curve and its role in robust estimation.  _Journal of the american statistical association_ , 69(346):383–393, 1974.
  * Hannah et al. [2011] L. A. Hannah, D. M. Blei, and W. B. Powell.  Dirichlet process mixtures of generalized linear models.  _Journal of Machine Learning Research_ , 12(6), 2011.
  * Hausman and Newey [1995] J. A. Hausman and W. K. Newey.  Nonparametric estimation of exact consumers surplus and deadweight loss.  _Econometrica: Journal of the Econometric Society_ , pages 1445–1476, 1995.
  * Hines et al. [2022] O. Hines, O. Dukes, K. Diaz-Ordaz, and S. Vansteelandt.  Demystifying statistical learning based on efficient influence functions.  _The American Statistician_ , pages 1–13, 2022.
  * Hong and Li [2018] H. Hong and J. Li.  The numerical delta method.  _Journal of Econometrics_ , 206(2):379–394, 2018\.
  * Hong et al. [2015] H. Hong, A. Mahajan, and D. Nekipelov.  Extremum estimation and numerical derivatives.  _Journal of Econometrics_ , 188(1):250–263, 2015\.
  * Huber [2004] P. J. Huber.  _Robust statistics_ , volume 523.  John Wiley & Sons, 2004.
  * Ichimura and Newey [2015] H. Ichimura and W. K. Newey.  The influence function of semiparametric estimators.  _arXiv preprint arXiv:1508.01378_ , 2015.
  * Ichimura and Newey [2022] H. Ichimura and W. K. Newey.  The influence function of semiparametric estimators.  _Quantitative Economics_ , 13(1):29–61, 2022\.
  * Ito et al. [2018] S. Ito, A. Yabe, and R. Fujimaki.  Unbiased objective estimation in predictive optimization.  In _International Conference on Machine Learning_ , pages 2176–2185. PMLR, 2018.
  * Jesson et al. [2021] A. Jesson, S. Mindermann, Y. Gal, and U. Shalit.  Quantifying ignorance in individual-level causal-effect estimates under hidden confounding.  In _International Conference on Machine Learning_ , pages 4829–4838. PMLR, 2021.
  * Jiang and Li [2016] N. Jiang and L. Li.  Doubly robust off-policy value evaluation for reinforcement learning.  _Proceedings of the 33rd International Conference on Machine Learning_ , 2016.
  * Jung et al. [2021a] Y. Jung, J. Tian, and E. Bareinboim.  Double machine learning density estimation for local treatment effects with instruments.  _Advances in Neural Information Processing Systems_ , 34:21821–21833, 2021a.
  * Jung et al. [2021b] Y. Jung, J. Tian, and E. Bareinboim.  Estimating identifiable causal effects through double machine learning.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 35, pages 12113–12122, 2021b.
  * Kallus and Uehara [2020] N. Kallus and M. Uehara.  Double reinforcement learning for efficient off-policy evaluation in markov decision processes.  _J. Mach. Learn. Res._ , 21(167):1–63, 2020.
  * Kallus and Uehara [2022] N. Kallus and M. Uehara.  Efficiently breaking the curse of horizon in off-policy evaluation with double reinforcement learning.  _Operations Research_ , 2022.
  * Kallus and Zhou [2020] N. Kallus and A. Zhou.  Confounding-robust policy evaluation in infinite-horizon reinforcement learning.  _Advances in Neural Information Processing Systems_ , 33:22293–22304, 2020.
  * Kallus and Zhou [2022] N. Kallus and A. Zhou.  Stateful offline contextual policy evaluation and learning.  In _International Conference on Artificial Intelligence and Statistics_ , pages 11169–11194. PMLR, 2022.
  * Kallus et al. [2019] N. Kallus, X. Mao, and M. Uehara.  Localized debiased machine learning: Efficient inference on quantile treatment effects and beyond.  _arXiv preprint arXiv:1912.12945_ , 2019.
  * Kandasamy et al. [2015] K. Kandasamy, A. Krishnamurthy, B. Poczos, L. Wasserman, et al.  Nonparametric von mises estimators for entropies, divergences and mutual informations.  _Advances in Neural Information Processing Systems_ , 28, 2015.
  * Kennedy [2016] E. H. Kennedy.  Semiparametric theory and empirical processes in causal inference.  In _Statistical causal inferences and their applications in public health research_ , pages 141–167. Springer, 2016.
  * Kennedy [2022] E. H. Kennedy.  Semiparametric doubly robust targeted double machine learning: a review.  _arXiv preprint arXiv:2203.06469_ , 2022.
  * Koh and Liang [2017] P. W. Koh and P. Liang.  Understanding black-box predictions via influence functions.  In _International conference on machine learning_ , pages 1885–1894. PMLR, 2017.
  * Lee and He [2018] D. Lee and N. He.  Stochastic primal-dual q-learning.  _arXiv preprint arXiv:1810.08298_ , 2018.
  * Newey [1994a] W. K. Newey.  The asymptotic variance of semiparametric estimators.  _Econometrica: Journal of the Econometric Society_ , pages 1349–1382, 1994a.
  * Newey [1994b] W. K. Newey.  Kernel estimation of partial means and a general variance estimator.  _Econometric Theory_ , 10(2):1–21, 1994b.
  * Oktay et al. [2020] D. Oktay, N. McGreivy, J. Aduol, A. Beatson, and R. P. Adams.  Randomized automatic differentiation.  _arXiv preprint arXiv:2007.10412_ , 2020.
  * Paternain et al. [2019] S. Paternain, L. Chamon, M. Calvo-Fullana, and A. Ribeiro.  Constrained reinforcement learning has zero duality gap.  _Advances in Neural Information Processing Systems_ , 32, 2019.
  * Pfanzagl [1990] J. Pfanzagl.  Estimation in semiparametric models.  In _Estimation in Semiparametric Models_ , pages 17–22. Springer, 1990.
  * Puterman [2014] M. L. Puterman.  _Markov decision processes: discrete stochastic dynamic programming_.  John Wiley & Sons, 2014.
  * Robins et al. [1994] J. M. Robins, A. Rotnitzky, and L. P. Zhao.  Estimation of regression coefficients when some regressors are not always observed.  _Journal of the American statistical Association_ , 89(427):846–866, 1994.
  * Rotnitzky et al. [2021] A. Rotnitzky, E. Smucler, and J. M. Robins.  Characterization of parameters with a mixed bias property.  _Biometrika_ , 108(1):231–238, 2021.
  * Sauer [2012] T. Sauer.  Numerical analysis, 2012.
  * Serrano and Neu [2020] J. B. Serrano and G. Neu.  Faster saddle-point optimization for solving large-scale markov decision processes.  In _Learning for Dynamics and Control_ , pages 413–423. PMLR, 2020\.
  * Tang et al. [2019] Z. Tang, Y. Feng, L. Li, D. Zhou, and Q. Liu.  Doubly robust bias reduction in infinite horizon off-policy estimation.  _arXiv preprint arXiv:1910.07186_ , 2019.
  * Thomas and Brunskill [2016] P. Thomas and E. Brunskill.  Data-efficient off-policy policy evaluation for reinforcement learning.  _Journal of Machine Learning Research_ , 2016.
  * Tsiatis [2006] A. A. Tsiatis.  Semiparametric theory and missing data.  2006\.
  * Ullah and Pagan [1999] A. Ullah and A. Pagan.  _Nonparametric econometrics_.  Cambridge university press Cambridge, 1999.
  * Van der Laan and Robins [2003] M. J. Van der Laan and J. M. Robins.  _Unified methods for censored longitudinal data and causality_ , volume 5.  Springer, 2003.
  * van der Laan et al. [2018] M. J. van der Laan, A. Bibaut, and A. R. Luedtke.  Cv-tmle for nonpathwise differentiable target parameters.  In _Targeted Learning in Data Science_ , pages 455–481. Springer, 2018.
  * Van der Vaart [2000] A. W. Van der Vaart.  _Asymptotic statistics_ , volume 3.  Cambridge university press, 2000.
  * Wang [2017] M. Wang.  Randomized linear programming solves the discounted markov decision problem in nearly-linear running time.  2017\.
  * Wilson et al. [2020] A. Wilson, M. Kasy, and L. Mackey.  Approximate cross-validation: Guarantees for model assessment and selection.  In _International Conference on Artificial Intelligence and Statistics_ , pages 4530–4540. PMLR, 2020.

Appendix: Empirical Gateaux Derivatives for Causal Inference

##  Appendix A Overview and background for influence functions

For a broad overview of influence functions and their role in causal
inference, see Tsiatis [2006], Hines et al. [2022], Kennedy [2022].

We focus on functionals that admit an asymptotically linear representation:

$\displaystyle\Psi(\tilde{P})-\Psi(P)$
$\displaystyle=-\int\phi(u;\tilde{P})dP(u)+R_{2}\left(\tilde{P},P\right),$

where the influence function, $\phi(u;\tilde{P})$, captures the first-order
bias in the functional. That is, we focus on functionals that are smooth
enough to admit a von Mises expansion with a mean-zero, finite-variance
influence function and a second-order remainder term
$R_{2}\left(P_{1},P\right)$.

##  Appendix B Proofs for Section 4

####  B.0.1 Supplementary lemmas

######  Lemma 2 (Intermediate decomposition).

$\displaystyle\frac{1}{\epsilon}\int\int
y\left(\frac{\tilde{p}_{\epsilon}(y,A=1,x)}{\tilde{p}_{\epsilon}(A=1,x)}-\frac{\tilde{p}(y,A=1,x)}{\tilde{p}(A=1,x)}\right)\tilde{p}(x)\mathrm{d}y\;\mathrm{d}x$
$\displaystyle\qquad=\int\frac{\tilde{p}(x)}{\tilde{p}_{\epsilon}(A=1,x)}\mathbb{I}\left[a_{i}=1\right]\tilde{\delta}_{(x_{i})}^{\lambda}(x)\left\\{\left(\int
y\;\tilde{\delta}_{y_{i}}^{\lambda}(y)\mathrm{d}y\right)-\mathbb{E}_{\tilde{p}}[Y\mid
A=1,x]\right\\}\mathrm{d}x.$

####  B.0.2 Proofs

###### Proof of Proposition 1.

We have:

$\displaystyle\frac{\Psi(\tilde{P}_{\epsilon_{i}})-\Psi(\tilde{P})}{\epsilon}$
$\displaystyle=\frac{1}{\epsilon}\int\left\\{\int
y\left[\frac{\tilde{p}_{\epsilon}(y,A=1,x)}{\tilde{p}_{\epsilon}(A=1,{x})}\left((1-\epsilon)\tilde{p}({x})+\epsilon\tilde{\delta}_{x_{{i}}}^{\lambda}\right)-\frac{\tilde{p}(y,A=1,x)}{\tilde{p}(A=1,x)}\tilde{p}(x)\right]\mathrm{d}y\right\\}\mathrm{d}x$
$\displaystyle=\left[\int\tilde{p}(x)\int
y\left(\frac{\tilde{p}_{\epsilon}(y,A=1,x)}{\tilde{p}_{\epsilon}(A=1,{x})}-\frac{\tilde{p}(y,A=1,x)}{\tilde{p}(A=1,x)}\right)\;\mathrm{d}y\;\mathrm{d}x\right]$
$\displaystyle\qquad+\cancel{\frac{1}{\epsilon}}\cancel{\epsilon}\int\int
y\left[\frac{\tilde{p}_{\epsilon}(y,A=1,x)}{\tilde{p}_{\epsilon}(A=1,{x})}\tilde{\delta}_{x_{{i}}}^{\lambda}-\frac{\tilde{p}(y,A=1,x)}{\tilde{p}(A=1,x)}\tilde{p}(x)\right]\mathrm{d}y\;\mathrm{d}x.$

The first term can be expanded using the decomposition of Lemma 2. We
recognize the second term as
$\left(\tilde{\mathbb{E}}_{\tilde{P}_{\epsilon}}[Y\mid
A=1,X=x_{{i}}]-\Psi(\tilde{P})\right)$ by definition of the perturbed
nuisances. We obtain:

$\displaystyle\frac{\Psi(\tilde{P}_{\epsilon_{i}})-\Psi(\tilde{P})}{\epsilon}$
$\displaystyle=(1-\epsilon)\int\tilde{p}(x)\int
y\left(\frac{\tilde{p}_{\epsilon}(y,A=1,x)}{\tilde{p}_{\epsilon}(A=1,x)}-\frac{\tilde{p}(y,A=1,x)}{\tilde{p}(A=1,x)}\right)\mathrm{d}y\;\mathrm{d}x$
$\displaystyle\qquad=\int\frac{\tilde{p}(x)}{\tilde{p}_{\epsilon}(A=1,x)}\mathbb{I}\left[a_{i}=1\right]\tilde{\delta}_{x_{i}}^{\lambda}(x)\left\\{\left(\int
y\;\tilde{\delta}_{y_{i}}^{\lambda}(y)\mathrm{d}y\right)-\mathbb{E}_{\tilde{P}}[Y\mid
A=1,x]\right\\}\;\mathrm{d}x$
$\displaystyle\qquad+\left(\tilde{\mathbb{E}}_{\tilde{P}_{\epsilon}}[Y\mid
A=1,X=x_{{i}}]-\Psi(\tilde{P})\right).$

∎

###### Proof of Lemma 2.

Verify by direct computation.

For simplicity we consider a single summand in the sum (the claim follows by
applying the same argument to every summand). Recall that
${\tilde{P}_{\epsilon}(A=1,X_{j})=(1-\epsilon)\tilde{p}(A=1,X_{j})+\epsilon\tilde{\mathbb{I}}\left[(o_{\backslash
y})_{\tilde{{i}}}\right]}$

$\displaystyle=\int\frac{\tilde{p}(x)}{\epsilon}\int
y\left[\frac{(1-\epsilon)\tilde{p}(y,A=1,x)+\epsilon\tilde{\delta}_{o_{i}}^{\lambda}}{\tilde{p}_{\epsilon}(A=1,x)}-\frac{\tilde{p}(y,A=1,x)}{\tilde{p}(A=1,x)}\right]\mathrm{d}y\;\mathrm{d}x$
$\displaystyle=\int\frac{\tilde{p}(x)}{\epsilon}\int
y\left[\frac{\left\\{\cancel{(1-\epsilon)}\tilde{p}(y,A=1,x)+\epsilon\tilde{\delta}_{o_{i}}^{\lambda}\right\\}\tilde{p}(A=1,x)}{\tilde{p}_{\epsilon}(A=1,x)\tilde{p}(A=1,x)}\right.$
$\displaystyle\qquad\qquad\qquad\qquad\left.-\frac{\tilde{p}(y,A=1,x)\left\\{\cancel{(1-\epsilon)}\tilde{p}(A=1,x)+\epsilon\tilde{\delta}_{(a_{{i}},x_{{i}})}^{\lambda}\right\\}}{\tilde{p}(A=1,x)\tilde{p}_{\epsilon}(A=1,x)}\right]\mathrm{d}y\;\mathrm{d}x$
$\displaystyle=\int\frac{\tilde{p}(x)}{\cancel{\epsilon}}\int
y\left[\cancel{\epsilon}\frac{\tilde{\delta}_{o_{i}}^{\lambda}}{\tilde{p}_{\epsilon}(A=1,x)}-\cancel{\epsilon}\frac{\tilde{p}(y,A=1,x)\tilde{\delta}_{(a_{{i}},x_{{i}})}^{\lambda}}{\tilde{p}(A=1,x)\tilde{p}_{\epsilon}(A=1,x)}\right]\mathrm{d}y\;\mathrm{d}x$
$\displaystyle=\int\frac{\tilde{p}(x)}{\tilde{p}_{\epsilon}(A=1,x)}\tilde{\delta}_{(a_{{i}},x_{{i}})}^{\lambda}\;\left\\{\int
y{\tilde{\delta}_{y_{i}}^{\lambda}}\mathrm{d}y-\mathbb{E}_{\tilde{p}}[Y\mid
A=1,x]\right\\}\mathrm{d}x$
$\displaystyle=\int\frac{\tilde{p}(x)}{\tilde{p}_{\epsilon}(A=1,x)}\mathbb{I}\left[a_{i}=1\right]\tilde{\delta}_{(x_{i})}^{\lambda}(x)\;\left\\{\int
y{\tilde{\delta}_{y_{i}}^{\lambda}}\mathrm{d}y-\mathbb{E}_{\tilde{p}}[Y\mid
A=1,x]\right\\}\mathrm{d}x$

∎

###### Proof of Lemma 1.

To bound $\mathbb{E}[(\tilde{e}_{\epsilon}(X)-e(X))^{2}],$ argue similarly and
observe that since $\tilde{e}_{\epsilon}=\tilde{p}_{\epsilon}(A=1\mid
x)=\tilde{p}_{\epsilon}(A=1,{x})/\tilde{p}(x),$

$\displaystyle\mathbb{E}[(\tilde{e}_{\epsilon}(X)-\tilde{e}(X))^{2}]$
$\displaystyle=\int(((1-\epsilon)\tilde{p}(A=1,x)+\epsilon\tilde{\delta}_{a,\tilde{x}}^{\lambda}(x))-\tilde{p}(A=1,x))^{2}/\tilde{p}(x)^{2}\mathop{}\\!\mathrm{d}x\cdot\tilde{p}(x)$
$\displaystyle=\epsilon^{2}\int(\tilde{\delta}_{a,\tilde{x}}^{\lambda}(x)-\tilde{p}(A=1,x))^{2}/\tilde{p}(x)\mathop{}\\!\mathrm{d}x$
$\displaystyle=\epsilon^{2}\int((\tilde{\delta}_{a,\tilde{x}}^{\lambda}(x)^{2}-2\tilde{\delta}_{a,\tilde{x}}^{\lambda}(x)\tilde{p}(A=1,x)+\tilde{p}(A=1,x)^{2})/\tilde{p}(x)\mathop{}\\!\mathrm{d}x$
$\displaystyle\leq\epsilon^{2}\int\nu(\tilde{\delta}_{a,\tilde{x}}^{\lambda}(x)^{2}+1)\mathop{}\\!\mathrm{d}x$
$\displaystyle=\epsilon^{2}(\nu\lambda^{-d}R(\tilde{\delta}^{\lambda})+\nu)=O(\epsilon^{2}\lambda^{-d})$

In the second-to-last inequality we lower bound $\tilde{p}(x)\geq\nu$ and
upper bound $\tilde{p}(A=1,x)\leq 1.$ In the last inequality we apply a change
of variables to identify the term that arises in analysis of variance of
kernel density estimates. We incur the $\lambda^{-d}$ dependence since we
assumed that multivariate kernels are product kernels. The roughness term
$R(\tilde{\delta}^{\lambda})=\int\tilde{\delta}^{2}{}(u)du$ depends on the
particular $\tilde{\delta}$ kernel function and is a kernel-dependent
constant. Hence the term scales overall as $\lambda^{-d}$ as in standard
analyses of variance of kernel density estimates; see e.g. Ullah and Pagan
[1999, pg. 23-26].

Next observe that:

$\displaystyle\lVert\tilde{\mu}_{\epsilon}(X)-\mu(X)\rVert$
$\displaystyle\leq\lVert\tilde{\mu}_{\epsilon}(X)-\tilde{\mu}(X)\rVert+\lVert\tilde{\mu}(X)-\mu(X)\rVert$
$\displaystyle=\textstyle\left(\int\left\\{\int\int
y\left(\frac{(1-\epsilon)\tilde{p}(y,A=1,x^{\prime})+\epsilon\tilde{\delta}_{o}^{\lambda}}{(1-\epsilon)\tilde{p}(A=1,x^{\prime})+\epsilon\tilde{\delta}_{a,x}^{\lambda}}-\frac{\tilde{p}(y,A=1,x^{\prime})}{\tilde{p}(A=1,x^{\prime})}\right)\mathrm{d}y\mathrm{d}x^{\prime}\right\\}^{2}p(x)\mathrm{d}x\right)^{\frac{1}{2}}+\lVert\tilde{\mu}(X)-\mu(X)\rVert$
$\displaystyle=\textstyle\left(\int\epsilon^{2}\left\\{\frac{\tilde{p}(x)}{\tilde{p}_{\epsilon}(A=1,x)}\mathbb{I}\left[a_{i}=1\right]\tilde{\delta}_{(x_{i})}^{\lambda}(x)\left\\{y_{i}-\mathbb{E}_{\tilde{p}}[Y\mid
A=1,x]\right\\}\right\\}^{2}p(x)\mathrm{d}x\right)^{\frac{1}{2}}+\lVert\tilde{\mu}(X)-\mu(X)\rVert$
$\displaystyle=O(\epsilon\lambda^{-d/2})+\lVert\tilde{\mu}(X)-\mu(X)\rVert$

where the first inequality comes from the triangle inequality, and the second
equality from similar arguments as in the previous bound. ∎

###### Proof of Theorem 1..

Let

$\tilde{e}(X)=\tilde{p}(A=1,X)/\tilde{p}(X)=\tilde{p}(A=1\mid
X),\tilde{\mu}(X)$

denote the propensity and conditional outcome nuisances arising from plug-in
estimation of $\tilde{P}$.

Let

${\tilde{e}_{\epsilon}(X)=\tilde{p}_{\epsilon}(A=1\mid
X)=\tilde{p}_{\epsilon}(A=1,X)}/{\tilde{p}(X)}$

denote the nuisance arising from plug-in estimation at the perturbed
distribution.

First observe by Corollary 1:

$\displaystyle\frac{\Psi(\tilde{P}_{\epsilon}^{{i}})-\Psi(\tilde{P})}{\epsilon}$
$\displaystyle\textstyle=(1-\epsilon)\frac{\mathbb{I}\left[A_{i}=1\right]}{\tilde{p}_{\epsilon}(A=1\mid
X_{i})}\left(Y_{i}-{\mathbb{E}_{\tilde{P}}[Y\mid
A=1,X_{i}]}\right)+\left({\mathbb{E}}_{\tilde{P}_{\epsilon}}[Y\mid
A=1,X_{{i}}]-\Psi(\tilde{P})\right)+O(\lambda^{J})$

This decomposition studies convergence of an AIPW estimator with possibly
biased but asymptotically consistent nuisances
$\tilde{\mu},\tilde{e}_{\epsilon}$ to the AIPW estimator with oracle nuisances
– a standard analysis – and an additional term that vanishes on the order of
$\epsilon$ (assuming strong overlap and boundedness). Define $\Gamma^{*}(O)$
as the doubly-robust score for AIPW with oracle nuisances $e^{*},\mu^{*}$, and
let $\Gamma^{*}_{j}$ denote its empirical counterpart at
$O_{j}=(X_{k},A_{k},Y_{k})$:

$\Gamma^{*}(O)={\mathbb{I}\left[A=1\right]\frac{Y-\mu^{*}(X)}{e^{*}(X)}+\mu^{*}(X)}.$

We will study the decomposition relative to the AIPW estimand with oracle
nuisance estimates, with $\Gamma(O;e,\mu)$ denoting the above score but with
nuisance functions $e,\mu$:

$\displaystyle\Psi_{n}(\tilde{P})-\mathbb{E}[\Gamma^{*}(O;e^{*},\mu^{*})]$
$\displaystyle=\textstyle(\Psi(\tilde{P})+\frac{1}{n}\sum_{{{i}}=1}^{n}{{\tilde{\phi}}}_{{i}})-\mathbb{E}[\Gamma^{*}(O;e^{*},\mu^{*})]$
$\displaystyle=\mathbb{E}_{n}[\Gamma(O;\tilde{e}_{\epsilon},\tilde{\mu})]-\mathbb{E}[\Gamma^{*}(O;e^{*},\mu^{*})]+\mathbb{E}_{n}[\tilde{\mu}_{\epsilon}(X)-\tilde{\mu}(X)]+O(\lambda^{\beta})$
$\displaystyle=(\mathbb{E}_{n}-\mathbb{E})[\Gamma(O;\tilde{e}_{\epsilon},\tilde{\mu})]+\mathbb{E}[\Gamma(O;\tilde{e}_{\epsilon},\tilde{\mu})-\Gamma^{*}(O;e^{*},\mu^{*})]+\mathbb{E}_{n}[\tilde{\mu}_{\epsilon}(X)-\tilde{\mu}(X)]+O(\lambda^{\beta}),$

where the second equality follows by Proposition 1, and the third by
adding/subtracting $\mathbb{E}[\Gamma(O;\tilde{e}_{\epsilon},\tilde{\mu})].$

The first term,
$(\mathbb{E}_{n}-\mathbb{E})[\Gamma(O;\tilde{e}_{\epsilon},\tilde{\mu})]$ is
the standard empirical process term of the decomposition: by the Donsker
assumption and classical stability results for Donsker classes [Van der Vaart,
2000].

$(\mathbb{E}_{n}-\mathbb{E})[\Gamma(O;\tilde{e}_{\epsilon},\tilde{\mu})]={(\mathbb{E}_{n}-\mathbb{E})[\Gamma(O;{e}^{*},\mu^{*})]+o_{p}(n^{-\frac{1}{2}}).}$

Hence, we have

$\Psi_{n}(\tilde{P})-\mathbb{E}[\Gamma^{*}(O;e^{*},\mu^{*})]=\mathbb{E}[\Gamma(O;\tilde{e}_{\epsilon},\tilde{\mu})-\Gamma^{*}(O;e^{*},\mu^{*})]+\mathbb{E}_{n}[\tilde{\mu}_{\epsilon}(X)-\tilde{\mu}(X)]+O(\lambda^{\beta})+o_{p}(n^{-\frac{1}{2}}).$
(10)

We obtain the product error conditions by standard decomposition of the
remaining first term.

$\displaystyle\mathbb{E}[\Gamma(O;\tilde{e}_{\epsilon},\tilde{\mu})-\Gamma^{*}(O;e^{*},\mu^{*})]$
$\displaystyle=\mathbb{E}\left[\left(\mathbb{I}\left[A=1\right]\frac{Y-\tilde{\mu}(X)}{\tilde{e}_{\epsilon}(X)}+\tilde{\mu}(X)\right)\right]-\mathbb{E}\left[\left(\mathbb{I}\left[A=1\right]\frac{Y-\mu^{*}(X)}{e^{*}(X)}+\mu^{*}(X)\right)\right]$
$\displaystyle=\mathbb{E}\left[\mathbb{I}\left[A=1\right]\left(\frac{1}{\tilde{e}_{\epsilon}(X)}-\frac{1}{e^{*}(X)}\right)(Y-\mu^{*}(X))\right]-\mathbb{E}\left[\left(\frac{1}{e^{*}(X)}-1\right)(\tilde{\mu}(X)-\mu^{*}(X))\right]$
$\displaystyle\qquad+\mathbb{E}\left[\mathbb{I}\left[A=1\right]\left(\frac{1}{\tilde{e}_{\epsilon}(X)}-\frac{1}{e^{*}(X)}\right)(\tilde{\mu}(X)-\mu^{*}(X))\right]$
$\displaystyle=\mathbb{E}\left[\mathbb{I}\left[A=1\right]\left(\frac{1}{\tilde{e}_{\epsilon}(X)}-\frac{1}{e^{*}(X)}\right)(\tilde{\mu}(X)-\mu^{*}(X))\right]$
$\displaystyle\leq\nu\|\tilde{e}_{\epsilon}(X)-{e^{*}(X)}\|_{2,P}\times\|\tilde{\mu}(X)-\mu^{*}(X)\|_{2,P}$

The second equality follows by adding/subtracting
$\mathbb{I}\left[A=1\right]/\tilde{e}_{\epsilon},\mathbb{I}\left[A=1\right]\tilde{\mu}/e^{*},\mathbb{I}\left[A=1\right]\mu^{*}/e^{*}$).
The third equality follows from iterated expectations. The last inequality
follows from the Cauchy-Schwarz inequality and overlap. Therefore, $\epsilon$
must satisfy a rate so that $\tilde{e}_{\epsilon}$ converges at typical rate-
double robustness fast rates, such as $n^{-1/4}$; under a suitable sample
splitting scheme.

For the last term of Equation 10, note that by Lemma 2, the numerical error
$\epsilon$ is multiplied by error rate of the outcome nuisance.

$\displaystyle\mathbb{E}[\tilde{\mu}_{\epsilon}(X)-\tilde{\mu}(X)]$
$\displaystyle=\epsilon\int\int\frac{\tilde{p}(x)}{\tilde{p}_{\epsilon}(A=1,x)}\mathbb{I}\left[a_{i}=1\right]\tilde{\delta}_{(x_{i})}^{\lambda}(x)\;\left\\{y_{i}-\mathbb{E}_{\tilde{p}}[Y\mid
A=1,x]\right\\}\mathrm{d}x$
$\displaystyle=\epsilon\;\mathbb{E}\left[\frac{\mathbb{I}\left[A=1\right]}{\tilde{p}_{\epsilon}(A=1\mid
X)}\;\left\\{Y-\mathbb{E}_{\tilde{p}}[Y\mid
A=1,X]\right\\}\right]+O(\lambda^{\beta})$
$\displaystyle\leq\|{\mathbb{I}\left[A=1\right]}/{\tilde{p}_{\epsilon}(A=1\mid
X)}\|_{2,P}\|{Y-\mathbb{E}_{\tilde{p}}[Y\mid
A=1,X]}\|_{2,P}+O(\lambda^{\beta})$
$\displaystyle=O(\lambda^{\beta})+O_{p}(n^{-r_{\mu}})\epsilon$

by application of Cauchy-Schwarz inequality and strong overlap. The claim
follows by application of Lemma 1 to deduce implications on $\epsilon$ from
standard product-error rate requirements on $\tilde{\mu},\tilde{e}$. ∎

##  Appendix C Proofs for Section 5

###  C.1 Empirical Gateaux derivative for multi-stage dynamic treatment
regimes

We study the estimation properties of the empirical Gateaux derivative for
multi-stage DTR from the probability density representation. In $T$-stage
dynamic treatment regimes, the causal quantity of interest is the mean
potential outcome $\mathbb{E}\left[Y^{\bar{a}}\right]$, where
$\bar{a}=(a_{0},\ldots,a_{T})$ is the (deterministic) treatment strategy.
Assuming $Y^{\bar{a}}$ is sequentially ignorable given treatment and covariate
history $(\bar{A}_{t},\bar{X}_{t})$ at each time $t$, this causal quantity can
be identified by the $g$-formula

$\mathbb{E}\left[Y^{\bar{a}}\right]=\int\mathbb{E}\left[Y\,|\,\bar{A}=\bar{a},\bar{X}=\bar{x}\right]\prod_{t=1}^{T}\tilde{p}(x_{t}\,|\,\bar{a}_{t-1},\bar{x}_{t-1})\mathop{}\\!\mathrm{d}\bar{x}$

where $\bar{A}=(A_{0},\ldots,A_{t}),\bar{X}=(X_{0},\ldots,X_{t})$.

To derive its influence function, we take the empirical Gateaux derivative by
considering the perturbation in the direction of
$o_{i}=(\bar{x}_{i},\bar{a}_{i},y_{i})$, where
$\bar{x}_{i}=(x_{i0},\ldots,x_{it})$ and $\bar{a}_{i}=(a_{i0},\ldots,a_{it})$.

$\displaystyle\Psi(\tilde{P}^{i}_{\epsilon})$ $\displaystyle=\int\left(\int
y\cdot\frac{\tilde{p}_{\epsilon}(y,\bar{A}=\bar{a},\bar{X}=\bar{x})}{\tilde{p}_{\epsilon}(\bar{A}=\bar{a},\bar{X}=\bar{x})}\mathop{}\\!\mathrm{d}y\right)\prod_{t=1}^{T}\frac{\tilde{p}_{\epsilon}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}_{\epsilon}(\bar{a}_{t-1},\bar{x}_{t-1}))}\cdot\tilde{p}_{\epsilon}(x_{0})\mathop{}\\!\mathrm{d}\bar{x}.$
(11)

Below we characterize how this empirical Gateaux derivative at the smoothed
distribution differs from the one at the (unsmoothed) estimated distribution.
The following result is analogous to Proposition 1 but is extended to the
dynamic treatment regimes.

######  Proposition 5.

$\displaystyle\frac{\Psi(\tilde{P}^{i}_{\epsilon})-\Psi(\tilde{P})}{\epsilon}$
$\displaystyle=\left[\int\mathbb{E}_{\tilde{p}_{\epsilon}}\left[Y\,|\,\bar{A}=\bar{a},\bar{X}_{1:T}=\bar{x}_{1:T},X_{0}=x_{i0}\right]\prod_{t=1}^{T}\tilde{p}(x_{t}\,|\,\bar{a}_{t-1},\bar{x}_{1:t-1},x_{i0})\mathop{}\\!\mathrm{d}\bar{x}-\Psi(\tilde{P})\right]$
$\displaystyle\quad+\int\tilde{p}(x_{0})\left\\{\left[\frac{\tilde{\delta}^{\lambda}_{(\bar{a}_{i},\bar{x}_{i})}}{\tilde{p}_{\epsilon}(\bar{A}=\bar{a},\bar{X}=\bar{x})}{}\left[\int
y\cdot\tilde{\delta}^{\lambda}_{y_{i}}\mathop{}\\!\mathrm{d}y-\mathbb{E}_{\tilde{p}}\left[Y\,|\,\bar{A}=\bar{a},\bar{X}=\bar{x}\right]\right]\right]\prod_{t=1}^{T}\frac{\tilde{p}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}(\bar{a}_{t-1},\bar{x}_{t-1}))}\right\\}\mathop{}\\!\mathrm{d}\bar{x}$
$\displaystyle\quad+\sum_{s=1}^{T}\int\tilde{p}(x_{0})\left\\{\left(\int
y\cdot\frac{\tilde{p}_{\epsilon}(y,\bar{A}=\bar{a},\bar{X}=\bar{x})}{\tilde{p}_{\epsilon}(\bar{A}=\bar{a},\bar{X}=\bar{x})}\mathop{}\\!\mathrm{d}y\right)\left[\left\\{\prod_{t=s+1}^{T}\frac{\tilde{p}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}(\bar{a}_{t-1},\bar{x}_{t-1}))}\prod_{t=1}^{s-1}\frac{\tilde{p}_{\epsilon}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}_{\epsilon}(\bar{a}_{t-1},\bar{x}_{t-1}))}\right.\right.\right.$
$\displaystyle\qquad\qquad\times\left.\left.\left.\frac{1}{\tilde{p}_{\epsilon}(\bar{a}_{s-1},\bar{x}_{s-1}))}\left[\tilde{\delta}^{\lambda}_{(\bar{a}_{i,s-1},\bar{x}_{i,s})})-\frac{\tilde{p}(x_{s},\bar{a}_{s-1},\bar{x}_{s-1})}{\tilde{p}(\bar{a}_{s-1},\bar{x}_{s-1}))}\tilde{\delta}^{\lambda}_{(\bar{a}_{i,s-1},\bar{x}_{i,s-1})}\right]\right\\}\right]\right\\}\mathop{}\\!\mathrm{d}\bar{x}$
$\displaystyle\quad+\epsilon\sum_{s=1}^{T}\int\mathbb{E}_{\tilde{p}_{\epsilon}}\left[Y\,|\,\bar{A}=\bar{a},\bar{X}_{1:T}=\bar{x}_{1:T},X_{0}=x_{i0}\right]$
$\displaystyle\qquad\qquad\times\left[\left\\{\prod_{t=s+1}^{T}\frac{\tilde{p}(x_{t},\bar{a}_{t-1},\bar{x}_{1:t-1},x_{i0})}{\tilde{p}(\bar{a}_{t-1},\bar{x}_{1:t-1},x_{i0}))}\prod_{t=1}^{s-1}\frac{\tilde{p}_{\epsilon}(x_{t},\bar{a}_{t-1},\bar{x}_{1:t-1},x_{i0})}{\tilde{p}_{\epsilon}(\bar{a}_{t-1},\bar{x}_{1:t-1},x_{i0}))}\frac{1}{\tilde{p}_{\epsilon}(\bar{a}_{s-1},\bar{x}_{1:s-1},x_{i0}))}\right.\right.$
$\displaystyle\qquad\qquad\qquad\qquad\times\left.\left.\left[\tilde{\delta}^{\lambda}_{(\bar{a}_{i,s-1},\bar{x}_{i,s})})-\frac{\tilde{p}(x_{s},\bar{a}_{s-1},\bar{x}_{1:s-1},x_{i0})}{\tilde{p}(\bar{a}_{s-1},\bar{x}_{1:s-1},x_{i0}))}\tilde{\delta}^{\lambda}_{(\bar{a}_{i,s-1},\bar{x}_{i,s-1})}\right]\right\\}\right]\mathop{}\\!\mathrm{d}\bar{x}$
(12)

For comparison, the exact efficient influence function under the non-Markovian
dynamic treatment regime is

$\displaystyle\left[\mathbb{E}_{\tilde{p}(\bar{x}_{1:T})}\left[\mathbb{E}_{\tilde{p}}\left[Y\,|\,\bar{A}=\bar{a},\bar{X}_{1:T}=\bar{x}_{1:T},X_{0}=x_{i0}\right]\right]-\Psi(\tilde{P})\right]$
$\displaystyle\quad+\frac{\mathbb{I}{(\bar{a}_{i})}}{\tilde{p}(\bar{A}=\bar{a}\,|\,\bar{X}=\bar{x})}{}\left[Y_{i}-\mathbb{E}_{\tilde{p}(\bar{x})}\left[\mathbb{E}_{\tilde{p}}\left[Y\,|\,\bar{A}=\bar{a},\bar{X}=\bar{x}\right]\right]\right]$
$\displaystyle\quad+\sum_{s=1}^{T}\frac{\mathbb{I}{(\bar{a}_{i,T-s},\bar{x}_{i,T-s+1})}}{\tilde{p}_{\epsilon}(\bar{a}_{T-s}\,|\,\bar{x}_{T-s})}\left[\int\mathbb{E}_{\tilde{p}}\left[Y\,|\,\bar{A}=\bar{a},\bar{X}=\bar{x}\right]{\prod_{t=T-s+1}^{T-1}\tilde{p}(x_{t}\,|\,\bar{a}_{t-1},\bar{x}_{t-1})\mathop{}\\!\mathrm{d}\bar{x}_{(T-s+1):(T-1)}}\right]$
$\displaystyle\quad-\sum_{s=1}^{T}\frac{\mathbb{I}{(\bar{a}_{i,T-s},\bar{x}_{i,T-s})}}{\tilde{p}(\bar{a}_{T-s}\,|\,\bar{x}_{T-s})}\left[\int\mathbb{E}_{\tilde{p}}\left[Y\,|\,\bar{A}=\bar{a},\bar{X}=\bar{x}\right]\prod_{t=T-s}^{T-1}\tilde{p}(x_{t}\,|\,\bar{a}_{t-1},\bar{x}_{t-1})\mathop{}\\!\mathrm{d}\bar{x}_{(T-s):(T-1)}\right]$

This proposition performs a decomposition of the empirical derivative. It
quantifies the additive bias due to the smoothed perturbation compared to the
(raw) perturbation at the observation $o_{i}$.

###### Proof.

We first decompose the empirical Gateaux derivative.

$\displaystyle\frac{\Psi(\tilde{P}^{i}_{\epsilon})-\Psi(\tilde{P})}{\epsilon}$
$\displaystyle=\frac{1}{\epsilon}\left[\int\left(\int
y\cdot\frac{\tilde{p}_{\epsilon}(y,\bar{A}=\bar{a},\bar{X}=\bar{x})}{\tilde{p}_{\epsilon}(\bar{A}=\bar{a},\bar{X}=\bar{x})}\mathop{}\\!\mathrm{d}y\right)\prod_{t=1}^{T}\frac{\tilde{p}_{\epsilon}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}_{\epsilon}(\bar{a}_{t-1},\bar{x}_{t-1}))}\cdot[(1-\epsilon)\tilde{p}(x_{0})+\epsilon\tilde{\delta}^{\lambda}_{x_{i0}}]\mathop{}\\!\mathrm{d}\bar{x}\right.$
$\displaystyle\quad\left.-\int\left(\int
y\cdot\frac{\tilde{p}(y,\bar{A}=\bar{a},\bar{X}=\bar{x})}{\tilde{p}(\bar{A}=\bar{a},\bar{X}=\bar{x})}\mathop{}\\!\mathrm{d}y\right)\prod_{t=1}^{T}\frac{\tilde{p}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}(\bar{a}_{t-1},\bar{x}_{t-1}))}\cdot\tilde{p}(x_{0})\mathop{}\\!\mathrm{d}\bar{x}\right]$
(14) $\displaystyle=\int\left[\left(\int
y\cdot\frac{\tilde{p}_{\epsilon}(y,\bar{A}=\bar{a},\bar{X}=\bar{x})}{\tilde{p}_{\epsilon}(\bar{A}=\bar{a},\bar{X}=\bar{x})}\mathop{}\\!\mathrm{d}y\right)\prod_{t=1}^{T}\frac{\tilde{p}_{\epsilon}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}_{\epsilon}(\bar{a}_{t-1},\bar{x}_{t-1}))}\cdot\tilde{\delta}^{\lambda}_{x_{i0}}\right.$
$\displaystyle\quad-\left.\left(\int
y\cdot\frac{\tilde{p}(y,\bar{A}=\bar{a},\bar{X}=\bar{x})}{\tilde{p}(\bar{A}=\bar{a},\bar{X}=\bar{x})}\mathop{}\\!\mathrm{d}y\right)\prod_{t=1}^{T}\frac{\tilde{p}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}(\bar{a}_{t-1},\bar{x}_{t-1}))}\cdot\tilde{p}(x_{0})\right]\mathop{}\\!\mathrm{d}\bar{x}$
$\displaystyle\quad+\frac{1}{\epsilon}\left[\int\tilde{p}(x_{0})\left[\left(\int
y\cdot\frac{\tilde{p}_{\epsilon}(y,\bar{A}=\bar{a},\bar{X}=\bar{x})}{\tilde{p}_{\epsilon}(\bar{A}=\bar{a},\bar{X}=\bar{x})}\mathop{}\\!\mathrm{d}y\right)\prod_{t=1}^{T}\frac{\tilde{p}_{\epsilon}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}_{\epsilon}(\bar{a}_{t-1},\bar{x}_{t-1}))}\right.\right.$
$\displaystyle\quad\left.\left.-\left(\int
y\cdot\frac{\tilde{p}(y,\bar{A}=\bar{a},\bar{X}=\bar{x})}{\tilde{p}(\bar{A}=\bar{a},\bar{X}=\bar{x})}\mathop{}\\!\mathrm{d}y\right)\prod_{t=1}^{T}\frac{\tilde{p}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}(\bar{a}_{t-1},\bar{x}_{t-1})}\right]\mathop{}\\!\mathrm{d}\bar{x}\right]$
(15)

$\displaystyle=\left[\int\mathbb{E}_{\tilde{p}_{\epsilon}}\left[Y\,|\,\bar{A}=\bar{a},\bar{X}_{T}=\bar{x}_{T}\right]\prod_{t=1}^{T}\tilde{p}_{\epsilon}(x_{t}\,|\,\bar{a}_{t-1},\bar{x}_{t-1})\tilde{p}(x_{0})\tilde{\delta}^{\lambda}_{x_{i0}}\mathop{}\\!\mathrm{d}\bar{x}-\Psi(\tilde{P})\right]$
$\displaystyle\quad+\frac{1}{\epsilon}\left[\int\tilde{p}(x_{0})\left[\left(\int
y\cdot\frac{\tilde{p}_{\epsilon}(y,\bar{A}=\bar{a},\bar{X}=\bar{x})}{\tilde{p}_{\epsilon}(\bar{A}=\bar{a},\bar{X}=\bar{x})}\mathop{}\\!\mathrm{d}y\right)\prod_{t=1}^{T}\frac{\tilde{p}_{\epsilon}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}_{\epsilon}(\bar{a}_{t-1},\bar{x}_{t-1}))}\right.\right.$
$\displaystyle\quad\left.\left.-\left(\int
y\cdot\frac{\tilde{p}(y,\bar{A}=\bar{a},\bar{X}=\bar{x})}{\tilde{p}(\bar{A}=\bar{a},\bar{X}=\bar{x})}\mathop{}\\!\mathrm{d}y\right)\prod_{t=1}^{T}\frac{\tilde{p}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}(\bar{a}_{t-1},\bar{x}_{t-1}))}\right]\mathop{}\\!\mathrm{d}\bar{x}\right]$
$\displaystyle=\left[\int\mathbb{E}_{\tilde{p}_{\epsilon}}\left[Y\,|\,\bar{A}=\bar{a},\bar{X}_{1:T}=\bar{x}_{1:T},X_{0}=x_{i0}\right]\prod_{t=1}^{T}\tilde{p}_{\epsilon}(x_{t}\,|\,\bar{a}_{t-1},\bar{x}_{1:t-1},x_{i0})\mathop{}\\!\mathrm{d}\bar{x}-\Psi(\tilde{P})\right]$
$\displaystyle\quad+\frac{1}{\epsilon}\left[\int\tilde{p}(x_{0})\left[\left(\int
y\cdot\frac{\tilde{p}_{\epsilon}(y,\bar{A}=\bar{a},\bar{X}=\bar{x})}{\tilde{p}_{\epsilon}(\bar{A}=\bar{a},\bar{X}=\bar{x})}\mathop{}\\!\mathrm{d}y\right)\prod_{t=1}^{T}\frac{\tilde{p}_{\epsilon}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}_{\epsilon}(\bar{a}_{t-1},\bar{x}_{t-1}))}\right.\right.$
$\displaystyle\quad\left.\left.-\left(\int
y\cdot\frac{\tilde{p}(y,\bar{A}=\bar{a},\bar{X}=\bar{x})}{\tilde{p}(\bar{A}=\bar{a},\bar{X}=\bar{x})}\mathop{}\\!\mathrm{d}y\right)\prod_{t=1}^{T}\frac{\tilde{p}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}(\bar{a}_{t-1},\bar{x}_{t-1}))}\right]\mathop{}\\!\mathrm{d}\bar{x}\right]$
(16)

Equation 14 is due to the definition of the smoothed perturbation. Equation 15
separates the effect of smoothed perturbation into one in the evaluation of
$y$ and the other in the evaluation of the density ratios. Equation 16 is due
to the definition of $\mathbb{E}_{\tilde{p}_{\epsilon}}\left[\cdot\right]$.

##### Calculating the second term of Equation 16.

We first decompose this term into two difference terms.

$\displaystyle\left(\int
y\cdot\frac{\tilde{p}_{\epsilon}(y,\bar{A}=\bar{a},\bar{X}=\bar{x})}{\tilde{p}_{\epsilon}(\bar{A}=\bar{a},\bar{X}=\bar{x})}\mathop{}\\!\mathrm{d}y\right)\prod_{t=1}^{T}\frac{\tilde{p}_{\epsilon}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}_{\epsilon}(\bar{a}_{t-1},\bar{x}_{t-1}))}$
(17) $\displaystyle\qquad\qquad-\left(\int
y\cdot\frac{\tilde{p}(y,\bar{A}=\bar{a},\bar{X}=\bar{x})}{\tilde{p}(\bar{A}=\bar{a},\bar{X}=\bar{x})}\mathop{}\\!\mathrm{d}y\right)\prod_{t=1}^{T}\frac{\tilde{p}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}(\bar{a}_{t-1},\bar{x}_{t-1}))}$
$\displaystyle=\left(\int
y\cdot\frac{\tilde{p}_{\epsilon}(y,\bar{A}=\bar{a},\bar{X}=\bar{x})}{\tilde{p}_{\epsilon}(\bar{A}=\bar{a},\bar{X}=\bar{x})}\mathop{}\\!\mathrm{d}y\right)\left[\prod_{t=1}^{T}\frac{\tilde{p}_{\epsilon}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}_{\epsilon}(\bar{a}_{t-1},\bar{x}_{t-1}))}-\prod_{t=1}^{T}\frac{\tilde{p}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}(\bar{a}_{t-1},\bar{x}_{t-1}))}\right]$
$\displaystyle\quad+\left[\left(\int
y\cdot\frac{\tilde{p}_{\epsilon}(y,\bar{A}=\bar{a},\bar{X}=\bar{x})}{\tilde{p}_{\epsilon}(\bar{A}=\bar{a},\bar{X}=\bar{x})}\mathop{}\\!\mathrm{d}y\right)-\left(\int
y\cdot\frac{\tilde{p}(y,\bar{A}=\bar{a},\bar{X}=\bar{x})}{\tilde{p}(\bar{A}=\bar{a},\bar{X}=\bar{x})}\mathop{}\\!\mathrm{d}y\right)\right]\prod_{t=1}^{T}\frac{\tilde{p}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}(\bar{a}_{t-1},\bar{x}_{t-1}))}.$
(18)

The second difference term can be decomposed in a similar way as in Lemma 2.

$\displaystyle\left(\int
y\cdot\frac{\tilde{p}_{\epsilon}(y,\bar{A}=\bar{a},\bar{X}=\bar{x})}{\tilde{p}_{\epsilon}(\bar{A}=\bar{a},\bar{X}=\bar{x})}\mathop{}\\!\mathrm{d}y\right)-\left(\int
y\cdot\frac{\tilde{p}(y,\bar{A}=\bar{a},\bar{X}=\bar{x})}{\tilde{p}(\bar{A}=\bar{a},\bar{X}=\bar{x})}\mathop{}\\!\mathrm{d}y\right)$
$\displaystyle=\int
y\cdot\left[\frac{\tilde{p}_{\epsilon}(y,\bar{A}=\bar{a},\bar{X}=\bar{x})}{\tilde{p}_{\epsilon}(\bar{A}=\bar{a},\bar{X}=\bar{x})}-\frac{\tilde{p}(y,\bar{A}=\bar{a},\bar{X}=\bar{x})}{\tilde{p}(\bar{A}=\bar{a},\bar{X}=\bar{x})}\right]\mathop{}\\!\mathrm{d}y$
$\displaystyle=\int
y\cdot\left[\frac{(1-\epsilon)\tilde{p}(y,\bar{A}=\bar{a},\bar{X}=\bar{x})+\epsilon\tilde{\delta}^{\lambda}_{o_{i}}}{\tilde{p}_{\epsilon}(\bar{A}=\bar{a},\bar{X}=\bar{x})}-\frac{\tilde{p}(y,\bar{A}=\bar{a},\bar{X}=\bar{x})}{\tilde{p}(\bar{A}=\bar{a},\bar{X}=\bar{x})}\right]\mathop{}\\!\mathrm{d}y$
$\displaystyle=\int
y\cdot\left[\frac{((1-\epsilon)\tilde{p}(y,\bar{A}=\bar{a},\bar{X}=\bar{x})+\epsilon\tilde{\delta}^{\lambda}_{o_{i}})\tilde{p}(\bar{A}=\bar{a},\bar{X}=\bar{x})}{\tilde{p}_{\epsilon}(\bar{A}=\bar{a},\bar{X}=\bar{x})\tilde{p}(\bar{A}=\bar{a},\bar{X}=\bar{x})}\right.$
(19)
$\displaystyle\quad\quad\left.-\frac{\tilde{p}(y,\bar{A}=\bar{a},\bar{X}=\bar{x})((1-\epsilon)\tilde{p}(\bar{A}=\bar{a},\bar{X}=\bar{x})+\epsilon\tilde{\delta}^{\lambda}_{(\bar{a}_{i},\bar{x}_{i})})}{\tilde{p}_{\epsilon}(\bar{A}=\bar{a},\bar{X}=\bar{x})\tilde{p}(\bar{A}=\bar{a},\bar{X}=\bar{x})}\right]\mathop{}\\!\mathrm{d}y$
$\displaystyle=\int
y\cdot\left[\frac{(\epsilon\tilde{\delta}^{\lambda}_{o_{i}})\tilde{p}(\bar{A}=\bar{a},\bar{X}=\bar{x})}{\tilde{p}_{\epsilon}(\bar{A}=\bar{a},\bar{X}=\bar{x})\tilde{p}(\bar{A}=\bar{a},\bar{X}=\bar{x})}-\frac{\tilde{p}(y,\bar{A}=\bar{a},\bar{X}=\bar{x})(\epsilon\tilde{\delta}^{\lambda}_{(\bar{a}_{i},\bar{x}_{i})})}{\tilde{p}_{\epsilon}(\bar{A}=\bar{a},\bar{X}=\bar{x})\tilde{p}(\bar{A}=\bar{a},\bar{X}=\bar{x})}\right]\mathop{}\\!\mathrm{d}y$
$\displaystyle=\frac{\epsilon\tilde{\delta}^{\lambda}_{(\bar{a}_{i},\bar{x}_{i})}}{\tilde{p}_{\epsilon}(\bar{A}=\bar{a},\bar{X}=\bar{x})}{}\int
y\cdot\left[(\tilde{\delta}^{\lambda}_{y_{i}})-\frac{\tilde{p}(y,\bar{A}=\bar{a},\bar{X}=\bar{x})}{\tilde{p}(\bar{A}=\bar{a},\bar{X}=\bar{x})}\right]\mathop{}\\!\mathrm{d}y$
$\displaystyle=\frac{\epsilon\tilde{\delta}^{\lambda}_{(\bar{a}_{i},\bar{x}_{i})}}{\tilde{p}_{\epsilon}(\bar{A}=\bar{a},\bar{X}=\bar{x})}{}\left[\int
y\cdot\tilde{\delta}^{\lambda}_{y_{i}}\mathop{}\\!\mathrm{d}y-\mathbb{E}_{\tilde{p}}\left[Y\,|\,\bar{A}=\bar{a},\bar{X}=\bar{x}\right]\right]$

We next analyze the first difference term by induction.

$\displaystyle\prod_{t=1}^{T}\frac{\tilde{p}_{\epsilon}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}_{\epsilon}(\bar{a}_{t-1},\bar{x}_{t-1}))}-\prod_{t=1}^{T}\frac{\tilde{p}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}(\bar{a}_{t-1},\bar{x}_{t-1}))}$
$\displaystyle=\prod_{t=1}^{T-1}\frac{\tilde{p}_{\epsilon}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}_{\epsilon}(\bar{a}_{t-1},\bar{x}_{t-1}))}\left[\frac{\tilde{p}_{\epsilon}(x_{T},\bar{a}_{T-1},\bar{x}_{T-1})}{\tilde{p}_{\epsilon}(\bar{a}_{T-1},\bar{x}_{T-1}))}-\frac{\tilde{p}(x_{T},\bar{a}_{T-1},\bar{x}_{T-1})}{\tilde{p}(\bar{a}_{T-1},\bar{x}_{T-1}))}\right]$
$\displaystyle\quad+\frac{\tilde{p}(x_{T},\bar{a}_{T-1},\bar{x}_{T-1})}{\tilde{p}(\bar{a}_{T-1},\bar{x}_{T-1}))}\left[\prod_{t=1}^{T-1}\frac{\tilde{p}_{\epsilon}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}_{\epsilon}(\bar{a}_{t-1},\bar{x}_{t-1}))}-\prod_{t=1}^{T-1}\frac{\tilde{p}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}(\bar{a}_{t-1},\bar{x}_{t-1}))}\right]$
(21)

Next we evaluate the key term in this induction
$\frac{\tilde{p}_{\epsilon}(x_{s},\bar{a}_{s-1},\bar{x}_{s-1})}{\tilde{p}_{\epsilon}(\bar{a}_{s-1},\bar{x}_{s-1}))}-\frac{\tilde{p}(x_{s},\bar{a}_{s-1},\bar{x}_{s-1})}{\tilde{p}(\bar{a}_{s-1},\bar{x}_{s-1}))}$.

$\displaystyle\frac{\tilde{p}_{\epsilon}(x_{s},\bar{a}_{s-1},\bar{x}_{s-1})}{\tilde{p}_{\epsilon}(\bar{a}_{s-1},\bar{x}_{s-1}))}-\frac{\tilde{p}(x_{s},\bar{a}_{s-1},\bar{x}_{s-1})}{\tilde{p}(\bar{a}_{s-1},\bar{x}_{s-1}))}$
$\displaystyle=\frac{((1-\epsilon)\tilde{p}(x_{s},\bar{a}_{s-1},\bar{x}_{s-1})+\epsilon\tilde{\delta}^{\lambda}_{o_{i}})\tilde{p}(\bar{a}_{s-1},\bar{x}_{s-1}))}{\tilde{p}_{\epsilon}(\bar{a}_{s-1},\bar{x}_{s-1}))\tilde{p}(\bar{a}_{s-1},\bar{x}_{s-1}))}-\frac{\tilde{p}(x_{s},\bar{a}_{s-1},\bar{x}_{s-1})((1-\epsilon)\tilde{p}(\bar{a}_{s-1},\bar{x}_{s-1}))+\epsilon\tilde{\delta}^{\lambda}_{o_{i}})}{\tilde{p}(\bar{a}_{s-1},\bar{x}_{s-1}))\tilde{p}_{\epsilon}(\bar{a}_{s-1},\bar{x}_{s-1}))}$
$\displaystyle=\frac{(\epsilon\tilde{\delta}^{\lambda}_{(\bar{a}_{i,s-1},\bar{x}_{i,s})})}{\tilde{p}_{\epsilon}(\bar{a}_{s-1},\bar{x}_{s-1}))}-\frac{\tilde{p}(x_{s},\bar{a}_{s-1},\bar{x}_{s-1})(\epsilon\tilde{\delta}^{\lambda}_{(\bar{a}_{i,s-1},\bar{x}_{i,s-1})})}{\tilde{p}(\bar{a}_{s-1},\bar{x}_{s-1}))\tilde{p}_{\epsilon}(\bar{a}_{s-1},\bar{x}_{s-1}))}$
$\displaystyle=\frac{\epsilon}{\tilde{p}_{\epsilon}(\bar{a}_{s-1},\bar{x}_{s-1}))}\left[\tilde{\delta}^{\lambda}_{(\bar{a}_{i,s-1},\bar{x}_{i,s})})-\frac{\tilde{p}(x_{s},\bar{a}_{s-1},\bar{x}_{s-1})}{\tilde{p}(\bar{a}_{s-1},\bar{x}_{s-1}))}\tilde{\delta}^{\lambda}_{(\bar{a}_{i,s-1},\bar{x}_{i,s-1})}\right]$
(22)

Thus we have the following inductive relationship in $T$,

$\displaystyle\prod_{t=1}^{T}\frac{\tilde{p}_{\epsilon}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}_{\epsilon}(\bar{a}_{t-1},\bar{x}_{t-1}))}-\prod_{t=1}^{T}\frac{\tilde{p}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}(\bar{a}_{t-1},\bar{x}_{t-1}))}$
$\displaystyle=\prod_{t=1}^{T-1}\frac{\tilde{p}_{\epsilon}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}_{\epsilon}(\bar{a}_{t-1},\bar{x}_{t-1}))}\frac{\epsilon}{\tilde{p}_{\epsilon}(\bar{a}_{T-1},\bar{x}_{T-1}))}\left[\tilde{\delta}^{\lambda}_{(\bar{a}_{i,T-1},\bar{x}_{i,T})})-\frac{\tilde{p}(x_{T},\bar{a}_{T-1},\bar{x}_{T-1})}{\tilde{p}(\bar{a}_{T-1},\bar{x}_{T-1}))}\tilde{\delta}^{\lambda}_{(\bar{a}_{i,T-1},\bar{x}_{i,T-1})}\right]$
$\displaystyle\quad+\frac{\tilde{p}(x_{T},\bar{a}_{T-1},\bar{x}_{T-1})}{\tilde{p}(\bar{a}_{T-1},\bar{x}_{T-1}))}\left[\prod_{t=1}^{T-1}\frac{\tilde{p}_{\epsilon}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}_{\epsilon}(\bar{a}_{t-1},\bar{x}_{t-1}))}-\prod_{t=1}^{T-1}\frac{\tilde{p}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}(\bar{a}_{t-1},\bar{x}_{t-1}))}\right]$
(23)
$\displaystyle=\sum_{s=T-1}^{T}\left\\{\prod_{t=s+1}^{T}\frac{\tilde{p}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}(\bar{a}_{t-1},\bar{x}_{t-1}))}\prod_{t=1}^{s-1}\frac{\tilde{p}_{\epsilon}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}_{\epsilon}(\bar{a}_{t-1},\bar{x}_{t-1}))}\frac{\epsilon}{\tilde{p}_{\epsilon}(\bar{a}_{s-1},\bar{x}_{s-1}))}\right.$
$\displaystyle\qquad\left.\qquad\times\left[\tilde{\delta}^{\lambda}_{(\bar{a}_{i,s-1},\bar{x}_{i,s})})-\frac{\tilde{p}(x_{s},\bar{a}_{s-1},\bar{x}_{s-1})}{\tilde{p}(\bar{a}_{s-1},\bar{x}_{s-1}))}\tilde{\delta}^{\lambda}_{(\bar{a}_{i,s-1},\bar{x}_{i,s-1})}\right]\right\\}$
$\displaystyle\quad+\prod_{t=T-1}^{T}\frac{\tilde{p}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}(\bar{a}_{t-1},\bar{x}_{t-1}))}\left[\prod_{t=1}^{T-2}\frac{\tilde{p}_{\epsilon}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}_{\epsilon}(\bar{a}_{t-1},\bar{x}_{t-1}))}-\prod_{t=1}^{T-2}\frac{\tilde{p}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}(\bar{a}_{t-1},\bar{x}_{t-1}))}\right]$
(24)
$\displaystyle=\sum_{s=1}^{T}\left\\{\prod_{t=s+1}^{T}\frac{\tilde{p}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}(\bar{a}_{t-1},\bar{x}_{t-1}))}\prod_{t=1}^{s-1}\frac{\tilde{p}_{\epsilon}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}_{\epsilon}(\bar{a}_{t-1},\bar{x}_{t-1}))}\frac{\epsilon}{\tilde{p}_{\epsilon}(\bar{a}_{s-1},\bar{x}_{s-1}))}\right.$
$\displaystyle\qquad\qquad\times\left.\left[\tilde{\delta}^{\lambda}_{(\bar{a}_{i,s-1},\bar{x}_{i,s})})-\frac{\tilde{p}(x_{s},\bar{a}_{s-1},\bar{x}_{s-1})}{\tilde{p}(\bar{a}_{s-1},\bar{x}_{s-1}))}\tilde{\delta}^{\lambda}_{(\bar{a}_{i,s-1},\bar{x}_{i,s-1})}\right]\right\\}$
(25)

These calculations imply that

$\displaystyle\frac{\Psi(\tilde{P}^{i}_{\epsilon})-\Psi(\tilde{P})}{\epsilon}$
$\displaystyle=\left[\int\mathbb{E}_{\tilde{p}_{\epsilon}}\left[Y\,|\,\bar{A}=\bar{a},\bar{X}_{1:T}=\bar{x}_{1:T},X_{0}=x_{i0}\right]\prod_{t=1}^{T}\tilde{p}_{\epsilon}(x_{t}\,|\,\bar{a}_{t-1},\bar{x}_{1:t-1},x_{i0})\mathop{}\\!\mathrm{d}\bar{x}-\Psi(\tilde{P})\right]$
$\displaystyle\quad+\frac{1}{\epsilon}\left[\int\tilde{p}(x_{0})\left[\left(\int
y\cdot\frac{\tilde{p}_{\epsilon}(y,\bar{A}=\bar{a},\bar{X}=\bar{x})}{\tilde{p}_{\epsilon}(\bar{A}=\bar{a},\bar{X}=\bar{x})}\mathop{}\\!\mathrm{d}y\right)\prod_{t=1}^{T}\frac{\tilde{p}_{\epsilon}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}_{\epsilon}(\bar{a}_{t-1},\bar{x}_{t-1}))}\right.\right.$
$\displaystyle\quad\left.\left.-\left(\int
y\cdot\frac{\tilde{p}(y,\bar{A}=\bar{a},\bar{X}=\bar{x})}{\tilde{p}(\bar{A}=\bar{a},\bar{X}=\bar{x})}\mathop{}\\!\mathrm{d}y\right)\prod_{t=1}^{T}\frac{\tilde{p}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}(\bar{a}_{t-1},\bar{x}_{t-1}))}\right]\mathop{}\\!\mathrm{d}\bar{x}\right]$
(26)
$\displaystyle=\left[\int\mathbb{E}_{\tilde{p}_{\epsilon}}\left[Y\,|\,\bar{A}=\bar{a},\bar{X}_{1:T}=\bar{x}_{1:T},X_{0}=x_{i0}\right]\prod_{t=1}^{T}\tilde{p}_{\epsilon}(x_{t}\,|\,\bar{a}_{t-1},\bar{x}_{1:t-1},x_{i0})\mathop{}\\!\mathrm{d}\bar{x}-\Psi(\tilde{P})\right]$
$\displaystyle\quad+\frac{1}{\epsilon}\int\tilde{p}(x_{0})\left\\{\left(\int
y\cdot\frac{\tilde{p}_{\epsilon}(y,\bar{A}=\bar{a},\bar{X}=\bar{x})}{\tilde{p}_{\epsilon}(\bar{A}=\bar{a},\bar{X}=\bar{x})}\mathop{}\\!\mathrm{d}y\right)\left[\prod_{t=1}^{T}\frac{\tilde{p}_{\epsilon}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}_{\epsilon}(\bar{a}_{t-1},\bar{x}_{t-1}))}-\prod_{t=1}^{T}\frac{\tilde{p}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}(\bar{a}_{t-1},\bar{x}_{t-1}))}\right]\right.$
$\displaystyle\quad+\left.\left[\left(\int
y\cdot\frac{\tilde{p}_{\epsilon}(y,\bar{A}=\bar{a},\bar{X}=\bar{x})}{\tilde{p}_{\epsilon}(\bar{A}=\bar{a},\bar{X}=\bar{x})}\mathop{}\\!\mathrm{d}y\right)-\left(\int
y\cdot\frac{\tilde{p}(y,\bar{A}=\bar{a},\bar{X}=\bar{x})}{\tilde{p}(\bar{A}=\bar{a},\bar{X}=\bar{x})}\mathop{}\\!\mathrm{d}y\right)\right]\prod_{t=1}^{T}\frac{\tilde{p}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}(\bar{a}_{t-1},\bar{x}_{t-1}))}\right\\}\mathop{}\\!\mathrm{d}\bar{x}$
(27)
$\displaystyle=\left[\int\mathbb{E}_{\tilde{p}_{\epsilon}}\left[Y\,|\,\bar{A}=\bar{a},\bar{X}_{1:T}=\bar{x}_{1:T},X_{0}=x_{i0}\right]\prod_{t=1}^{T}\tilde{p}_{\epsilon}(x_{t}\,|\,\bar{a}_{t-1},\bar{x}_{1:t-1},x_{i0})\mathop{}\\!\mathrm{d}\bar{x}-\Psi(\tilde{P})\right]$
$\displaystyle\quad+\frac{1}{\epsilon}\int\tilde{p}(x_{0})\left\\{\left(\int
y\cdot\frac{\tilde{p}_{\epsilon}(y,\bar{A}=\bar{a},\bar{X}=\bar{x})}{\tilde{p}_{\epsilon}(\bar{A}=\bar{a},\bar{X}=\bar{x})}\mathop{}\\!\mathrm{d}y\right)\left[\prod_{t=1}^{T}\frac{\tilde{p}_{\epsilon}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}_{\epsilon}(\bar{a}_{t-1},\bar{x}_{t-1}))}-\prod_{t=1}^{T}\frac{\tilde{p}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}(\bar{a}_{t-1},\bar{x}_{t-1}))}\right]\right\\}\mathop{}\\!\mathrm{d}\bar{x}$
$\displaystyle\quad+\frac{1}{\epsilon}\int\tilde{p}(x_{0})\left\\{\left[\frac{\epsilon\tilde{\delta}^{\lambda}_{(\bar{a}_{i},\bar{x}_{i})}}{\tilde{p}_{\epsilon}(\bar{A}=\bar{a},\bar{X}=\bar{x})}{}\left[\int
y\cdot\tilde{\delta}^{\lambda}_{y_{i}}\mathop{}\\!\mathrm{d}y-\mathbb{E}_{\tilde{p}}\left[Y\,|\,\bar{A}=\bar{a},\bar{X}=\bar{x}\right]\right]\right]\prod_{t=1}^{T}\frac{\tilde{p}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}(\bar{a}_{t-1},\bar{x}_{t-1}))}\right\\}\mathop{}\\!\mathrm{d}\bar{x}$
(28)
$\displaystyle=\left[\int\mathbb{E}_{\tilde{p}_{\epsilon}}\left[Y\,|\,\bar{A}=\bar{a},\bar{X}_{1:T}=\bar{x}_{1:T},X_{0}=x_{i0}\right]\prod_{t=1}^{T}\tilde{p}_{\epsilon}(x_{t}\,|\,\bar{a}_{t-1},\bar{x}_{1:t-1},x_{i0})\mathop{}\\!\mathrm{d}\bar{x}-\Psi(\tilde{P})\right]$
$\displaystyle\quad+\frac{1}{\epsilon}\int\tilde{p}(x_{0})\left\\{\left(\int
y\cdot\frac{\tilde{p}_{\epsilon}(y,\bar{A}=\bar{a},\bar{X}=\bar{x})}{\tilde{p}_{\epsilon}(\bar{A}=\bar{a},\bar{X}=\bar{x})}\mathop{}\\!\mathrm{d}y\right)\left[\sum_{s=1}^{T}\left\\{\prod_{t=s+1}^{T}\frac{\tilde{p}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}(\bar{a}_{t-1},\bar{x}_{t-1}))}\prod_{t=1}^{s-1}\frac{\tilde{p}_{\epsilon}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}_{\epsilon}(\bar{a}_{t-1},\bar{x}_{t-1}))}\right.\right.\right.$
$\displaystyle\qquad\qquad\times\left.\left.\left.\frac{\epsilon}{\tilde{p}_{\epsilon}(\bar{a}_{s-1},\bar{x}_{s-1}))}\left[\tilde{\delta}^{\lambda}_{(\bar{a}_{i,s-1},\bar{x}_{i,s})})-\frac{\tilde{p}(x_{s},\bar{a}_{s-1},\bar{x}_{s-1})}{\tilde{p}(\bar{a}_{s-1},\bar{x}_{s-1}))}\tilde{\delta}^{\lambda}_{(\bar{a}_{i,s-1},\bar{x}_{i,s-1})}\right]\right\\}\right]\right\\}\mathop{}\\!\mathrm{d}\bar{x}$
$\displaystyle\quad+\frac{1}{\epsilon}\int\tilde{p}(x_{0})\left\\{\left[\frac{\epsilon\tilde{\delta}^{\lambda}_{(\bar{a}_{i},\bar{x}_{i})}}{\tilde{p}_{\epsilon}(\bar{A}=\bar{a},\bar{X}=\bar{x})}{}\left[\int
y\cdot\tilde{\delta}^{\lambda}_{y_{i}}\mathop{}\\!\mathrm{d}y-\mathbb{E}_{\tilde{p}}\left[Y\,|\,\bar{A}=\bar{a},\bar{X}=\bar{x}\right]\right]\right]\prod_{t=1}^{T}\frac{\tilde{p}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}(\bar{a}_{t-1},\bar{x}_{t-1}))}\right\\}\mathop{}\\!\mathrm{d}\bar{x}$
(29)

$\displaystyle=\left[\int\mathbb{E}_{\tilde{p}_{\epsilon}}\left[Y\,|\,\bar{A}=\bar{a},\bar{X}_{1:T}=\bar{x}_{1:T},X_{0}=x_{i0}\right]\prod_{t=1}^{T}\tilde{p}_{\epsilon}(x_{t}\,|\,\bar{a}_{t-1},\bar{x}_{1:t-1},x_{i0})\mathop{}\\!\mathrm{d}\bar{x}-\Psi(\tilde{P})\right]$
$\displaystyle\quad+\int\tilde{p}(x_{0})\left\\{\left(\int
y\cdot\frac{\tilde{p}_{\epsilon}(y,\bar{A}=\bar{a},\bar{X}=\bar{x})}{\tilde{p}_{\epsilon}(\bar{A}=\bar{a},\bar{X}=\bar{x})}\mathop{}\\!\mathrm{d}y\right)\left[\sum_{s=1}^{T}\left\\{\prod_{t=s+1}^{T}\frac{\tilde{p}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}(\bar{a}_{t-1},\bar{x}_{t-1}))}\prod_{t=1}^{s-1}\frac{\tilde{p}_{\epsilon}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}_{\epsilon}(\bar{a}_{t-1},\bar{x}_{t-1}))}\right.\right.\right.$
$\displaystyle\qquad\qquad\times\left.\left.\left.\frac{1}{\tilde{p}_{\epsilon}(\bar{a}_{s-1},\bar{x}_{s-1}))}\left[\tilde{\delta}^{\lambda}_{(\bar{a}_{i,s-1},\bar{x}_{i,s})})-\frac{\tilde{p}(x_{s},\bar{a}_{s-1},\bar{x}_{s-1})}{\tilde{p}(\bar{a}_{s-1},\bar{x}_{s-1}))}\tilde{\delta}^{\lambda}_{(\bar{a}_{i,s-1},\bar{x}_{i,s-1})}\right]\right\\}\right]\right\\}\mathop{}\\!\mathrm{d}\bar{x}$
$\displaystyle\quad+\int\tilde{p}(x_{0})\left\\{\left[\frac{\tilde{\delta}^{\lambda}_{(\bar{a}_{i},\bar{x}_{i})}}{\tilde{p}_{\epsilon}(\bar{A}=\bar{a},\bar{X}=\bar{x})}{}\left[\int
y\cdot\tilde{\delta}^{\lambda}_{y_{i}}\mathop{}\\!\mathrm{d}y-\mathbb{E}_{\tilde{p}}\left[Y\,|\,\bar{A}=\bar{a},\bar{X}=\bar{x}\right]\right]\right]\prod_{t=1}^{T}\frac{\tilde{p}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}(\bar{a}_{t-1},\bar{x}_{t-1}))}\right\\}\mathop{}\\!\mathrm{d}\bar{x}$
(30)

##### Calculating the first term of Equation 30.

The first term of Equation 30 can be decomposed as follows.

$\displaystyle\int\mathbb{E}_{\tilde{p}_{\epsilon}}\left[Y\,|\,\bar{A}=\bar{a},\bar{X}_{1:T}=\bar{x}_{1:T},X_{0}=x_{i0}\right]\prod_{t=1}^{T}\tilde{p}_{\epsilon}(x_{t}\,|\,\bar{a}_{t-1},\bar{x}_{1:t-1},x_{i0})\mathop{}\\!\mathrm{d}\bar{x}$
$\displaystyle=\int\mathbb{E}_{\tilde{p}_{\epsilon}}\left[Y\,|\,\bar{A}=\bar{a},\bar{X}_{1:T}=\bar{x}_{1:T},X_{0}=x_{i0}\right]\prod_{t=1}^{T}\tilde{p}(x_{t}\,|\,\bar{a}_{t-1},\bar{x}_{1:t-1},x_{i0})\mathop{}\\!\mathrm{d}\bar{x}$
(31)
$\displaystyle\qquad+\int\mathbb{E}_{\tilde{p}_{\epsilon}}\left[Y\,|\,\bar{A}=\bar{a},\bar{X}_{1:T}=\bar{x}_{1:T},X_{0}=x_{i0}\right]$
$\displaystyle\qquad\qquad\times\left[\prod_{t=1}^{T}\frac{\tilde{p}_{\epsilon}(x_{t},\bar{a}_{t-1},\bar{x}_{1:t-1},x_{i0})}{\tilde{p}_{\epsilon}(\bar{a}_{t-1},\bar{x}_{1:t-1},x_{i0}))}-\prod_{t=1}^{T}\frac{\tilde{p}(x_{t},\bar{a}_{t-1},\bar{x}_{1:t-1},x_{i0})}{\tilde{p}(\bar{a}_{t-1},\bar{x}_{1:t-1},x_{i0}))}\right]\mathop{}\\!\mathrm{d}\bar{x}$
(32)

Following a similar calculation as in Equation 25, we have that

$\displaystyle\prod_{t=1}^{T}\frac{\tilde{p}_{\epsilon}(x_{t},\bar{a}_{t-1},\bar{x}_{1:t-1},x_{i0})}{\tilde{p}_{\epsilon}(\bar{a}_{t-1},\bar{x}_{1:t-1},x_{i0}))}-\prod_{t=1}^{T}\frac{\tilde{p}(x_{t},\bar{a}_{t-1},\bar{x}_{1:t-1},x_{i0})}{\tilde{p}(\bar{a}_{t-1},\bar{x}_{1:t-1},x_{i0}))}$
(33)
$\displaystyle=\sum_{s=1}^{T}\left\\{\prod_{t=s+1}^{T}\frac{\tilde{p}(x_{t},\bar{a}_{t-1},\bar{x}_{1:t-1},x_{i0})}{\tilde{p}(\bar{a}_{t-1},\bar{x}_{1:t-1},x_{i0}))}\prod_{t=1}^{s-1}\frac{\tilde{p}_{\epsilon}(x_{t},\bar{a}_{t-1},\bar{x}_{1:t-1},x_{i0})}{\tilde{p}_{\epsilon}(\bar{a}_{t-1},\bar{x}_{1:t-1},x_{i0}))}\frac{\epsilon}{\tilde{p}_{\epsilon}(\bar{a}_{s-1},\bar{x}_{1:s-1},x_{i0}))}\right.$
$\displaystyle\qquad\qquad\times\left.\left[\tilde{\delta}^{\lambda}_{(\bar{a}_{i,s-1},\bar{x}_{i,s})})-\frac{\tilde{p}(x_{s},\bar{a}_{s-1},\bar{x}_{1:s-1},x_{i0})}{\tilde{p}(\bar{a}_{s-1},\bar{x}_{1:s-1},x_{i0}))}\tilde{\delta}^{\lambda}_{(\bar{a}_{i,s-1},\bar{x}_{i,s-1})}\right]\right\\}$
(34)

Thus the first term of Equation 30 is

$\displaystyle\int\mathbb{E}_{\tilde{p}_{\epsilon}}\left[Y\,|\,\bar{A}=\bar{a},\bar{X}_{1:T}=\bar{x}_{1:T},X_{0}=x_{i0}\right]\prod_{t=1}^{T}\tilde{p}_{\epsilon}(x_{t}\,|\,\bar{a}_{t-1},\bar{x}_{1:t-1},x_{i0})\mathop{}\\!\mathrm{d}\bar{x}$
$\displaystyle=\int\mathbb{E}_{\tilde{p}_{\epsilon}}\left[Y\,|\,\bar{A}=\bar{a},\bar{X}_{1:T}=\bar{x}_{1:T},X_{0}=x_{i0}\right]\prod_{t=1}^{T}\tilde{p}(x_{t}\,|\,\bar{a}_{t-1},\bar{x}_{1:t-1},x_{i0})\mathop{}\\!\mathrm{d}\bar{x}$
(35)
$\displaystyle\qquad+\int\mathbb{E}_{\tilde{p}_{\epsilon}}\left[Y\,|\,\bar{A}=\bar{a},\bar{X}_{1:T}=\bar{x}_{1:T},X_{0}=x_{i0}\right]$
(36)
$\displaystyle\qquad\qquad\times\left[\sum_{s=1}^{T}\left\\{\prod_{t=s+1}^{T}\frac{\tilde{p}(x_{t},\bar{a}_{t-1},\bar{x}_{1:t-1},x_{i0})}{\tilde{p}(\bar{a}_{t-1},\bar{x}_{1:t-1},x_{i0}))}\prod_{t=1}^{s-1}\frac{\tilde{p}_{\epsilon}(x_{t},\bar{a}_{t-1},\bar{x}_{1:t-1},x_{i0})}{\tilde{p}_{\epsilon}(\bar{a}_{t-1},\bar{x}_{1:t-1},x_{i0}))}\right.\right.$
$\displaystyle\qquad\qquad\times\left.\left.\frac{\epsilon}{\tilde{p}_{\epsilon}(\bar{a}_{s-1},\bar{x}_{1:s-1},x_{i0}))}\left[\tilde{\delta}^{\lambda}_{(\bar{a}_{i,s-1},\bar{x}_{i,s})})-\frac{\tilde{p}(x_{s},\bar{a}_{s-1},\bar{x}_{1:s-1},x_{i0})}{\tilde{p}(\bar{a}_{s-1},\bar{x}_{1:s-1},x_{i0}))}\tilde{\delta}^{\lambda}_{(\bar{a}_{i,s-1},\bar{x}_{i,s-1})}\right]\right\\}\right]\mathop{}\\!\mathrm{d}\bar{x}$
(37)

Plugging this calculation into
$\frac{\Psi(\tilde{P}^{i}_{\epsilon})-\Psi(\tilde{P})}{\epsilon}:$

$\displaystyle\frac{\Psi(\tilde{P}^{i}_{\epsilon})-\Psi(\tilde{P})}{\epsilon}$
$\displaystyle=\left[\int\mathbb{E}_{\tilde{p}_{\epsilon}}\left[Y\,|\,\bar{A}=\bar{a},\bar{X}_{1:T}=\bar{x}_{1:T},X_{0}=x_{i0}\right]\prod_{t=1}^{T}\tilde{p}_{\epsilon}(x_{t}\,|\,\bar{a}_{t-1},\bar{x}_{1:t-1},x_{i0})\mathop{}\\!\mathrm{d}\bar{x}-\Psi(\tilde{P})\right]$
$\displaystyle\quad+\int\tilde{p}(x_{0})\left\\{\left(\int
y\cdot\frac{\tilde{p}_{\epsilon}(y,\bar{A}=\bar{a},\bar{X}=\bar{x})}{\tilde{p}_{\epsilon}(\bar{A}=\bar{a},\bar{X}=\bar{x})}\mathop{}\\!\mathrm{d}y\right)\left[\sum_{s=1}^{T}\left\\{\prod_{t=s+1}^{T}\frac{\tilde{p}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}(\bar{a}_{t-1},\bar{x}_{t-1}))}\prod_{t=1}^{s-1}\frac{\tilde{p}_{\epsilon}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}_{\epsilon}(\bar{a}_{t-1},\bar{x}_{t-1}))}\right.\right.\right.$
$\displaystyle\qquad\qquad\times\left.\left.\left.\frac{1}{\tilde{p}_{\epsilon}(\bar{a}_{s-1},\bar{x}_{s-1}))}\left[\tilde{\delta}^{\lambda}_{(\bar{a}_{i,s-1},\bar{x}_{i,s})})-\frac{\tilde{p}(x_{s},\bar{a}_{s-1},\bar{x}_{s-1})}{\tilde{p}(\bar{a}_{s-1},\bar{x}_{s-1}))}\tilde{\delta}^{\lambda}_{(\bar{a}_{i,s-1},\bar{x}_{i,s-1})}\right]\right\\}\right]\right\\}\mathop{}\\!\mathrm{d}\bar{x}$
$\displaystyle\quad+\int\tilde{p}(x_{0})\left\\{\left[\frac{\tilde{\delta}^{\lambda}_{(\bar{a}_{i},\bar{x}_{i})}}{\tilde{p}_{\epsilon}(\bar{A}=\bar{a},\bar{X}=\bar{x})}{}\left[\int
y\cdot\tilde{\delta}^{\lambda}_{y_{i}}\mathop{}\\!\mathrm{d}y-\mathbb{E}_{\tilde{p}}\left[Y\,|\,\bar{A}=\bar{a},\bar{X}=\bar{x}\right]\right]\right]\prod_{t=1}^{T}\frac{\tilde{p}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}(\bar{a}_{t-1},\bar{x}_{t-1}))}\right\\}\mathop{}\\!\mathrm{d}\bar{x}$
(38)
$\displaystyle=\left[\int\mathbb{E}_{\tilde{p}_{\epsilon}}\left[Y\,|\,\bar{A}=\bar{a},\bar{X}_{1:T}=\bar{x}_{1:T},X_{0}=x_{i0}\right]\prod_{t=1}^{T}\tilde{p}(x_{t}\,|\,\bar{a}_{t-1},\bar{x}_{1:t-1},x_{i0})\mathop{}\\!\mathrm{d}\bar{x}-\Psi(\tilde{P})\right]$
$\displaystyle\quad+\int\tilde{p}(x_{0})\left\\{\left[\frac{\tilde{\delta}^{\lambda}_{(\bar{a}_{i},\bar{x}_{i})}}{\tilde{p}_{\epsilon}(\bar{A}=\bar{a},\bar{X}=\bar{x})}{}\left[\int
y\cdot\tilde{\delta}^{\lambda}_{y_{i}}\mathop{}\\!\mathrm{d}y-\mathbb{E}_{\tilde{p}}\left[Y\,|\,\bar{A}=\bar{a},\bar{X}=\bar{x}\right]\right]\right]\prod_{t=1}^{T}\frac{\tilde{p}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}(\bar{a}_{t-1},\bar{x}_{t-1}))}\right\\}\mathop{}\\!\mathrm{d}\bar{x}$
$\displaystyle\quad+\sum_{s=1}^{T}\int\tilde{p}(x_{0})\left\\{\left(\int
y\cdot\frac{\tilde{p}_{\epsilon}(y,\bar{A}=\bar{a},\bar{X}=\bar{x})}{\tilde{p}_{\epsilon}(\bar{A}=\bar{a},\bar{X}=\bar{x})}\mathop{}\\!\mathrm{d}y\right)\left[\left\\{\prod_{t=s+1}^{T}\frac{\tilde{p}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}(\bar{a}_{t-1},\bar{x}_{t-1}))}\prod_{t=1}^{s-1}\frac{\tilde{p}_{\epsilon}(x_{t},\bar{a}_{t-1},\bar{x}_{t-1})}{\tilde{p}_{\epsilon}(\bar{a}_{t-1},\bar{x}_{t-1}))}\right.\right.\right.$
$\displaystyle\qquad\qquad\times\left.\left.\left.\frac{1}{\tilde{p}_{\epsilon}(\bar{a}_{s-1},\bar{x}_{s-1}))}\left[\tilde{\delta}^{\lambda}_{(\bar{a}_{i,s-1},\bar{x}_{i,s})})-\frac{\tilde{p}(x_{s},\bar{a}_{s-1},\bar{x}_{s-1})}{\tilde{p}(\bar{a}_{s-1},\bar{x}_{s-1}))}\tilde{\delta}^{\lambda}_{(\bar{a}_{i,s-1},\bar{x}_{i,s-1})}\right]\right\\}\right]\right\\}\mathop{}\\!\mathrm{d}\bar{x}$
$\displaystyle\quad+\epsilon\sum_{s=1}^{T}\int\mathbb{E}_{\tilde{p}_{\epsilon}}\left[Y\,|\,\bar{A}=\bar{a},\bar{X}_{1:T}=\bar{x}_{1:T},X_{0}=x_{i0}\right]$
(39)
$\displaystyle\qquad\qquad\times\left[\left\\{\prod_{t=s+1}^{T}\frac{\tilde{p}(x_{t},\bar{a}_{t-1},\bar{x}_{1:t-1},x_{i0})}{\tilde{p}(\bar{a}_{t-1},\bar{x}_{1:t-1},x_{i0}))}\prod_{t=1}^{s-1}\frac{\tilde{p}_{\epsilon}(x_{t},\bar{a}_{t-1},\bar{x}_{1:t-1},x_{i0})}{\tilde{p}_{\epsilon}(\bar{a}_{t-1},\bar{x}_{1:t-1},x_{i0}))}\right.\right.$
$\displaystyle\quad\qquad\qquad\times\left.\left.\frac{1}{\tilde{p}_{\epsilon}(\bar{a}_{s-1},\bar{x}_{1:s-1},x_{i0}))}\left[\tilde{\delta}^{\lambda}_{(\bar{a}_{i,s-1},\bar{x}_{i,s})})-\frac{\tilde{p}(x_{s},\bar{a}_{s-1},\bar{x}_{1:s-1},x_{i0})}{\tilde{p}(\bar{a}_{s-1},\bar{x}_{1:s-1},x_{i0}))}\tilde{\delta}^{\lambda}_{(\bar{a}_{i,s-1},\bar{x}_{i,s-1})}\right]\right\\}\right]\mathop{}\\!\mathrm{d}\bar{x}$
(40)

∎

###  C.2 Proofs for infinite-horizon evaluation

###### Proof of Proposition 3.

Background on sensitivity analysis of linear programs to coefficient matrix
perturbations. Our proof strategy broadly follows from that of Freund [1985],
which studies perturbations to linear programs of the constraint coefficient
matrix. However our perturbation is entrywise upon the transition probability
(a conditional transition probability); hence nonlinear in $\epsilon$; we
compute the entrywise matrix derivative directly using an alternative
technique described in the proof of Freund [1985, Thm.1 ]. But, the overall
argument for sensitivity analysis subject to left-hand-side perturbations
based on active set identification is similar, so we recount the argument here
as a preface.

Namely, note that for $\epsilon$ small enough, the active basis remains the
same when solving with respect to unperturbed or perturbed distributions,
$P_{0}$ or $P_{\epsilon}$. Consider the primal program in standard form,

$d(\epsilon)=\max\\{c^{\top}x\colon A^{\epsilon}x=b\\}.$

Freund [1985] considers perturbations where $A^{\epsilon}=A+\epsilon G$.
Letting $\beta$ denote an index vector for the optimal basis, note that

$d(\epsilon)=c_{\beta}x^{*}(\epsilon)=c_{\beta}(A^{\epsilon}_{\beta})^{-1}b_{\beta}$

by nondegeneracy and properties of the basic feasible solution for linear
programming. Freund [1985] notes that the derivative
$d^{\prime}(\epsilon)=-\overline{\pi}G\overline{x}$ (where
$(\overline{\pi},\overline{x})$ are the dual and primal optimal solutions,
respectively) can be obtained from the well-known matrix inverse fact

$\frac{dM^{-1}(\epsilon)}{d\epsilon}=-M^{-1}(\epsilon)(\frac{dM(\epsilon)}{d\epsilon})M^{-1}(\epsilon).$
(41)

We simply compute $\frac{dM^{-1}(\epsilon)}{d\epsilon}$ for the perturbations
to conditional probabilities, rather than linear perturbations restricted to
the form $A+\epsilon G$, for a matrix $G$ independent of $\epsilon$.

Perturbations for the optimal policy.  Define the $(s,s^{\prime})$ entry of
the perturbed transition matrix, in the direction of a generic observation
$(\tilde{s},\tilde{a},\tilde{s}^{\prime}),$ as
$P^{\epsilon}_{a}(s,s^{\prime})=\frac{(1-\epsilon)p(s,a,s^{\prime})+\epsilon\mathbb{I}\left[s,a,s^{\prime}\right]}{(1-\epsilon)d(s,a)+\epsilon\mathbb{I}\left[s,a\right]}.$

Define the perturbed version of Equation 6:

$\Psi_{D}(P_{\epsilon})=\underset{V}{\operatorname{min}}\left\\{(1-\gamma)\sum_{s}\left((1-\epsilon)\mu_{0}(s)+\epsilon\mathbb{I}\left[\tilde{s}\right]\right)V_{\epsilon}^{*}(s)\colon\left(I-\gamma
P_{a}^{\epsilon}\right){V}-{r}_{a}\geq 0,\quad\forall
a\in\mathcal{A}\right\\}$ (42)

Let $V_{\epsilon}^{*}(s)\in\arg\min\Psi_{D}(P_{\epsilon})$ achieve the optimal
solution of $\Psi_{D}(P_{\epsilon})$.

We will repeatedly use the fact that
${\frac{d}{d\epsilon}\\{(1-\epsilon)p(o)+\epsilon\mathbb{I}\left[o\right]\\}=\epsilon(\mathbb{I}\left[o\right]-p(o)).}$

Then

$\displaystyle\frac{d}{d\epsilon}\Psi_{D}(P_{\epsilon})\Big{|}_{\epsilon=0}$
$\displaystyle=(1-\gamma)\sum_{s}\frac{d}{d\epsilon}\left\\{(1-\epsilon)\mu_{0}(s)+\epsilon\mathbb{I}\left[\tilde{s}\right]\right\\}\Big{|}_{\epsilon=0}V^{*}(s)+(1-\gamma)\sum_{s}\mu_{0}(s)\frac{d}{d\epsilon}\\{V_{\epsilon}^{*}(s)\\}\Big{|}_{\epsilon=0}$
$\displaystyle=(1-\gamma)\sum_{s}(\mathbb{I}\left[s=\tilde{s}\right]-\mu_{0}(s))V^{*}(s)+(1-\gamma)\sum_{s}\mu_{0}(s)\frac{d}{d\epsilon}\\{V_{\epsilon}^{*}(s)\\}\Big{|}_{\epsilon=0}.$

where $V^{*}(s)$ is the unperturbed solution. Note that, letting $P_{a^{*}}$
denote the transition matrix under $a(s)\sim\pi^{*}$, we have that
$V_{\epsilon}^{*}=(I-\gamma P_{a^{*}}^{\epsilon})^{-1}r_{a^{*}}.$

Entrywise, where we overload notation for $a^{*}=a^{*}(s)$ (the optimal action
at $s$, i.e. under active basis), and
$\mathbb{I}\left[\tilde{s},\tilde{a},\tilde{s^{\prime}}\right]=\mathbb{I}\left[(s,a^{*},s^{\prime})=(\tilde{s},\tilde{a},\tilde{s^{\prime}})\right],$
etc:

$\displaystyle\frac{d(I-\gamma
P^{\epsilon}_{a^{*}})_{(s,s^{\prime})}}{d\epsilon}\Big{|}_{\epsilon=0}$
$\displaystyle=-\gamma\;\frac{d}{d\epsilon}\left\\{\frac{(1-\epsilon)p(s,a^{*},s^{\prime})+\epsilon\mathbb{I}\left[\tilde{s},\tilde{a},\tilde{s}^{\prime}\right]}{(1-\epsilon)d(s,a^{*})+\epsilon\mathbb{I}\left[\tilde{s},\tilde{a}\right]}\right\\}\Big{|}_{\epsilon=0}$
$\displaystyle=\gamma\left(\frac{\mathbb{I}\left[\tilde{s},\tilde{a},\tilde{s}^{\prime}\right]-p(s,a^{*},s^{\prime})}{d(s,a^{*})}-p(s,a^{*},s^{\prime})\frac{\mathbb{I}\left[\tilde{s},\tilde{a}\right]-d(s,a^{*})}{d(s,a^{*})^{2}}\right)$
$\displaystyle=\gamma\left(\frac{\mathbb{I}\left[\tilde{s},\tilde{a},\tilde{s}^{\prime}\right]}{d(s,a^{*})}-p(s^{\prime}\mid
s,a^{*})\frac{\mathbb{I}\left[\tilde{s},\tilde{a}\right]}{d(s,a^{*})}\right)$

Hence, applying Equation 41,

$\displaystyle\frac{d}{d\epsilon}\Psi_{D}(P_{\epsilon})\Big{|}_{\epsilon=0}$
$\displaystyle=(1-\gamma)\sum_{s}(\mathbb{I}\left[s=\tilde{s}\right]-\mu_{0}(s))V^{*}(s)+(1-\gamma)\sum_{s}\mu_{0}(s)\frac{d}{d\epsilon}\\{V_{\epsilon}^{*}(s)\\}\Big{|}_{\epsilon=0}$
$\displaystyle=(1-\gamma)\sum_{s}(\mathbb{I}\left[s=\tilde{s}\right]-\mu_{0}(s))V^{*}(s)$
$\displaystyle\qquad+(1-\gamma){\mu_{0}}(I-\gamma
P_{a^{*}})^{-1}\left(\frac{d(I-\gamma
P^{\epsilon}_{a^{*}})}{d\epsilon}\Big{|}_{\epsilon=0}\right)(I-\gamma
P_{a^{*}})^{-1}{r}.$

Note that, as in Freund [1985], $(1-\gamma){\mu_{0}}(I-\gamma
P_{a^{*}})^{-1}=\mu^{*}$ is the dual (unperturbed) optimal solution of Example
3 and $(I-\gamma P_{a^{*}})^{-1}{r}$ is the primal (unperturbed) optimal
solution $V^{*}$. Therefore,

$\displaystyle\frac{d}{d\epsilon}\Psi_{D}(P_{\epsilon})\Big{|}_{\epsilon=0}$
(43) $\displaystyle=(1-\gamma)(V^{*}(\tilde{s})-\sum_{s}\mu_{0}(s)V^{*}(s))$
$\displaystyle\qquad-\gamma\sum_{s}\mu^{*}(s,a)\sum_{s^{\prime}}V^{*}(s^{\prime})\left(\frac{\mathbb{I}\left[\tilde{s},\tilde{a},\tilde{s}^{\prime}\right]-p(s,a^{*},s^{\prime})}{d(s,a^{*})}-p(s,a^{*},s^{\prime})\frac{\mathbb{I}\left[\tilde{s},\tilde{a}\right]-d(s,a^{*})}{d(s,a^{*})^{2}}\right)$
(44)
$\displaystyle=(1-\gamma)(V^{*}(\tilde{s})-\sum_{s}\mu_{0}(s)V^{*}(s))-\gamma\frac{\mu^{*}(\tilde{s},\tilde{a})}{d(\tilde{s},\tilde{a})}\left(V(\tilde{s}^{\prime})-(P^{\pi}V^{*})(\tilde{s})\right)$

where
$(P^{\pi}V^{*})(\tilde{s})=\sum_{s^{\prime}}p(s^{\prime}\mid\tilde{s},a)V^{*}(\tilde{s})$
is the Bellman operator evaluated starting from the state $\tilde{s}$.

Further note that
$\gamma(P^{\pi}V^{*})(\tilde{s})=(V(\tilde{s})-r(\tilde{s},\tilde{a}))$, (i.e.
a Bellman residual term), and
$\Psi_{D}(P)={(1-\gamma)\sum_{s}\mu_{0}(s)V^{*}(s))}.$

Hence

$\frac{d}{d\epsilon}\Psi_{D}(P_{\epsilon})\Big{|}_{\epsilon=0}=(1-\gamma)V^{*}(\tilde{s})-\frac{\mu^{*}(\tilde{s},\tilde{a})}{d(\tilde{s},\tilde{a})}\left(r(\tilde{s},\tilde{a})+\gamma
V(\tilde{s}^{\prime})-V(\tilde{s})\right)-\Psi_{D}(P)$ (45)

∎

###### Proof of Proposition 4.

We analyze the error introduced by finite differencing. Consider the dual
formulation of Equation 6, and the perturbed formulation of Equation 42. Fix
$\epsilon$ small enough such that the active basis remains the same under $P$
and $P_{\epsilon}$; let $P^{\pi}$ denote the submatrix of the transition
matrix $P$ (stacked for the optimization problem) corresponding to the active
basis (equivalently, optimal policy).

Define the entries of $G_{\epsilon}$ such that $P_{\epsilon}=P+\epsilon
G_{\epsilon},$ i.e. (analogous to Lemma 2)

$(G_{\epsilon}^{a})_{s,s^{\prime}}=\frac{\mathbb{I}\left[\tilde{s},\tilde{a},\tilde{s}^{\prime}\right]}{d_{\epsilon}(s,a)}-\frac{\mathbb{I}\left[\tilde{s},\tilde{a}\right]P(s^{\prime}\mid
s,a)}{d_{\epsilon}(s,a)}$

Let $V_{\epsilon}^{\pi}(s),V^{\pi}(s)$ be the active subvector of the optimal
solution achieving the optimal values of $\Psi(P_{\epsilon}),\Psi(P)$,
respectively.

Then

$\displaystyle\epsilon^{-1}(\Psi_{D}(P_{\epsilon})-\Psi_{D}(P))$
$\displaystyle\textstyle=\epsilon^{-1}(1-\gamma)\left(\sum_{s}\left((1-\epsilon)\mu_{0}(s)+\epsilon\mathbb{I}\left[\tilde{s}\right]\right)V_{\epsilon}^{\pi}(s)-\sum_{s}\mu_{0}(s)V^{\pi}(s)\right)$
$\displaystyle\textstyle=\epsilon^{-1}(1-\gamma)\left(\sum_{s}\epsilon\left(\mathbb{I}\left[\tilde{s}\right]-\mu_{0}(s)\right)V_{\epsilon}^{\pi}(s)+\sum_{s}\mu_{0}(s)(V_{\epsilon}^{\pi}(s)-V^{\pi}(s))\right)$
$\displaystyle\textstyle=(1-\gamma)\left(\sum_{s}\left(\mathbb{I}\left[\tilde{s}\right]-\mu_{0}(s)\right)V_{\epsilon}^{\pi}(s)+\epsilon^{-1}\sum_{s}\mu_{0}(s)(V_{\epsilon}^{\pi}(s)-V^{\pi}(s))\right)$
(46)

The optimal solutions $V_{\epsilon}^{\pi},V^{\pi}$ differ in the perturbation
matrix; inverting the basic feasible solution, where invertibility is assured
by nondegeneracy, gives that:

$\displaystyle V_{\epsilon}^{\pi}=(I-\gamma
P^{\pi}_{\epsilon})^{-1}r_{\pi},\qquad V^{\pi}=(I-\gamma
P^{\pi})^{-1}r_{\pi}.$

We analyze $V_{\epsilon}^{\pi}-V^{\pi}=((I-\gamma
P^{\pi}_{\epsilon})^{-1}-(I-\gamma P^{\pi})^{-1})r_{\pi}$ via a general result
about perturbed matrices. Note that for two matrices $A,B$, clearly
$(A+B)^{-1}=A^{-1}-A^{-1}B(A+B)^{-1}.$ Applying this recursively we obtain the
power series expansion

$(A+B)^{-1}=A^{-1}-A^{-1}\sum_{k=0}^{\infty}(B)^{k}(A^{-1})^{k}.$

Applying this result to our setting with $A=(I-\gamma P_{\epsilon}^{\pi})$ and
$B=-\gamma\epsilon G_{\epsilon}$, we obtain

$V_{\epsilon}^{\pi}-V^{\pi}=\epsilon\gamma(I-\gamma
P^{\pi})^{-1}G_{\epsilon}^{\pi}(I-\gamma P^{\pi})^{-1}r_{\pi}+O(\epsilon^{2})$
(47)

The statement follows by substituting in the above expression, evaluating
based on definition of the entries of $G_{\epsilon}^{\pi}$, and substituting
the above and Equation 46 into Equation 44 of the proof of Proposition 3 and