---
abstract: 'We develop a method for predicting the performance of reinforcement learning and bandit algorithms, given historical data that may have been generated by a different algorithm. Our estimator has the property that its prediction converges in probability to the true performance of a counterfactual algorithm at the fast $\sqrt{N}$ rate, as the sample size $N$ increases. We also show a correct way to estimate the variance of our prediction, thus allowing the analyst to quantify the uncertainty in the prediction. These properties hold even when the analyst does not know which among a large number of potentially important state variables are really important. These theoretical guarantees make our estimator safe to use. We finally apply it to improve advertisement design by a major advertisement company. We find that our method produces smaller mean squared errors than state-of-the-art methods.'
bibliography:
- 'reference.bib'
---

Introduction
============

Interactive reinforcement learning (RL) and bandit systems (e.g. personalized education and medicine, ad/news/recommendation/search platforms) produce log data valuable for evaluating and redesigning the systems. For example, the logs of a news recommendation system record which news article was presented and whether the user read it, giving the system designer a chance to make its recommendation more relevant.

Exploiting log data is, however, more difficult than conventional supervised machine learning: the result of each log is only observed for the action chosen by the system (e.g. the presented news) but not for all the other actions the system could have taken. Moreover, the log entries are biased in that the logs over-represent actions favored by the system.

A potential solution to this problem is an A/B test that compares the performance of counterfactual systems. However, A/B testing counterfactual systems is often technically or managerially infeasible, since deploying a new policy is time- and money-consuming, and entails a risk of failure.

This leads us to the problem of *counterfactual (off-policy, offline) evaluation and learning*, where one aims to use batch data collected by a logging policy to estimate the value of a counterfactual policy or algorithm without employing it. Such evaluation allows us to compare the performance of counterfactual policies to decide which policy should be deployed in the field. This alternative approach thus solves the above problem with the naive A/B test approach. Key prior studies include @li2010contextual [@Strehl2010; @li2011unbiased; @Li2012; @Swaminathan2015; @Swaminathan2015b; @wang2016optimal; @swaminathan2017off; @Dimakopoulou2017; @narita2018efficient] for bandit algorithms, and @precup2000eligibility [@precup2001temporal; @bottou2013counterfactual; @thomas2015confidence; @Jiang16; @munos2016safe; @Thomas16; @gu2017interpolated; @liu2018representation; @Farajtabar2018MoreRD; @hanna2018importance; @Irpan2019OffPolicyEV; @kallus2019intrinsically] for RL algorithms.

**Method.** For off-policy evaluation with log data of RL feedback, this paper develops and empirically implements a novel technique with desirable theoretical properties. To do so, we consider a class of RL algorithms, including contextual bandit algorithms as important special cases. This class includes most of the widely-used algorithms such as (deep) Q-learning, Actor Critique, contextual $\epsilon$-Greedy and Thompson Sampling, as well as their non-contextual analogs and random A/B testing. We allow the logging policy to be an unknown function of numerous potentially important state variables. This feature is salient in real-world applications. We also allow the evaluation target policy to be degenerate, again a key feature of real-life situations.

             ![image](figure1.png){width="15cm"}

*Notes*: This figure shows root mean squared error of the predicted CTRs of the evaluation policy compared to its actual CTR. Root mean squared error is normalized by the actual CTRs for confidentiality reasons. We obtain these estimates by the inverse probability weighting estimator, Doubly Robust estimator, and Double Machine Learning estimator (our proposal) using the true logging policy (propensity score) or estimated one. These objects are defined and analyzed in Section 3. Bars indicate 95% confidence intervals based on our asymptotic variance estimators.

We consider an offline estimator for the expected reward from a counterfactual policy. Our estimator integrates a well-known doubly robust (DR) estimator (@rotnitzky1995semiparametric and modern studies cited above) with “double/debiased machine learning" developed in econometrics and statistics [@Chernozhukov2018; @Chernozhukov2018LR]. Building upon these prior studies, we show the following result:

> **Theoretical Result.** Our estimator is “$\sqrt{N}$-consistent” in the sense that its prediction converges in probability to the true performance of a counterfactual policy at the $\sqrt{N}$ rate as the sample size $N$ increases. Our estimator is also “asymptotically normal." We provide a consistent estimator of its asymptotic variance, thus allowing for measuring statistical uncertainty in our prediction.

For many special cases, including all contextual bandit algorithms, our estimator is shown to have lowest variance in a wide class of estimators, achieving variance reduction relative to standard estimators. These theoretical properties suggest that our estimator is safe to use. In contrast, existing estimators lack or have not been shown to have these theoretical guarantees.

**Application.** We empirically apply our estimator to evaluate and optimize the design of online advertisement formats. Our application is based on proprietary data provided by CyberAgent Inc., the second largest Japanese advertisement company with about 5 billion USD market capitalization (as of February 2020). This company uses randomly chosen bandit algorithms to determine the visual design of advertisements assigned to users. This A/B test of randomly choosing an algorithm produces logged data and the ground truth for the performance of alternative algorithms.

We use this data to examine the performance of our proposed method. Specifically, we use the log data from an algorithm to predict the the click through rates (CTR) of another algorithm, and then assess the accuracy of our prediction by comparing it with the ground truth. This exercise shows the following:

> **Empirical Result.** Our estimator produces smaller mean squared errors than widely-used benchmark methods in the spirit of @Jiang16 and @Thomas16.

This result is reported in Figure \[fig:result2\], where the mean squared errors using our estimator (Double Machine Learning; colored red) are lower than others using existing estimators (Inverse Probability Weighting or Doubly Robust). This improvement is statistically significant at the 5% level. Importantly, this result holds regardless of whether we know the data-generating logging policy or not; these two scenarios correspond to the figure’s two panels, respectively. This finding shows that our estimator can substantially reduce bias and uncertainty we face in real-world decision-making.

This improved performance motivates us to use our estimator to optimize the advertisement design for maximizing the CTR. We estimate how much the CTR would be improved by a counterfactual policy of choosing the best action (advertisement) for each context (user characteristics). This exercise produces the following bottomline: Our estimator predicts the hypothetical policy to statistically significantly improve the CTR by 30% (compared to the logging policy) in one of the three campaigns we analyze. Our approach thus generates valuable managerial conclusions.

Setup
=====

Data Generating Process
-----------------------

We consider historical data of trajectories that follow a Markov Decision Process (MDP) as a mathematical description of RL and bandit algorithms. Specifically, we closely follow the setups of @Jiang16, @Thomas16 and @Farajtabar2018MoreRD.

An MDP is given by ${\cal M}=\langle {\cal S},{\cal A},P_{S_0},P_S,P_R\rangle$, where $\cal{S}$ is the state space, $\cal{A}$ is the action space, $P_{S_0}: {\cal S}\rightarrow [0,1]$ is the initial state distribution, $P_S: {\cal S}\times {\cal A}\rightarrow \Delta({\cal S})$ is the transition function with $P_S(s'|s,a)$ being the probability of seeing state $s'$ after taking action $a$ given state $s$, and $P_R:{\cal S}\times {\cal A}\times \mathbb{R}\rightarrow [0,1]$ be the conditional distribution of the immediate reward with $P_R(\cdot|s,a)$ being the immediate reward distribution conditional on the state and action being $(s,a)$. Given $P_R$, we define the mean reward function $\mu:{\cal S}\times {\cal A}\rightarrow \mathbb{R}$ as $\mu(s,a)=\int r dP_R(r|s,a)$, and define the reward variance function $\sigma_R^2:{\cal S}\times {\cal A}\rightarrow \mathbb{R}_{+}$ as $\sigma^2_R(s,a)=\int (r-\mu(s,a))^2dP_R(r|s,a)$. For simplicity, we assume that the state space ${\cal S}$ and the action space ${\cal A}$ are finite.

We call a function $\pi: {\cal S}\rightarrow \Delta({\cal A})$ a [*policy*]{}, which assigns each state $s\in {\cal S}$ a distribution over actions, where $\pi(a|s)$ is the probability of taking action $a$ when the state is $s$. Let $H=(S_0,A_0,R_0,...,S_T,A_T,R_T)$ be a [*trajectory*]{}, where $S_t$, $A_t$ and $R_t$ are the state, the action, and the reward in step $t$, respectively, and $T$ is the number of steps. We say that a trajectory $H$ is generated by a policy $\pi$, or $H\sim \pi$ in short if $H$ is generated by the following process:

-   The initial state $S_0$ is drawn from the initial distribution $P_{S_0}$. Given $S_0$, the action $A_0$ is randomly chosen based on $\pi(\cdot|S_0)$. The reward $R_0$ is drawn from the conditional reward distribution $P_R(\cdot|S_0,A_0)$.

-   For $t=1,...,T$, the state $S_t$ is determined based on the transition function $P_S(\cdot|S_{t-1},A_{t-1})$. Given $S_t$, the action $A_t$ is randomly chosen based on $\pi(\cdot|S_t)$. The reward $R_t$ is drawn from the conditional reward distribution $P_R(\cdot|S_t,A_t)$.

Suppose that we observe [*historical data*]{} $\{H_i\}_{i=1}^N$ where trajectories are independently generated by a fixed [*behavior policy*]{} $\pi_b$, i.e., $H_i\sim \pi_b$ independently across $i$. The historical data is a collection of iid trajectories. Importantly, we allow the components of the data generating process ${\cal M}$ and $\pi_b$ to vary with the sample size $N$. Specifically, let ${\cal M}_N$ and $\pi_{bN}$ be the MDP and the behavior policy, respectively, when the sample size is $N$, and let $\mathbb{P}_N$ denote the resulting probability distribution of $H_i$. $\mathbb{P}_N$ is allowed to vary with $N$ in a way that the functions $P_{S_0N}$, $P_{RN}$ and $\pi_{bN}$ are high dimensional relative to sample size $N$ even when $N$ is large. In some RL problems, for example, there are a large number of possible states. To capture the feature that the number of states $|{\cal S}_N|$ is potentially large relative to sample size $N$, we may consider a sequence of $\mathbb{P}_N$ such that $|{\cal S}_N|$ is increasing with $N$. For the sake of notational simplicity, we make implicit the dependence of ${\cal M}$ and $\pi_{b}$ on $N$.

We assume that we know the state space ${\cal S}$ and the action space ${\cal A}$ but know none of the functions $P_{S_0}$, $P_S$ and $P_R$. In some environments, we know the function $\pi_b$ or observe the probability vector $(\pi_b(a|S_{it}))_{a\in {\cal A}, t=1,...,T}$ for every trajectory $i$ in the historical data. Our approach is usable regardless of the availability of such knowledge on the behavior policy.

Examples
--------

This data generating process allows for many popular RL and bandit algorithms, as the following examples illustrate.

\[ex:deepQ\] In each round $t$, given state $s_t$, a Q Learning algorithm picks the best action based on the estimated Q-value of each actions, $Q(s, a)$, which estimates the expected cumulative reward from taking action $a$ (following the state and the policy). Choice probabilities can be determined with an $\epsilon$-Greedy or soft-max rule, for instance. In the case where the soft-max rule is employed, the probability of taking each action is as follows: $$\pi(a|s_t) = \frac{ \exp ( Q(s_t,a)) }{\sum_{a' \in A} \exp ( Q(s_t,a')) }.$$ Deep Q Learning algorithms estimate Q-value functions through deep learning methods.

\[ex:AC\] An Actor Critic is a hybrid method of value-based approach such as Q-learning and policy-based method such as REINFORCE. This algorithm has two components called Actor and Critic. Critic estimates the value function and Actor updates the policy using the value of Critic. In each round t, we pick the best action according to the value of Actor with some probability. As in Deep $Q$ Learning algorithms, we can use $\epsilon$-Greedy and soft-max for determining an action.

Contextual bandit algorithms are also important examples. When $T=0$, a trajectory takes the form of $H=(S_0,A_0,R_0)$. Regarding $S_0$ as a context, it is possible to consider $\{H_i\}_{i=1}^N$ as batch data generated by a contextual bandit algorithm. In additional examples below, the algorithms use past data to estimate the mean reward function $\mu$ and the reward variance function $\sigma^2_R$. Let $\hat\mu$ and $\hat \sigma^2_R$ denote any given estimators of $\mu$ and $\sigma^2_R$, respectively.

\[ex:e-greedy\] When the context is $s_0$, we choose the best action based on $\hat \mu(s_0,a)$ with probability $1-\epsilon$ and choose an action uniformly at random with probability $\epsilon$: $$\pi(a|s_0)= \begin{cases}
    1-\epsilon + \frac{\epsilon}{|{\cal A}|}& \ \ \ \text{if $a={\mathop{\rm argmax}\limits}_{a'\in {\cal A}}\hat \mu(s_0,a')$}\\
    \frac{\epsilon}{|{\cal A}|} & \ \ \ \text{otherwise}.
    \end{cases}$$

\[ex:Thompson\] When the context is $s_0$, we sample the potential reward $r_0(a)$ from the normal distribution $N(\hat \mu(s_0,a),\hat \sigma^2_R(s_0,a))$ for each action, and choose the action with the highest sampled potential reward, ${\mathop{\rm argmax}\limits}_{a'\in {\cal A}}r_0(a')$. As a result, this algorithm chooses actions with the following probabilities: $$\pi(a|s_0) = \Pr\{a={\mathop{\rm argmax}\limits}_{a'\in {\cal A}}r_0(a')\},$$ where $(r_0(a))_{a\in{\cal A}}\sim N(\hat \mu(s_0),\hat \Sigma(s_0))$, $\hat \mu(s_0)=(\hat \mu(s_0,a))_{a\in{\cal A}}$, and $\hat \Sigma(s_0)$ is the diagonal matrix whose diagonal entries are $(\sigma^2_R(s_0,a))_{a\in{\cal A}}$.

Prediction Target
-----------------

With the historical data $\{H_i\}_{i=1}^N$, we are interested in estimating the discounted value of the [*evaluation policy*]{} $\pi_e$, which might be different from $\pi_b$: $$V^{\pi_e}\coloneqq\mathbb{E}_{H\sim \pi_e}[\sum_{t=0}^T \gamma^t R_t],$$ where $\gamma\in [0,1]$ is the discount factor.

Estimator and Its Properties {#section:estimation}
============================

The estimation of $V^{\pi_e}$ involves estimation of the behavior policy $\pi_b$ (if unknown), the transition function $P_S$, and the mean reward function $\mu$. These functions may be high dimensional and complex. To handle the issue, we use the double/debiased machine learning (DML) method developed by @Chernozhukov2018.

Before presenting our proposed estimator, we introduce some notations. $H_t^{s,a}=(S_0,A_0,...,S_t,A_t)$ is a trajectory of the state and action up to step $t$. $\rho_t^{\pi_e}: ({\cal S}\times {\cal A})^t\rightarrow \mathbb{R}_+$ is an importance weight function with $$\rho_t^{\pi_e}(H_t^{s,a})\coloneqq \prod_{t'=0}^t \frac{\pi_e(A_{t'}|S_{t'})}{\pi_b(A_{t'}|S_{t'})}$$ being the probability of the first $t$ steps of $H$ under the evaluation policy $\pi_e$ divided by its probability under the behavior policy $\pi_b$. Viewing $\rho_t^{\pi_e}$ as a function of $\pi_b$, define $\rho_t^{\pi_e}(H_t^{s,a};\tilde \pi_b)$ as the value of $\rho_t^{\pi_e}(H_t^{s,a})$ defined as above with the true behavior policy $\pi_b$ replaced with a candidate function $\tilde \pi_b$, i.e., $$\rho_t^{\pi_e}(H_t^{s,a};\tilde \pi_b)\coloneqq \prod_{t'=0}^t\frac{\pi_e(A_{t'}|S_{t'})}{\tilde \pi_b(A_{t'}|S_{t'})}.$$ We can think of $\rho_t^{\pi_e}(H_t^{s,a};\tilde \pi_b)$ as the estimated importance weight function when we use $\tilde \pi_b$ as the estimate of $\pi_b$. By definition, $\rho_t^{\pi_e}(H_t^{s,a};\pi_b)=\rho_t^{\pi_e}(H_t^{s,a})$, where the left-hand side is $\rho_t^{\pi_e}(H_t^{s,a};\tilde \pi_b)$ evaluated at the true $\pi_b$ and the right-hand side is the true importance weight function.

Finally, let $q_t^{\pi_e}:{\cal S}\times {\cal A}\rightarrow \mathbb{R}$ be the action-value function under policy $\pi_e$ at step $t$, where $q_t^{\pi_e}(s,a)\coloneqq\mathbb{E}_{H\sim \pi_e}[\sum_{t'=t}^T\gamma^{t'-t}R_{t'}|S_t=s,A_t=a]$. Using the transition function $P_S$ and the mean reward function $\mu$, $q_t^{\pi_e}$ can be obtained recursively: let $$\begin{aligned}
    q_T^{\pi_e}(s,a) &= \mu(s,a) \label{eq:Q_T}\end{aligned}$$ and for $t=0,...,T-1$, $$\begin{aligned}
    &q_t^{\pi_e}(s,a)\nonumber\\
    =&\mu(s,a)+\gamma\sum_{(s',a')}P_S(s'|s,a)\pi_e(a'|s')q_{t+1}^{\pi_e}(s',a'). \label{eq:Q_t}\end{aligned}$$ Our estimator is based on the following expression of $V^{\pi_e}$ [@Thomas16]: $$\begin{aligned}
V^{\pi_e}=&~\mathbb{E}_{H\sim \pi_b}[\psi(H;\pi_b,\{q_t^{\pi_e}\}_{t=0}^T)] \label{eq:value},\end{aligned}$$ where for any candidate tuple $\tilde \eta=(\tilde \pi_b,\{\tilde q_t^{\pi_e}\}_{t=0}^T)$, we define $$\begin{aligned}
&\psi(H;\tilde \eta)\\
=&~\sum_{t=0}^T \gamma^t \{\rho_t^{\pi_e}(H_{t}^{s,a};\tilde \pi_b)(R_{t}-\tilde q_t^{\pi_e}(S_{t},A_{t}))\nonumber\\
&~~~~~~~+\rho_{t-1}^{\pi_e}(H_{t-1}^{s,a};\tilde \pi_b)\sum_{a\in{\cal A}}\pi_e(a|S_{t})\tilde q_t^{\pi_e}(S_{t},a)\},\end{aligned}$$ where we set $\rho_{-1}^{\pi_e}=0$. We construct our estimator as follows.

1.  Take a $K$-fold random partition $(I_k)_{k=1}^K$ of trajectory indices $\{1,...,N\}$ such that the size of each fold $I_k$ is $n=N/K$. Also, for each $k=1,...,K$, define $I_k^c\coloneqq \{1,...,N\}\setminus I_k$.

2.  For each $k=1,...,K$, construct estimators $\hat \pi_{b,k}$ (if $\pi_b$ is unknown), $\hat\mu_{k}$ and $\hat P_{S,k}$ of $\pi_b$, $\mu$ and $P_S$ using the subset of data $\{H_i\}_{i\in I_k^c}$. We then construct estimator $\{\hat q^{\pi_e}_{t,k}\}_{t=0}^T$ of $\{q^{\pi_e}_t\}_{t=0}^T$ by plugging $\hat\mu_{k}$ and $\hat P_{S,k}$ into the recursive formulation (\[eq:Q\_T\]) and (\[eq:Q\_t\]).

3.  Given $\hat\eta_k=(\hat\pi_{b,k},\{\hat q_{t,k}^{\pi_e}\}_{t=0}^T)$, $k=1,...,K$, the DML estimator $\hat V^{\pi_e}_{{\rm DML}}$ is given by $$\begin{aligned}
        \hat V^{\pi_e}_{{\rm DML}}
        =&\frac{1}{K}\sum_{k=1}^K \frac{1}{n}\sum_{i\in I_k}\psi(H_{i};\hat \eta_k)\\
        =& \frac{1}{K}\sum_{k=1}^K \frac{1}{n}\sum_{i\in I_k}\sum_{t=0}^T \gamma^t \{\rho_t^{\pi_e}(H_{it}^{s,a};\hat \pi_{b,k})(R_{it}-\\
        &\hat q_{t,k}^{\pi_e}(S_{it},A_{it}))+\rho_{t-1}^{\pi_e}(H_{it-1}^{s,a};\hat \pi_{b,k})\\
        &\sum_{a\in{\cal A}}\pi_e(a|S_{it})\hat q_{t,k}^{\pi_e}(S_{it},a)\}.
        \end{aligned}$$

Possible estimation methods for $\pi_b$, $\mu$ and $P_S$ in Step 2 are (i) classical nonparametric methods such as kernel and series estimation, (ii) off-the-shelf machine learning methods such as random forests, lasso, neural nets, and boosted regression trees, and (iii) existing methods developed in the off-policy policy evaluation literature such as representation balancing MDPs [@liu2018representation]. These methods, especially (ii) and (iii), are usable even when the analyst does not know the relevant state variables and there are a large number of potentially important state variables. This DML estimator differs from the DR estimator developed by @Jiang16 and @Thomas16 in that we use the cross-fitting procedure, as explained next.

[**A. Cross-Fitting.**]{} The above method uses a sample-splitting procedure called [*cross-fitting*]{}, where we split the data into $K$ folds, take the sample analogue of (\[eq:value\]) using one of the folds ($I_k$) with $\pi_b$ and $\{ q_t^{\pi_e}\}_{t=0}^T$ estimated from the remaining folds ($I_k^c$) plugged in, and average the estimates over the $K$ folds to produce a single estimator. Cross-fitting has two advantages. First, if we used instead the whole sample both for estimating $\pi_b$ and $\{ q_t^{\pi_e}\}_{t=0}^T$ and for computing the final estimate of $V^{\pi_e}$ (the “full-data” variant of the DR estimator of @Thomas16), substantial bias might arise due to overfitting [@Chernozhukov2018; @newey2018crossfitting]. Cross-fitting removes the potential bias by making the estimate of $V^{\pi_e}$ independent of the estimates of $\pi_b$ and $\{ q_t^{\pi_e}\}_{t=0}^T$. Thanks to this, the DML estimator has properties such as $\sqrt{N}$-consistency and asymptotic normality under milder condtions than those necessary without sample splitting.

Second, a standard sample splitting procedure uses a half of the data to construct estimates of $\pi_b$ and $\{ q_t^{\pi_e}\}_{t=0}^T$ and the other half to compute the estimate of $V^{\pi_e}$ (the DR estimator of @Jiang16 and the “half-data” variant of the DR estimator of @Thomas16).[^1] In contrast, cross-fitting swaps the roles of the main fold ($I_k$) and the rest ($I_k^c$) so that all trajectories are used for the final estimate, which enables us to make efficient use of data.

[**B. Neyman Orthogonality.**]{} There is another key ingredient important for DML to have desirable properties. The DML estimator is constructed by plugging in the estimates of $\pi_b$ and $\{ q_t^{\pi_e}\}_{t=0}^T$, which may be severely biased due to regularization if they are estimated with machine learning methods. However, the DML estimator is expected to be robust to the bias, since $\psi$ satisfies the [*Neyman orthogonality*]{} condition [@Chernozhukov2018]. The condition requires that for any candidate tuple $\tilde \eta=(\tilde \pi_b,\{\tilde q_t^{\pi_e}\}_{t=0}^T)$, $$\left. \frac{\partial \mathbb{E}_{H\sim \pi_b}[\psi(H;\eta+r(\tilde \eta-\eta))]}{\partial r}\right|_{r=0}=0,$$ where $\eta=(\pi_b,\{q_t^{\pi_e}\}_{t=0}^T)$ is the tuple of the true functions (see the supplementary material for the proof that DML satisfies this). Intuitively, the Neyman orthogonality condition means that the right-hand side of (\[eq:value\]) is locally insensitive to the value of $\pi_b$ and $\{q_t^{\pi_e}\}_{t=0}^T$. As a result, plugging in noisy estimates of $\pi_b$ and $\{q_t^{\pi_e}\}_{t=0}^T$ does not strongly affect the estimate of $V^{\pi_e}$.

In contrast, the well-known Inverse Probability Weighting (IPW) estimator is based on the following expression of $V^{\pi_e}$: $V^{\pi_e}=\mathbb{E}_{H\sim \pi_b}[\psi_{IPW}(H;\pi_b)]$, where for any candidate $\tilde \pi_b$, we define $\psi_{IPW}(H;\tilde \pi_b)=\sum_{t=0}^T \gamma^t \rho_t^{\pi_e}(H_{t}^{s,a};\tilde \pi_b)R_{t}$. The Neyman orthogonality condition does not hold for IPW: for some $\tilde \pi_b\neq \pi_b$, $$\left. \frac{\partial \mathbb{E}_{H\sim \pi_b}[\psi_{IPW}(H;\pi_b+r(\tilde \pi_b-\pi_b))]}{\partial r}\right|_{r=0}\neq0.$$ Therefore, IPW is not robust to bias in the estimate of $\pi_b$.

$\sqrt{N}$-consistency and Asymptotic Normality
-----------------------------------------------

Let $$\sigma^2 =\mathbb{E}_{H\sim \pi_b}[(\psi(H;\eta)-V^{\pi_e})^2]$$ be the variance of $\psi(H;\eta)$. To derive the properties of the DML estimator, we make the following assumption. Recall that $\mathbb{P}_N$, the probability distribution of $H_i$, is allowed to vary with the sample size $N$, and that any characteristics of the distribution may implicitly depend on $N$.

\[as:regularity\]

1.  \[as:variance\] There exists a constant $c_0>0$ such that $c_0\le \sigma^2 <\infty$ for all $N$.

2.  \[as:nuisance-est\]The estimator $\hat \eta_k=(\hat \pi_{b,k},\{\hat q_{t,k}^{\pi_e}\}_{t=0}^T)$ belongs to a set ${\cal T}_N$ with probability approaching one, where ${\cal T}_N$ contains the true $\eta=(\pi_b,\{q_{t}^{\pi_e}\}_{t=0}^T)$ and satisfies the following:

    (i) \[as:moment\] There exist constants $q>2$ and $c_1>0$ such that $\sup_{\tilde \eta\in {\cal T}_N}(\mathbb{E}_{H\sim \pi_b}[(\psi(H;\tilde\eta)-V^{\pi_e})^q])^{1/q}\le c_1$ for all $N$.

    (ii) \[as:dif\] $
                \sup_{\tilde \eta\in {\cal T}_N}(\mathbb{E}_{H\sim \pi_b}[(\psi(H;\tilde \eta)-\psi(H;\eta))^2])^{1/2} = o(1).
                $

    (iii) \[as:second\_d\] $
                \sup_{r\in(0,1),\tilde \eta\in {\cal T}_N} \left|\dfrac{\partial^2 \mathbb{E}_{H\sim \pi_b}[\psi(H;\eta+r(\tilde \eta-\eta))]}{\partial r^2}\right|\\=o(1/\sqrt{N}).
                $

Assumption \[as:regularity\] \[as:variance\] assumes that the variance of $\psi(H;\eta)$ is nonzero and finite. Assumption \[as:regularity\] \[as:nuisance-est\] states that the estimator $(\hat \pi_{b,k},\{\hat q_{t,k}^{\pi_e}\}_{t=0}^T)$ belongs to the set ${\cal T}_N$, a shrinking neighborhood of the true $(\pi_b,\{q_{t}^{\pi_e}\}_{t=0}^T)$, with probability approaching one. It requires that $(\hat \pi_{b,k},\{\hat q_{t,k}^{\pi_e}\}_{t=0}^T)$ converge to $(\pi_b,\{q_{t}^{\pi_e}\}_{t=0}^T)$ at a sufficiently fast rate so that the three rate conditions in Assumption \[as:regularity\] \[as:nuisance-est\] are satisfied.

The following proposition establishes the $\sqrt{N}$-consistency and asymptotic normality of $\hat V^{\pi_e}_{{\rm DML}}$ and provides a consistent estimator for the asymptotic variance.

\[prop:step-DR\] Suppose that Assumption \[as:regularity\] holds. Then, $$\sqrt{N}\sigma^{-1}(\hat V^{\pi_e}_{{\rm DML}}-V^{\pi_e})\rightsquigarrow N(0,1),$$ where the symbol $\rightsquigarrow$ denotes convergence in distribution. Moreover, let $$\begin{aligned}
\hat \sigma^2&=\frac{1}{K}\sum_{k=1}^K \frac{1}{n}\sum_{i\in I_k}(\psi(H_i;\hat \eta_k)-\hat V^{\pi_e}_{{\rm DML}})^2\end{aligned}$$ be a sample analogue of $\sigma^2$. We then have $$\hat \sigma^2=\sigma^2 +o_p(1)$$ $$\sqrt{N}\hat \sigma^{-1}(\hat V^{\pi_e}_{{\rm DML}}-V^{\pi_e})\rightsquigarrow N(0,1).$$

The above convergence results hold under any sequence of probability distributions $\{\mathbb{P}_N\}_{N\ge 1}$ as long as Assumption \[as:regularity\] holds. Therefore, our approach is usable, for example, in the case where there are a growing number of possible states, that is, $|{\cal S}|$ is increasing with $N$.

A Special Case: Contextual Bandits
----------------------------------

Suppose that the algorithms are contextual bandits, as in Examples \[ex:e-greedy\] and \[ex:Thompson\]. When $T=0$, the DML estimator becomes $$\begin{aligned}
    \hat V^{\pi_e}_{DML} =& \frac{1}{K}\sum_{k=1}^K \frac{1}{n}\sum_{i\in I_k} \{\frac{\pi_e(A_{i0}|S_{i0})}{\hat\pi_{b,k}(A_{i0}|S_{i0})} (R_{i0}-\hat\mu_k(S_{i0},A_{i0}))\\
    &~~~~~~~~~~~~~~~~~~+\sum_{a\in{\cal A}}\pi_e(a|S_{i0})\hat\mu_k(S_{i0},a)\},\end{aligned}$$ where $(\hat \pi_{b,k},\hat \mu_k)$ is the estimator of $(\pi_b,\mu)$ using the subset of data $\{H_i\}_{i\in I_k^c}$. This estimator is the same as the DR estimator of @Dudik2014 except that we use the cross fitting procedure. Proposition \[prop:step-DR\] has the following implication for the contextual bandit case.

Suppose that $T=0$ and that Assumption \[as:regularity\] holds. Then, $$\sqrt{N}\sigma^{-1}_{CB}(\hat V^{\pi_e}_{DML}-V^{\pi_e})\rightsquigarrow N(0,1),$$ where $$\begin{aligned}
        \sigma^2_{CB}&=\mathbb{E}_{S_0\sim P_{S_0}}[\sum_{a\in{\cal A}}\frac{\pi_e(a|S_0)^2}{\pi_b(a|S_0)}\sigma^2_R(S_0,a)\\
        &~~~~~~+(\sum_{a\in{\cal A}}\pi_e(a|S_0)\mu(S_0,a)-V^{\pi_e})^2].
    \end{aligned}$$

The variance expression coincides with the “semiparametric efficiency bound" obtained by @narita2018efficient, where the semiparametric efficiency bound is the smallest possible asymptotic variance among all consistent and asymptotically normal estimators. Hence $\hat V^\pi_{DML}$ is the lowest variance estimator at a given $\mathbb{P}_N=\mathbb{P}$.

Experiment
==========

                                       ![image](sample_creatives.png){width="10cm"}

  ------------- ------------------- -------------------
                                    
  Method           relative-RMSE       relative-RMSE
  IPW            0.72531 (0.03077)   0.72907 (0.03061)
  DR             0.65470 (0.02900)   0.65056 (0.02880)
  DML            0.56623 (0.02952)   0.56196 (0.02922)
  Sample size                       
                                    
  ------------- ------------------- -------------------

*Notes*: This table shows the relative root mean squared error of the predicted CTRs of the evaluation policy compared to its actual CTR. The standard errors of these root mean squared errors are in parentheses.

                            ![image](figure2.png){width="12.5cm"}

*Notes*: This figure shows estimates of the expected CTRs of the logging policy and a counterfactual policy of choosing the best action for each context. Three panels correspond to three campaigns. CTRs are multiplied by a constant for confidentiality reasons. We obtain these estimates by our proposed estimator (the DML estimator using the estimated propensity scores). Bars indicate 95% confidence intervals based on our asymptotic variance estimators developed in Proposition \[prop:step-DR\].

We apply our estimator described in Section \[section:estimation\] to empirically evaluate the design of online advertisements. This application uses proprietary data provided by CyberAgent Inc., which we described in the introduction. This company uses bandit algorithms to determine the visual design of advertisements assigned to user impressions in a mobile game. Some examples of ad designs are shown in Figure \[fig:ads\_sample\].

Our data are logged data from a 7-days A/B test on several ad “campaigns," where each campaign randomly uses either a multi-armed bandit (MAB) algorithm or a contextual bandit (CB) algorithm for each user impression. Both of the two algorithms are updated every 6 hours. As the A/B test lasts for 7 days, there are 28 time periods in total. The MAB policy and CB policy stay constant across rounds within each time period.

In the notation of our theoretical framework, reward $R_0$ is a click while action $A_0$ is one of the possible individual advertisement designs. If the algorithm is a contextual bandit, context $S_0$ is user and ad characteristics used by the algorithm. Context $S_0$ is high dimensional and has tens of thousands of possible values.

We utilize this data to examine the performance of our proposed method. For each campaign and time period, we regard the MAB algorithm as the behavior policy $\pi_b$ and the CB algorithm as the evaluation policy $\pi_e$. We use the log data from user impressions assigned to the MAB algorithm to estimate the value of the CB policy by our method. We also compute the actual value of the CB policy using the log data from user impressions assigned to the CB algorithm. We then compare the predicted value with the actual one.

We consider the relative root mean squared error (relative-RMSE) as a performance metric. Let $C$ and $T$ denote the number of campaigns ($C=4$) and the number of time periods ($T=28$), respectively. Let $N_{c,t}$ denote the number of user impressions assigned to the CB algorithm for campaign $c$ in period $t$. Let $\hat V^{\pi_e}_{c,t}$ and $\bar V^{\pi_e}_{c,t}$ be the estimated value and the actual value (the click rate) for campaign $c$ in period $t$. We define relative-RMSE as follows: $$\begin{aligned}
&relative-RMSE \\
=&~ \sqrt{\frac{1}{\sum_{c=1}^C\sum_{t=1}^T N_{c,t}}\sum_{c=1}^C\sum_{t=1}^T N_{c,t}\left(\frac{\hat V^{\pi_e}_{c,t}-\bar V^{\pi_e}_{c,t}}{\bar V^{\pi_e}_{c,t}}\right)^2}.\end{aligned}$$ As the actual click rate $\bar V^{\pi_e}_{c,t}$ varies across campaigns and time periods, we normalize the prediction error $\hat V^{\pi_e}_{c,t}-\bar V^{\pi_e}_{c,t}$ by dividing it by the actual value $\bar V^{\pi_e}_{c,t}$ to equally weight every campaign-time combination.

We calculate the standard error (statistical uncertainty) of relative-RMSE by a bootstrap-like procedure. This procedure is based on normal approximation of the distributions of $\hat V^{\pi_e}_{c,t}$ and $\bar V^{\pi_e}_{c,t}$: for each $(c,t)$, $\hat V^{\pi_e}_{c,t}\sim N(V^{\pi_e}_{c,t},\frac{\hat \sigma^{2,ope}_{c,t}}{N^{ope}_{c,t}})$ and $\bar V^{\pi_e}_{c,t}\sim N(V^{\pi_e}_{c,t},\frac{\hat \sigma^{2,online}_{c,t}}{N_{c,t}})$, where $V_{c,t}^{\pi_e}$ is the true value of policy $\pi_e$, $\hat \sigma^{2,ope}_{c,t}$ is the estimator for the asymptotic variance of $\hat V^{\pi_e}_{c,t}$ given in Proposition \[prop:step-DR\], $N^{ope}_{c,t}$ is the number of impressions used to estimate $\hat V^{\pi_e}_{c,t}$, and $\hat \sigma^{2,online}_{c,t}=\bar V^{\pi_e}_{c,t}(1-\bar V^{\pi_e}_{c,t})$ is the sample variance of the click indicator among the impressions assigned to the CB algorithm.

The standard error is computed as follows. First, we compute $\hat\sigma^{2,ope}_{c,t}$ and $\hat\sigma^{2,online}_{c,t}$ for each $(c,t)$. Second, we draw $\hat V^{\pi_e,sim}_{c,t}$ and $\bar V^{\pi_e,sim}_{c,t}$ independently from $N(\hat V^{\pi_e}_{c,t},\frac{\hat \sigma^{2,ope}_{c,t}}{N^{ope}_{c,t}})$ and $N(\hat V^{\pi_e}_{c,t},\frac{\hat \sigma^{2,online}_{c,t}}{N_{c,t}})$ for every $(c,t)$, and calculate the relative-RMSE using the draws $\{(\hat V^{\pi_e,sim}_{c,t},\bar V^{\pi_e,sim}_{c,t})\}_{c=1,...,C, t=1,...,T}$. We then repeat the second step 100,000 times, and compute the standard deviation of the simulated relative-RMSEs.

We examine the performance of three off policy evaluators, Inverse Probability Weighting (IPW; @Strehl2010, [-@Strehl2010]), Doubly Robust (DR), and Double Machine Learning (DML; our proposal) with $2$-fold cross-fitting. IPW and DR are given by $$\begin{aligned}
\hat V^{\pi_e}_{{\rm IPW}} =& \frac{1}{N}\sum_{i=1}^N \frac{\pi_e(A_{i0}|S_{i0})}{\hat\pi_{b}(A_{i0}|S_{i0})} R_{i0},\\
\hat V^{\pi_e}_{{\rm DR}} =& \frac{1}{N}\sum_{i=1}^N \{\frac{\pi_e(A_{i0}|S_{i0})}{\hat\pi_{b}(A_{i0}|S_{i0})} (R_{i0}-\hat\mu(S_{i0},A_{i0})\\
&~~~~~~~~~~~~~~~~~~+\sum_{a\in{\cal A}}\pi_e(a|S_{i0})\hat\mu(S_{i0},a)\},\end{aligned}$$ where $(\hat \pi_{b},\hat \mu)$ is the estimator of $(\pi_b,\mu)$ using the [*whole*]{} sample. $\hat V^{\pi_e}_{{\rm IPW}}$ does not satisfy the Neyman orthogonality nor use the cross-fitting. $\hat V^{\pi_e}_{{\rm DR}}$ satisfies the Neyman orthogonality but does not use the cross-fitting. As we formalize in the theory part, the DML estimator is expected to perform better than the other two estimators.

For both DR and DML, we use LightGBM [@ke2017lightgbm] as a reward estimator $\hat\mu$. For DR, we split the whole data into training and validation sets, train the reward model with the training set, and tune hyperparameters with the validation set. We then use the tuned model and the whole sample to obtain the reward estimator $\hat\mu$. For DML, we hold out a fold ($I_k$) and apply the same procedure to the rest of the data ($I_k^c$) to obtain the reward estimator $\hat\mu_k$. We repeat this for each of the $K(=2)$ folds.

For the propensity score, we use the true one or an estimated one. We compute the true propensity score by Monte Carlo simulation of the beta distribution used in Thompson Sampling the MAB algorithm uses. Our estimated propensity score is the empirical share of a selected arm in the log.

We present our key empirical result in Table \[ca-bandit-result\] and Figure \[fig:result2\]. Regardless of whether the true propensity score is available or not, DML outperforms IPW by more than 20 % and outperforms DR by more than 10 % in terms of relative-RMSE. These differences are statistically significant at 5% level. This result supports the value of our proposed estimator.

Finally, we use our proposed method to measure the performance of a new counterfactual policy of choosing the ad design predicted to maximize the CTR. We obtain the counterfactual policy by training a click prediction model with LightGBM on the data from $t = 1$ to $t = T-1$. In the training, we set the number of leaves to 20 and the learning rate to 0.01, and decide the number of boost rounds by cross-validation. We then use our estimator to predict its performance on the data where $t = T$. This final period $t=T$ contains three campaigns. The resulting off-policy evaluation results are reported in Figure \[table2\]. The results show that the counterfactual policy performs better than the existing algorithms in two of the three campaigns, with a statistically significant improvement in a campaign.

Conclusion
==========

This paper proposes a new off-policy evaluation method, by marrying the doubly robust estimator with double/debiased machine learning. Our estimator has two features. First, unlike the IPW estimator, it is robust to the bias in the estimates of the behavior policy and of the action-value function (Neyman orthogonality). Second, we use a sample-splitting procedure called cross-fitting. This removes overfitting bias that would arise without sample splitting but still makes full use of data, which makes our estimator better than DR estimators. Theoretically, we show that our estimator is $\sqrt{N}$-consistent and asymptotically normal with a consistent variance estimator, thus allowing for correct statistical inference. Our experiment shows that our estimator outperforms the standard DR and IPW estimators in terms of the root mean squared error. This result not only demonstrates the capability of our estimator to reduce prediction errors, but also suggests the more general possibility that the two features of our estimator (Neyman orthogonality and cross-fitting) may improve many variants of the DR estimator such as MAGIC [@Thomas16], SWITCH [@wang2016optimal] and MRDR [@Farajtabar2018MoreRD].

Appendices {#appendices .unnumbered}
==========

Lemmas {#lemma}
======

\[lemma:step-IS:R\] For $t=0,...,T$, $\mathbb{E}_{H\sim \pi_b}[\rho_t^{\pi_e}(H_t^{s,a})R_t]=\mathbb{E}_{H\sim \pi_e}[R_t]$.

Let $P^{\pi}(h_t^{s,a})$ denote the probability of observing trajectory $h_t^{s,a}=(s_0,a_0,...,s_t,a_t)$ when $H\sim \pi$ for some policy $\pi$. Under our data generating process, $$P^{\pi}(h_0^{s,a})=P_{S_0}(s_0)\pi(a_0|s_0),$$ and for $t\ge 1$, $$\begin{aligned}
    P^{\pi}(h_t^{s,a})=&P_{S_0}(s_0)\pi(a_0|s_0)P_S(s_1|s_0,a_0)\\
    &\cdots P_S(s_t|s_{t-1},a_{t-1})\pi(a_t|s_t).
    \end{aligned}$$ Hence, $P^{\pi_b}(h_t^{s,a})\rho_t^{\pi_e}(h_t^{s,a})=P^{\pi_e}(h_t^{s,a})$ for any $t=0,...,T$. We then have that $$\begin{aligned}
    &\mathbb{E}_{H\sim \pi_b}[\rho_t^{\pi_e}(H_t^{s,a})R_t]\\
    =~&\mathbb{E}_{H\sim \pi_b}[\rho_t^{\pi_e}(H_t^{s,a})\mathbb{E}_{H\sim \pi_b}[R_t|H_t^{s,a}]]\\
    =~&\mathbb{E}_{H\sim \pi_b}[\rho_t^{\pi_e}(H_t^{s,a})\mu(S_t,A_t)]\\
    =~&\sum_{h_{t}^{s,a}}P^{\pi_b}(h_t^{s,a})\rho_t^{\pi_e}(h_t^{s,a})\mu(s_t,a_t)\\
    =~&\sum_{h_{t}^{s,a}}P^{\pi_e}(h_t^{s,a})\mu(s_t,a_t)\\
    =~&\mathbb{E}_{H\sim \pi_e}[\mu(S_t,A_t)]\\
    =~&\mathbb{E}_{H\sim \pi_e}[\mathbb{E}_{H\sim \pi_e}[R_t|H_t^{s,a}]]\\
    =~&\mathbb{E}_{H\sim \pi_e}[R_t],
    \end{aligned}$$ where we use $\mathbb{E}_{H\sim \pi_b}[R_t|H_t^{s,a}]=\mathbb{E}_{H\sim \pi_e}[R_t|H_t^{s,a}]=\mu(S_t,A_t)$, and we use $P^{\pi_b}(h_t^{s,a})\rho_t^{\pi_e}(h_t^{s,a})=P^{\pi_e}(h_t^{s,a})$ for the fourth equality.

\[lemma:step-IS\] $\mathbb{E}_{H\sim \pi_b}[\sum_{t=0}^T\gamma^t\rho_t^{\pi_e}(H_t^{s,a})R_t]=V^{\pi_e}$.

This immediately follows from Lemma \[lemma:step-IS:R\] and the definition of $V^{\pi_e}$.

\[lemma:step-IS:S-A\] For $t=0,...,T$ and for any measurable function $g_t:({\cal S}\times{\cal A})^t\rightarrow \mathbb{R}$, $\mathbb{E}_{H\sim \pi_b}[\rho_t^{\pi_e}(H_t^{s,a})g_t(H_t^{s,a})]=\mathbb{E}_{H\sim \pi_e}[g_t(H_t^{s,a})]$.

We have that $$\begin{aligned}
    &\mathbb{E}_{H\sim \pi_b}[\rho_t^{\pi_e}(H_t^{s,a})g_t(H_t^{s,a})]\\
    =~&\sum_{h_{t}^{s,a}}P^{\pi_b}(h_t^{s,a})\rho_t^{\pi_e}(h_t^{s,a})g(h_t^{s,a})\\
    =~&\sum_{h_{t}^{s,a}}P^{\pi_e}(h_t^{s,a})g(h_t^{s,a})\\
    =~&\mathbb{E}_{H\sim \pi_e}[g_t(H_t^{s,a})],
    \end{aligned}$$ where we use $P^{\pi_b}(h_t^{s,a})\rho_t^{\pi_e}(h_t^{s,a})=P^{\pi_e}(h_t^{s,a})$ for the second equality.

\[lemma:step-IS:S\] For $t=1,...,T$ and for any measurable function $g_t:({\cal S}\times{\cal A})^{t-1}\times {\cal S}\rightarrow \mathbb{R}$, $\mathbb{E}_{H\sim \pi_b}[\rho_{t-1}^{\pi_e}(H_{t-1}^{s,a})g_t(H_{t-1}^{s,a},S_t)]=\mathbb{E}_{H\sim \pi_e}[g_t(H_{t-1}^{s,a},S_t)]$.

We have that $$\begin{aligned}
    &\mathbb{E}_{H\sim \pi_b}[\rho_{t-1}^{\pi_e}(H_{t-1}^{s,a})g_t(H_{t-1}^{s,a},S_t)]\\
    =~&\sum_{(h_{t-1}^{s,a},s_t)}P^{\pi_b}(h_{t-1}^{s,a})P_S(s_t|s_{t-1},a_{t-1})\\
    &~~~~~~~~~~~~~~~~\times\rho_{t-1}^{\pi_e}(h_{t-1}^{s,a})g(h_{t-1}^{s,a},s_t)\\
    =~&\sum_{(h_{t-1}^{s,a},s_t)}P^{\pi_e}(h_{t-1}^{s,a})P_S(s_t|s_{t-1},a_{t-1})g(h_{t-1}^{s,a},s_t)\\
    =~&\mathbb{E}_{H\sim \pi_e}[g_t(H_{t-1}^{s,a},S_t)],
    \end{aligned}$$ where we use $P^{\pi_b}(h_{t-1}^{s,a})\rho_{t-1}^{\pi_e}(h_{t-1}^{s,a})=P^{\pi_e}(h_{t-1}^{s,a})$ for the second equality.

Proof of Proposition 1 {#proof}
======================

We use Theorems 3.1 and 3.2 of @Chernozhukov2018 for the proof. We verify that $\mathbb{E}_{H\sim\pi_b}[\psi(H;\eta)]=V^{\pi_e}$ and that $\psi$ satisfies the Neyman orthogonality condition. For the first part, $$\begin{aligned}
&\mathbb{E}_{H\sim\pi_b}[\psi(H;\eta)]\\
=~&\sum_{t=0}^T \gamma^t \mathbb{E}_{H\sim\pi_b}[\rho_t^{\pi_e}(H_t^{s,a})(R_t-q_t^{\pi_e}(S_t,A_t))\\
&~~~~~~~~~~+\rho_{t-1}^{\pi_e}(H_{t-1}^{s,a})\sum_{a\in{\cal A}}\pi_e(a|S_t)q_t^{\pi_e}(S_t,a)]\\
=~&V^{\pi_e} -\sum_{t=0}^T \gamma^t \{\mathbb{E}_{H\sim\pi_e}[q_t^{\pi_e}(S_t,A_t)]\\
&~~~~~~~~~~~~~~~~~~~~-\mathbb{E}_{H\sim\pi_e}[\sum_{a\in{\cal A}}\pi_e(a|S_t)q_t^{\pi_e}(S_t,a)]\}\\
=~&V^{\pi_e}-\sum_{t=0}^T \gamma^t \{\mathbb{E}_{H\sim\pi_e}[q_t^{\pi_e}(S_t,A_t)]\\
&~~~~~~~~~~~~~~~~~~~~-\mathbb{E}_{H\sim\pi_e}[q_t^{\pi_e}(S_t,A_t)]\}\\
=~&V^{\pi_e},\end{aligned}$$ where we use Lemmas \[lemma:step-IS\], \[lemma:step-IS:S-A\] and \[lemma:step-IS:S\] for the second equality.

We now show that $\psi$ satisfies the Neyman orthogonality condition. Let $$D\rho_t^{\pi_e}(H_t^{s,a})[\tilde \pi_b]\coloneqq \left. \frac{\partial \rho_t^{\pi_e}(H_t^{s,a};\pi_b+r(\tilde \pi_b-\pi_b))}{\partial r}\right|_{r=0}$$ for any candidate $\tilde \pi_b$. Note that by Lemmas \[lemma:step-IS:S-A\] and \[lemma:step-IS:S\], $$\begin{aligned}
&\mathbb{E}_{H\sim \pi_b}[-\rho_t^{\pi_e}(H_t^{s,a})\tilde q_t^{\pi_e}(S_t,A_t)\\
&+\rho_{t-1}^{\pi_e}(H_{t-1}^{s,a})\sum_{a\in{\cal A}}\pi_e(a|S_t)\tilde q_t^{\pi_e}(S_t,a)]\\
=~&\mathbb{E}_{H\sim \pi_e}[-\tilde q_t^{\pi_e}(S_t,A_t)+\sum_{a\in{\cal A}}\pi_e(a|S_t)\tilde q_t^{\pi_e}(S_t,a)]\\
=~&\mathbb{E}_{H\sim \pi_e}[-\tilde q_t^{\pi_e}(S_t,A_t)]+\mathbb{E}_{H\sim \pi_e}[\tilde q_t^{\pi_e}(S_t,A_t)]\\
=~&0\end{aligned}$$ for $t=0,...,T$. We then have that for any candidate $\tilde \eta=(\tilde \pi_b,\{\tilde q_{t}^{\pi_e}\}_{t=0}^T)$, $$\begin{aligned}
&\left. \frac{\partial \mathbb{E}_{H\sim\pi_b}[\psi(H;\eta+r(\tilde \eta-\eta))]}{\partial r}\right|_{r=0}\\
=~&\sum_{t=0}^T \gamma^t\mathbb{E}_{H\sim \pi_b}[D\rho_t^{\pi_e}(H_t^{s,a})[\tilde \pi_b](R_t-q_t^{\pi_e}(S_t,A_t))\\
&~~~~~~~~~~~~~-\rho_t^{\pi_e}(H_t^{s,a})\tilde q_t^{\pi_e}(S_t,A_t)\\
&~~~~~~~~~~~~~+D\rho_{t-1}^{\pi_e}(H_{t-1}^{s,a})[\tilde \pi_b]\sum_{a\in{\cal A}}\pi_e(a|S_t)q_t^{\pi_e}(S_t,a)\\
&~~~~~~~~~~~~~+\rho_{t-1}^{\pi_e}(H_{t-1}^{s,a})\sum_{a\in{\cal A}}\pi_e(a|S_t)\tilde q_t^{\pi_e}(S_t,a)]\\
=~&\sum_{t=0}^T \gamma^t\mathbb{E}_{H\sim \pi_b}[D\rho_t^{\pi_e}(H_t^{s,a})[\tilde \pi_b](R_t-q_t^{\pi_e}(S_t,A_t))\\
&~~~~~~~~~~~~~+D\rho_{t-1}^{\pi_e}(H_{t-1}^{s,a})[\tilde \pi_b]\sum_{a\in{\cal A}}\pi_e(a|S_t)q_t^{\pi_e}(S_t,a)]\\
=~&\gamma^T\mathbb{E}_{H\sim \pi_b}[D\rho_T^{\pi_e}(H_T^{s,a})[\tilde \pi_b](R_T-q_T^{\pi_e}(S_T,A_T))]\\
&+\sum_{t=0}^{T-1} \gamma^t\mathbb{E}_{H\sim \pi_b}[D\rho_t^{\pi_e}(H_t^{s,a})[\tilde \pi_b](R_t-q_t^{\pi_e}(S_t,A_t)\\
&~~~~~~~~~~~~+\gamma\sum_{a\in{\cal A}}\pi_e(a|S_{t+1})q_{t+1}^{\pi_e}(S_{t+1},a))].\end{aligned}$$ Since $\mathbb{E}_{H\sim \pi_b}[R_T|H_T^{s,a}]=\mu(S_T,A_T)$ and $q_{T_e}^\pi(S_T,A_T)=\mu(S_T,A_T)$, the first term is zero by the law of iterated expectations. The second term is also zero, since for $t=0,...,T-1$, $$\begin{aligned}
&\mathbb{E}_{H\sim \pi_b}[D\rho_t^{\pi_e}(H_t^{s,a})[\tilde \pi_b]\\
&~~~~~~\times(R_t+\gamma\sum_{a\in{\cal A}}\pi_e(a|S_{t+1})q_{t+1}^{\pi_e}(S_{t+1},a))]\\
=~&\mathbb{E}_{H\sim \pi_b}[D\rho_t^{\pi_e}(H_t^{s,a})[\tilde \pi_b]\\
&\times\mathbb{E}_{H\sim \pi_b}[R_t+\gamma\sum_{a\in{\cal A}}\pi_e(a|S_{t+1})q_{t+1}^{\pi_e}(S_{t+1},a)|H_t^{s,a}]]\\
=~&\mathbb{E}_{H\sim \pi_b}[D\rho_t^{\pi_e}(H_t^{s,a})[\tilde \pi_b](\mu(S_t,A_t)\\
&+\gamma\sum_{s\in {\cal S}}P_S(s|S_t,A_t)\sum_{a\in{\cal A}}\pi_e(a|s)q_{t+1}^{\pi_e}(s,a))]\\
=&\mathbb{E}_{H\sim \pi_b}[D\rho_t^{\pi_e}(H_t^{s,a})[\tilde\pi_b]q_t^{\pi_e}(S_t,A_t)],\end{aligned}$$ where we use the recursive formulation of $q_t^\pi$ for the last equality.

The convergence results then follow from Theorems 3.1 and 3.2 of @Chernozhukov2018.

[^1]: In their experiments, @Jiang16 also implement cross-fitting as one variant of the DR estimator. Their DR estimator with cross-fitting is the same as our DML estimator if we plug the true behavior policy $\pi_b$ into DML. However, we (i) allow the behavior policy to be unknown, (ii) allow for having a large number of potentially important state variables, and (iii) present statistical properties of the estimator while they do not.
