{'fraction_non_alphanumeric': 0.06353083789442751, 'fraction_numerical': 0.019983102523601366, 'mean_word_length': 3.8191715347849176, 'pattern_counts': {'":': 0, '<': 16, '<?xml version=': 0, '>': 45, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Deep Jump Learning for Off-Policy Evaluation in Continuous Treatment Settings

We consider off-policy evaluation (OPE) in continuous treatment settings, such as personalized dose-finding. In OPE, one aims to estimate the mean outcome under a new treatment decision rule using historical data generated by a different decision rule. Most existing works on OPE focus on discrete treatment settings. To handle continuous treatments, we develop a novel estimation method for OPE using deep jump learning. The key ingredient of our method lies in adaptively discretizing the treatment space using deep discretization, by leveraging deep learning and multi-scale change point detection. This allows us to apply existing OPE methods in discrete treatments to handle continuous treatments. Our method is further justified by theoretical results, simulations, and a real application to Warfarin Dosing.


Introduction
Individualization proposes to leverage omni-channel data to meet individual needs. Individualized decision making plays a vital role in a wide variety of applications. Examples include individualized treatment regime in precision medicine (Qian and Murphy, 2011;Chakraborty and Moodie, 2013;Wang et al., 2018;Chen et al., 2020), customized pricing strategy in economics (Qiang and Bayati, 2016;Turvey, 2017), personalized recommendation system in marketing (McInerney et al., 2018), etc. Prior to adopting any decision rule in practice, it is crucial to know the impact of implementing such a policy. In medical and public-policy domains, it is risky to apply a treatment decision rule or policy online to estimate its mean outcome (see, e.g., Murphy et al., 2001;Hirano et al., 2003;Li et al., 2011). Off-policy evaluation (OPE) thus attracts a lot of attention by estimating the mean outcome under a new decision rule (or policy), i.e., the value function, using the offline data generated by a different decision rule.
First, the validity of these methods requires the mean outcome to be a smooth function over the treatment space. This assumption could be violated in applications such as dynamic pricing, where the expected demand for a product has jump discontinuities as a function of the charged price (den Boer and Keskin, 2020). Specifically, a product could attract a new segment of customers if the seller lowers the price below a certain threshold.
Thus, there will be a sudden increase in demand by a small price reduction, yielding a discontinuous demand function. Second, these kernel-based methods typically use a single bandwidth parameter. This is sub-optimal in cases where the second-order derivative of the conditional mean function has an abrupt change in the treatment space; see Section 3.1 for details. Addressing these limitations requires the development of new policy evaluation tools and theory.
Our contributions are summarized as follows. Methodologically, we develop a deep jump learning (DJL) method by integrating deep learning (LeCun et al., 2015), multi-scale change point detection (see e.g., Niu et al., 2016, for an overview), and the doubly-robust value estimators in discrete domains. Our method does not require kernel bandwidth selection. It does not suffer from the limitations of kernel-based methods. The key ingredient of our method lies in adaptively discretizing the treatment space using deep discretization. This allows us to apply the IPW or DR methods to derive the value estimator.
The discretization addresses the first limitation of kernel-based methods, allowing us to handle discontinuous value functions. The adaptivity addresses the second limitation of kernel-based methods. Specifically, it guarantees the optimality of the proposed method in cases where the second-order derivative of the conditional mean function has an abrupt change in the treatment space. Theoretically, we derive the convergence rate of the value estimator under DJL for any policy of interest, allowing the conditional mean outcome to be either a continuous or piecewise function of the treatment; see Theorems 1 and 2 for details. Under the piecewise model assumption, the rate of convergence is faster than kernel-based OPE methods. Under the continuous model assumption, kernel-based estimators converge at a slower rate when the bandwidth undersmoothes or oversmoothes the data. Proofs of our theorems rely on establishing the uniform rate of convergence of deep learning estimators; see Lemma D.1 in the supplementary article. We expect this result to be of general interest in contributing to the line of work on developing theories for deep learning methods (see e.g., Imaizumi and Fukumizu, 2019;Schmidt-Hieber et al., 2020;Farrell et al., 2021). Empirically, we show the proposed deep jump learning outperforms existing state-of-the-art OPE methods in both simulations and a real data application to warfarin dosing.

Preliminaries
We first formulate the OPE problem in continuous domains. We next review some related literature on the DR value estimator, kernel based evaluation methods, and multi-scale change point detection.

Problem Formulation
The observed offline datasets can be summarized into denotes the feature-treatment-outcome triplet for the ith individual and n denotes the total sample size. We assume these data triplets are independent copies of the population variables (X, A, Y ). Let X ∈ R p and A denote the p dimensional feature and treatment (or action) space, respectively. We focus on the setting where A is onedimensional, as in personalized dose finding and dynamic pricing. A decision rule or policy π : X → A determines the treatment to be assigned given the observed feature. We use b to denote the propensity score, also known as the behavior policy, that generates the observed data. Specifically, b(•|x) denotes the probability density function of A given X = x. Define the expected outcome function conditional on the feature-treatment pair as Q(x, a) = E(Y |X = x, A = a).
As standard in the OPE and the causal inference literature (see e.g., Chen et al., 2016;Kallus and Zhou, 2018), we assume the stable unit treatment value assumption, no unmeasured confounders assumption, and the positivity assumption are satisfied. The positivity assumption requires b to be uniformly bounded away from zero. The latter two assumptions are automatically satisfied in randomized studies. These three assumptions guarantee that a policy's value is estimable from the observed data. Specifically, for a target policy π of interest, its value can be represented by V (π) = E[Q{X, π(X)}].
Our goal is to estimate the value V (π) based on the observed data.

Doubly Robust Estimator and Kernel-Based Evaluation
For discrete treatments, Dudík et al. (2011) proposed a DR estimator of V (π) by where I(•) denotes the indicator function, Q and b(a|x) denote some estimators for the conditional mean function Q and the propensity score b, respectively. The second term inside the bracket corresponds to an augmentation term. Its expectation equals zero when Q = Q. The purpose of adding this term is to offer additional protection against potential model misspecification of Q. Such an estimator is doubly-robust in the sense that its consistency relies on either the estimation model of Q or b to be correctly specified. It can be semi-parametrically efficient whereas the importance sampling or direct method are not.
By setting Q = 0, equation 1 reduces to the IPW estimator.
In continuous treatment domains, the indicator function I{A i = π(X i )} equals zero almost surely. Consequently, naively applying equation 1 yields a plug-in estimator n i=1 Q{X i , π(X i )}/n. To address this concern, the kernel-based methods proposed to replace the indicator function in equation 1 with a kernel function Here, the bandwidth h represents a trade-off. The variance of the resulting value estimator decays with h. Yet, its bias increases with h. More specifically, it follows from Theorem 1 of Kallus and Zhou (2018) that the leading term of the bias is equal to To ensure the term in equation 14 decays to zero as h goes to 0, it requires the expected second derivative of the function Q(x, a) exists, and thus Q(x, a) needs to be a smooth function of a. However, as commented in the introduction, this assumption could be violated in certain applications.

Multi-Scale Change Point Detection
To adaptively discretize the treatment space, we leverage ideas from multi-scale change point detection literature. The change point analysis considers an ordered sequence of data, Y 1:n = {Y 1 , · · · , Y n }, with unknown change point locations, τ = {τ 1 , · · · , τ K } for some unknown integer K. Here, τ i is an integer between 1 and n − 1 inclusive, and satisfies τ i < τ j for i < j. These change points split the data into K + 1 segments. Assume there are sufficiently many data points lying within each segment such that the expected reward can be consistently estimated. Within each segment, the expected outcome is a constant function; see the left panel of Figure 1 for details. A number of methods have been proposed on estimating change points (see for example, Boysen et al., 2009;Killick et al., 2012;Frick et al., 2014;Fryzlewicz, 2014Fryzlewicz, , 2020, and the references therein), by minimizing a penalized objective function: where C is a cost function that measures the goodness-of-the-fit of the constant function within each segment and γ n K penalizes the number of change points. We remark that all the above cited works focused on either models without features or linear models. Our proposal goes beyond these works in that we consider models with features and use deep neural networks (DNN) to capture the complex relationship between the outcome and features.

Deep Jump Learning
In Section 3.1, we use a toy example to demonstrate the limitation of kernel-based methods.
We present the main idea of our algorithm in Section 3.2. Details are given in Section 3.3. For simplicity, we set the action space A = [0, 1]. Define a discretization D for the for some 0 = τ 0 < τ 1 < τ 2 < · · · < τ K−1 < τ K = 1 and some integer K ≥ 1. The union of these intervals covers A. We use J(D) to denote the set of change point locations, i.e., {τ 1 , · · · , τ K−1 }. We use |D| to denote the number of intervals in D and |I| to denote the length of any interval I.

Toy Example
As discussed in the introduction, existing kernel-based methods use a single bandwidth to construct the value estimator. Ideally, the bandwidth h in the kernel K[{A i −π(X i )}/h] shall vary with π(X i ) to improve the accuracy of the value estimator. To elaborate this, consider the function Q(x, a) = 10 max(a 2 − 0.25, 0) log(x + 2) for any x, a ∈ [0, 1]. By definition, Q(x, a) is smooth over the entire feature-treatment space. However, it has different patterns when the treatment belongs to different intervals. Specifically, for a ∈ [0, 0.5], Q(x, a) is constant as a function of a. For a ∈ (0.5, 1], Q(x, a) depends quadratically in a. See the middle panel of Figure 1 for details.
Consider the target policy π( Similarly, denote the corresponding kernel-based value estimators by is a constant function of a ∈ [0, 0.5], its second-order derivative ∂ 2 Q(x, a)/∂a 2 equals zero. In view of equation 14, when π(x) ≤ 0.5, the bias of V (1) (π; h) will be small even with a sufficiently large h. As such, a large h is preferred to reduce the variance of V (1) (π; h). When π(x) > 0.5, a small h is preferred to reduce the bias of V (2) (π; h). A simulation study is provided to demonstrate the drawback of the kernel-based methods. Specifically, we set X, A ∼ Uniform[0, 1] and generate Y |X, A ∼ N {Q(X, A), 1}. We apply the kernel-based methods with a Gaussian kernel to estimate V (1) (π) and V (2) (π) with the sample size n = 300 over 100 replications.
See Table 1 for details of the absolute error and standard deviation of V (1) (π; h) and V (2) (π; h) with two different bandwidths h = 0.4 and 1. It can be seen that due to the use of a single bandwidth, the kernel-based estimator suffers from either a large absolute error or a large variance. Middle panel: the oracle conditional mean function Q on the feature-treatment space for the toy example. Right panel: the green curve presents the oracle Q{x, π(x)} under target policy π(x) = x in the toy example; and the red curve is the fitted mean value by DJL and the pink dash line corresponds to the 95% confidence bound. Table 1: The absolute error and the standard deviation (in parentheses) of the estimated values for V (1) and V (2) , using DJL and kernel-based methods, for target policy π(x) = x in the toy example.

Methods Indicator Deep Jump Learning
To overcome this limitation, we propose to adaptively discretize the treatment space into a union of disjoint intervals such that within each interval I, the conditional mean function Q can be well-approximated by some functions q I that depend on features but not on the treatment (constant in a), i.e., Q(•, a) ≈ I∈D {I(a ∈ I)q I (•)}. By the discretization, one can apply the IPW or DR methods to evaluate the value. The advantage of adaptive discretization is illustrated in the right panel of Figure 1, where we apply the proposed DJL method to the toy example. See details of the proposed method and its implementation in Sections 3.2 and 3.3. When a ≤ 0.5, Q(x, a) is constant in a. It is likely that our procedure will not further split the interval [0, 0.5]. Consequently, the corresponding DR estimator for V (1) (π) will not suffer from large variance. When a > 0.5, our procedure will split (0.5, 1] into a series of sub-intervals, approximating Q by a step function. This guarantees the resulting DR estimator for V (2) (π) will not suffer from large bias. Consequently, the proposed value estimator achieves a smaller mean squared error than kernel-based estimators.
See Table 1 for details.

The Main Idea
We consider the following two model assumptions, which cover a variety of scenarios in practice.
Model 1: Piecewise function. Suppose for some partition D 0 of [0, 1] and a collection of functions {q I,0 } I∈D 0 .
Model 2: Continuous function. Suppose Q is a continuous function of a and x.
Model 1 covers the dynamic pricing example we mentioned in the introduction. In our simulation studies in Section 5.1, the underlying model is set to be a piecewise function in Scenarios 1 and 2. Model 2 covers the personalized dose-finding example, Scenarios 3 and 4 in our simulation studies, as well as the real data section in Section 5.2. We next detail the proposed method, which will work when either Model 1 or 2 holds.
Motivated by Model 1, our goal is to identify an optimal discretization D such that for each interval I ∈ D, Q(x, a) is approximately a constant function of a ∈ I. Specifically, under Model 1, we assume the function Q(x, a) is a piecewise function on the action space.
Within each segment I, the function Q(x, a) is a constant function of a, but can be any function of the features x. In other words, Q(x, a 1 ) = Q(x, a 2 ) for any a 1 , a 2 ∈ I. Thus, we denote the function Q(x, a) at each segment I as q I (x), which yields a piecewise function Q(x, a) = I q I (x)I(a ∈ I), as stated in equation 4. In the real applications, the true function Q(x, a) could be either a continuous function, or a piecewise function. As such, we propose to approximate the underlying unknown function Q(x, a) by these piecewise functions of a using the proposed DJL method. Such approximation allows us to derive the DR estimator based on D. The bias and variance of the resulting estimator are largely affected by the number of intervals in D. Specifically, if | D| is too small, then the piecewise approximation is not accurate, leading to a biased estimator. If | D| is too large, then D will contain many short intervals, and the resulting estimator might suffer from a large variance.
To this end, we develop a data-adaptive method to compute D. We first divide the treatment space A into m disjoint intervals: [0, 1/m), [1/m, 2/m), . . . , [(m − 1)/m, 1]. We require the integer m to diverge with the sample size n, such that the conditional mean function Q can be well-approximated by a piecewise function on these intervals. Note that these m initial intervals is not equal to D, but only serve as the initial candidate intervals.
Yet, D will be constructed by adaptively combining some of these intervals. We find in our numerical studies that the size of the final partition | D| is usually much less than m (see Table 6 in Appendix B for more details). In practice, we recommend to set the initial number of intervals m to be proportional to the sample size n, i.e., m = n/c for some constant c > 0. The performance of the resulting value estimator is not overly sensitive to the choice of c.
We define B(m) as the set of discretizations D such that each interval I ∈ D corresponds to a union of some of the m initial intervals. Each discretization D ∈ B(m) is associated with a set of functions {q I } I∈D . We model these q I using DNNs, to capture the complex dependence between the outcome and features. When Q(•, a) is well-approximated by for some regularization parameter γ n and some function class of DNNs Q I . Here, the penalty term γ n |D| in equation 5 controls the total number of intervals in D, as in multi-scale change point detection. A large γ n results in few intervals in D and a potential large bias of the value estimator, whereas a small γ n procedures a large number of intervals in D, leading to a noisy value estimator. The theoretical order of γ n is detailed in Section 4. In practice, we use cross-validation to select γ n . We refer to this step as deep discretization. Details of solving equation 5 are given in Section 3.3.
Given D and { q I : I ∈ D}, we apply the DR estimator in equation 1 to derive the value estimate for any target policy of interest π, i.e., where b I (x) is some estimator of the generalized propensity score function Pr(A ∈ I|X = x).
We call this method as the deep jump learning. We remark that the proposed method yields a consistent value estimator allowing the function Q to be either a continuous or piecewise function of the treatment. Under different model assumptions, we derive the corresponding rate of convergence of our method in Section 4.

The Complete Algorithm for Deep Jump Learning
We present the details for DJL in this section. To further reduce the bias of the value estimator in equation 6, we employ a data splitting and cross-fitting strategy (Chernozhukov et al., 2017). That is, we use different subsets of data samples to estimate the discretization D and to construct the value estimator. Our algorithm consists of three steps: data splitting, deep discretization, and cross-fitting. We detail each of these steps below.
Step 1: data splitting. We divide all n samples into L disjoint subsets of equal size, where L denotes the indices of samples in the th subset for = 1, · · · , L. Let L c = {1, 2, · · · , n} − L as the complement of L . Data splitting allows us to use one part of the data, i.e., L c , to train machine learning models for the conditional mean function and propensity score function, and the remaining part, i.e., L , to estimate the value. We aggregate the resulting estimates over different subsets to get full efficiency, as summarized in the third step.
Step 2: deep discretization. For each = 1, · · · , L, we propose to apply deep discretization to compute a discretization D ( ) and { q ( ) I : I ∈ D ( ) } by solving a version of equation 5 using the data subset in L c only. We next present the computational details for solving this optimization. Our algorithm employs the pruned exact linear time method (Killick et al., 2012) to identify the change points with a cost function that involves DNN training. Specifically, for any interval I, define q ( ) I as the minimizer of arg min where |L c | denotes the number of samples in L c . Define the cost function C ( ) (I) as the minimum value of the objective function equation 7, i.e, ; number of initial intervals m; penalty term γ n ; target policy π. Local: Bellman function Bell ∈ R m ; partitions D; DNN functions { q I , b I : I ∈ D}; a vector τ ∈ N m ; a set of candidate point lists R. Output: the value estimator for target policy V (π). I. Split all n samples into L subsets as {L 1 , · · · , L L }; V (π) ← 0; II. Initialization: 1. Set even segment on the action space with m pieces: Create a function to calculate cost C with inputs (l, r): If C(l, r) == N U LL: Return C(l, r); III. For = 1, · · · , L: 1. Set the training dataset as L c = {1, 2, · · · , n} − L ; 2. Bell(0) ← −γ n ; D = [0, 1]; τ ← N ull; R(0) ← {0}; 3. Apply the pruned exact linear time method to get partitions: For v * = 1, . . . , m: (iv) r ← l; l ← τ (r); 6. Evaluation using testing dataset L : return V (π)/n .

Algorithm 1: Deep Jump Learning
Computation of D ( ) relies on dynamic programming (Friedrich et al., 2008). For any integer Our algorithm recursively updates the Bellman cost function for v * = 1, 2, · · · by where R v * is the candidate change points list. For a given v, the right-hand-side of equation 8   To solve equation 8, we search the optimal change point location v that minimizes Bell(v * ). This requires deep learning to estimate q ( ) This procedure is iterated to compute Bell(v * ) and τ (v * ) for v * = 1, . . . , m, to find the best change-point set for interval [0, 1]. The optimal partition D ( ) is determined by the values stored in τ . Specifically, we initialize D ( ) = [τ (m)/m, 1], r = m and recursively update D ( ) by setting D ( ) ← D ( ) ∪ [τ (r)/m, r/m) and r ← τ (r), as in dynamic programming (Friedrich et al., 2008).
Step 3: cross-fitting. For each interval in the estimated optimal partition D ( ) , let b ( ) I (x) denote some estimator for the propensity score Pr(A ∈ I|X = x). In a randomized study, the density function b(a|x) is known to us and we set b ( ) To deal with data from observational studies, we estimate the generalized propensity score with deep learning using the training dataset L c as b proposed value estimator for V (π) is to aggregate over V for = 1, · · · , L via cross-fitting, given by, Note the samples used to construct V inside bracket are independent from those to estimate

Theory
We investigate the theoretical properties of the proposed DJL method. All the proofs are provided in the supplementary article. Without loss of generality, assume the support X = [0, 1] p . To simplify the analysis, we focus on the case where the behavior policy b is known to us, which automatically holds for data from randomized studies. We focus on the setting where the conditional mean function Q is a smooth function of the features; see A1 below. Specifically, define the class of β-smooth functions, also known as Hölder smooth functions with exponent β, as for some constant c > 0, where β denotes the largest integer that is smaller than β and ∆ α denotes the differential operator ∆ α denote the differential operator: When β is an integer, β-smoothness essentially requires a function to have bounded derivatives up to the βth order. The Hölder smoothness assumption is commonly imposed in the current literature (see e.g., Farrell et al., 2021), which is a special example of the function classes that can be learned by neural nets. Meanwhile, the proposed DJL method is valid when Q(x, a) is a nonsmooth function of x as well (see e.g., Imaizumi and Fukumizu, 2019). Our theory thus can be further generalized to any function class that can be learned by neural nets at a certain rate. We introduce the following conditions.

Properties under Model 1
We first consider Model 1 where the function Q(x, a) takes the form of equation 4. As commented, this assumption holds in applications such as dynamic pricing. Without loss of generality, assume q I 1 ,0 = q I 2 ,0 for any two adjacent intervals I 1 , I 2 ∈ D 0 . This guarantees In addition, for any policy π such that for any We make a few remarks. First, the result in (i) imply that deep discretization correctly identifies the number of change points. The result in (ii) imply that any change point in D 0 can be consistently identified. In particular, J( D ( ) ) corresponds to a subset of {1/m, 2/m, · · · , (m − 1)/m}. For any true change point τ in D 0 , there will be a change point in D ( ) that approaches τ at a rate of n −2β/(2β+p) up to some logarithmic factors. Second, it can be seen from the proof of Theorem 1 that the two error terms O{n −2β/(2β+p) log 8 n} and O(n −1/2 ) in (iii) correspond to the bias and standard deviation of the proposed value estimator, respectively. When 2β > p, the bias term is negligible. A Wald-type confidence interval can be constructed to infer V (π). The assumption 2β > p allow the deep learning estimator to converge at a rate faster than n −1/4 . Such a condition is commonly imposed in the literature for inferring the average treatment effect (see e.g., Chernozhukov et al., 2017;Farrell et al., 2021). When β < p, i.e., the underlying conditional mean function Q is not smooth enough, the proposed estimator suffers from a large bias and might converge at a rate that is slower than the usual parametric rate. This concern can be addressed by employing the A-learning method (see e.g., Murphy, 2003;Schulte et al., 2014;Shi et al., 2018). The A-learning method is more robust and requires weaker conditions to achieve the parametric rate. Specifically, it only requires the difference Q(x, 1) − Q(x, 0) to belongs to Φ(β, c). This is weaker than requiring both Q(x, 1) and Q(x, 0) to belongs to Φ(β, c). Third, to ensure the consistency of the proposed value estimator, we require that the distribution of the random variable π(X) induced by the target policy does not have point-masses at the change point locations. This condition is mild. For nondynamic policies where π(X) = π 0 almost surely, it requires π 0 / ∈ J(D 0 ). We remark that the set J(D 0 ) has a zero Lebesgue measure on [0, 1]. For dynamic policies, it automatically holds when π(X) has a bounded density on [0, 1].

Properties under Model 2
We next consider Model 2 where the function Q(x, a) is continuous in the treatment space.
for all a 1 , a 2 ∈ [0, 1], x ∈ X , and some constant L > 0. Assume (A1) and (A2), and m is proportional to n and γ n is proportional to max{n −3/5 , n −2β/(2β+p) log 9 n}. Then for any target policy π, When 4β > 3p, the convergence rate is given by O p (n −1/5 ). We remark that the above upper bound is valid for any target policy π. convergence poses some unique challenges in deriving the results of Theorems 1 and 2. We need to control the initial number of the intervals m to be proportional to n and the order of penalty term γ n , so that uniform convergence rate can be established across all intervals.
To address this difficulty, we derive the tail inequality to bound the rate of convergence of the DNN estimator and use the Bonferroni's correction to establish the uniform rate of convergence.

Comparison with Kernel-Based Methods
To simplify the analysis, we assume the kernel function is symmetric, the nuisance function estimators Q and b are set to their oracle values Q and b, and that 4β > 3p. Suppose Model 1 holds. In Appendix C, we show that the convergence rate of kernel-based methods is given by O p (n −1/3 ) with optimal bandwidth selection. In contrast, the proposed estimator converges at a faster rate of O p (n −1/2 ). Suppose Model 2 holds. In Appendix C, we show that the convergence rate of kernel-based methods is given by kernel-based estimators converge at a slower rate when the bandwidth undersmoothes or oversmoothes the data. In addition, as we have commented in Section 3.1, in cases where the second-order derivative of Q has an abrupt change in the treatment space, kernel-based methods suffer from either a large bias, or a large variance. Specifically, when h is either much larger than n −1/5 or much smaller than n −3/5 , our estimator converges at a faster rate. Kernel-based estimators could converge at a faster rate when Q has a uniform degree of smoothness over the entire treatment space and the optimal bandwidth parameter is correctly identified.

Simulation Studies
In this section, we investigate the finite sample performance of the proposed method on the simulated and real datasets, in comparison to three kernel-based methods. The computing infrastructure used is a virtual machine in the AWS Platform with 72 processor cores and 144GB memory.

Simulation Settings
Simulated data are generated from the following model: . Consider the following different scenarios: S3 (toy example): Q(x, a) = 10 max{a 2 − 0.25, 0} log(x (1) + 2); The function Q(x, a) is a piecewise function of a under Scenarios 1 and 2, and is continuous under Scenarios 3 (toy example considered in Section 3.1) and 4. We set the target policy to be the optimal policy that achieves the highest possible mean outcome. The dimension of the features is fixed to p = 20. We consider four choices of the sample size, corresponding to n = 50, 100, 200 or 300.
We compare the proposed DJL method with three kernel-based methods (Kallus and In addition, we estimate the propensity score function using MLP as well. We set m = n/10 to achieve a good balance between the absolute error and the computational cost (see Figure 3 in Appendix B for details). The averaged computational time are summarized in Table 3 with additional results under large sample sizes n = 1000 ∼ 10000 in Table 4, in Appendix B. Overall, it takes a few minutes (less than 1 min for n = 50 and 14 mins for n = 300) to implement DJL, whereas the runtime of Kallus and Zhou (2018)'s method is 365 mins for sample size n = 50 and over 48 hours for n = 300. Thus, as suggested in Kallus and Zhou (2018), to implement their method, we first compute h * using data with sample size n 0 = 50.
To accommodate data with different n, we adjust h * by setting h * {n 0 /n} 0.2 . To implement Colangelo and Lee (2020)'s estimator, we consider a list of bandwidths suggested in their paper, given by h = cσ A n −0.2 with c ∈ {0.5, 0.75, 1.0, 1.5} and σ A is the standard deviation of the treatment. We then manually select the best bandwidth such that the resulting value estimator achieves the smallest mean squared error. The kernel-based method (SLOPE) by Su, Srinath and Krishnamurthy (2020) adopted the Lepski's method for bandwidth selection. In their implementation, they used the IPW estimator to evaluate the value. For a fair comparison, we replace it with DR to make the resulting estimator more efficient.
The average estimated value and its standard deviation over 100 replicates are illustrated in Figure 2 for different methods, with detailed values reported in Table 5 in Appendix B.
In addition, we provide the size of the final estimated partition under DJL in Table 6 in Appendix B, which is much smaller than m in most cases. It can be seen from Figure 2 that the proposed DJL method is very efficient and outperforms all competing methods in almost all cases. We note that the proposed method performs reasonably well even when the sample size is small (n = 50). In contrast, kernel-based methods fail to accurately estimate the value even in some cases when n = 300. Among the three kernel-based OPE approaches, we observe that the method developed by Su, Srinath and Krishnamurthy (2020) performs better in general. A potential limitation of our method is that it takes a longer computational time than the method of Colangelo and Lee (2020). To speed up the dynamic programming algorithm, for instance, the total variation or group-fused-lasso-type penalty can be used as a surrogate of the L 0 penalty to reduce the computational complexity (see e.g., Harchaoui and Lévy-Leduc, 2010). We use the min-max normalization to convert the range of the dose level A into [0, 1].
To compare among different methods, we calibrate the dataset to generate simulated outcomes. This allows us to use simulated data to calculate the bias and variance of each value estimator. Specifically, we first estimate the function Q(x, a) via the MLP regressor using the whole dataset. The goodness-of-the-fit of the fitted model under the MLP regressor is reported in Table 7 in Appendix B. We next use the fitted function Q(X, A) to simulate the data. For a given sample size N , we first randomly sample N feature- Next, for each j, we generate the outcome y j according to

This yields a simulated dataset
{(x j , a j , y j ) : 1 ≤ j ≤ N }. We are interested in evaluating the mean outcome under the optimal policy as π (X) ≡ arg max a∈[0,1] Q(X, a).
We apply the proposed DJL method and the three kernel-based methods to the calibrated Warfarin dataset. Absolute errors, standard deviations, and mean squared errors of the estimated values under the optimal policy are reported in Table 2 over 20 replicates with different evaluation methods. It can be observed from Table 2  For the last part, the construction of the final value estimator based on L-fold cross fitting, which repeats the above two steps L times. Therefore, by putting the above results together, the total expected computational complexity of the proposed DJL is O{L * m * 2 * ne(d − 1)w 2 }. Note that the computation for the last part (i.e., cross-fitting) can be easily implemented in parallel computing, and thus the total expected computational complexity of the proposed DJL can be reduced to O{m * 2 * ne(d − 1)w 2 }.

B Additional Experimental Results
We include additional experimental results in this section. First, the number of initial intervals m represents a trade-off between the estimation bias and the computational cost, as illustrated in Figure 3. In practice, we recommend to set m = n/10. When n is small, the performance of the resulting value estimator is not overly sensitive to the choice of c as long as c is not too large. See the left panel of Figure 3 for details. When n is large, we further investigate the computational capacity of the proposed method by setting m = n/10 for large sample sizes and report the corresponding computational time in Table 4. We use Scenario 1 and consider the sample size chosen from n ∈ {1000, 2000, 5000, 10000} for illustration. It turns out that such a choice of c can still handle datasets with a few thousand observations. Here, we use parallel computing to process each fold, as our algorithm employs data splitting and cross-fitting. This largely facilitates the computation, leading to shorter computation time compared to those listed in Table 3. Finally, when n is extremely large, setting m = n/10 might be computationally intensive. In addition to parallel computing, there are some other techniques we can use to handle datasets with large sample size. For instance, in the change-point literature, Lu et al. (2017) proposed an intelligence sampling method to identify multiple change points with long time series data. Their method would not lose much statistical efficiency, but is much more computationally efficient. It is possible to adopt such an intelligence sampling method to our setting for adaptive discretization.
This would enable our method to handle large datasets.
Define a policy π such that the density function of π(X) equals We aim to show for such Q and π, the best possible convergence rate of kernel-based estimator is n −1/3 .
We first consider its variance. Suppose the conditional variance of Y |A, X is uniformly bounded away from 0. Similar to Theorem 1 of Colangelo and Lee (2020) We next consider its bias. Since the behavior policy is known, the bias is equal to Using the change of variable a = ht + π(X), the bias equals Consider any 0 < h ≤ for some sufficiently small > 0. The bias is then equal to Under the symmetric condition on the kernel function, the above quantity is equal to 2 3 1/2 Consequently, the bias is lower bounded by O (1) Model 1 holds, we can show the bias equals Similarly, we can show the right-hand-side is lower bounded by O(1)h. This implies that the convergence rate is at least O(1)(n −1 h −1 + h) under Model 2.

D Technical Proof
Throughout the proof, we use c, C, c 0 ,c, c * , etc., to denote some universal constants whose values are allowed to change from place to place. Let O i = {X i , Y i } denote the data summarized from the ith observation. For any two positive sequences {a n } n and {b n } n .
The notation a n b n means that there exists some universal constant c > 1 such that c −1 b n ≤ a n ≤ cb n for any n. The notation a n ∝ b n means that there exists some universal constant c > 0 such that a n ≤ cb n for all n.
Proofs of Theorems 1 and 2 rely on Lemmas D.1, D.2 and D.3. In particular, Lemma D.1 establishes the uniform convergence rate of q ( ) I for any I whose length is no shorter than o(γ n ) and belongs to the set of intervals: for some integers i 1 and i 2 that satisfy 0 ≤ i 1 < i 2 < m} ∪ {[i 3 /m, 1] : for some integers i 3 that satisfy 0 ≤ i 3 < m}.
To state this lemma, we first introduce some notations. For any such interval I, define the function q I,0 (x) = E(Y |A ∈ I, X = x). It is immediate to see that the definition of q I,0 here is consistent with the one defined in equation 4 for any I ⊆ D 0 .
Here, the expectation in equation 11 is taken with respect to a testing sample X.
Lemma D.2 Assume either conditions in Theorem 1 or 2 are satisfied. Then there exists some constantC > 0 such that the followings hold with probability at least 1 − O(n −2 ): For any 1 ≤ ≤ L, I ∈ I(m) and |I| ≥ cγ n , We first present the proofs for these three lemmas. Next we present the proofs for Theorems 1 and 2.

D.1 Proof of Lemma D.1
The number of folds L is bounded. It suffices to derive the uniform convergence rate for each . By definition, q ( ) I is the minimizer of the least square loss, arg min q∈Q I i∈L c I(A i ∈ I)|Y i − q(X i )| 2 . It follows that i∈L c for any q and I.
The first term on the right-hand-side measures the approximation bias of the class of deep neural networks. Since E[I(A ∈ I){Y − q I,0 (X)}|X] = 0, the second term corresponds to the stochastic error. The rest of the proof is divided into three parts. In Part 1, we bound the approximation error. In Part 2, we bound the stochastic error. Finally, we combine these two parts together to derive the uniform convergence rate for q ( ) I . Part 1. Under the given condition, we have Q(•, a) ∈ Φ(β, c), b(a|•) ∈ Φ(β, c) for some c > 0 and any a. We now argue that there exists some constant C > 0 such that q I,0 ∈ Φ(β, C) for any I. This can be proven based on the relation that Specifically, we have that sup x |q I,0 (x)| ≤ sup a,x |Q(x, a)| ≤ c. Suppose β ≤ 1. For any x 1 , x 2 ∈ X , consider the difference |q I,0 (x 1 ) − q I,0 (x 2 )|. Under the positivity assumption, we have inf a,x b(a|x) ≥ c * for some c * > 0. It follows that |q I,0 (x 1 ) − q I,0 (x 2 )| ≤ I |Q(x 1 , a) − Q(x 2 , a)|b(a|x 1 )da I b(a|x 1 )da Consequently, q I,0 ∈ Φ(β, c + 2c 2 /c 2 * ). Suppose β > 1. Then both Q(•, a) and b(a|•) are β -differentiable. By changing the order of integration and differentiation, we can show that q I,0 (x) is β -differentiable as well. As an illustration, when β < 2, we have β = 1. According to the chain rule, we have Moreover, using similar arguments in proving q I,0 ∈ Φ(β, c + 2c 2 /c 2 * ) when β < 1, we can show that all the partial derivatives of q I,0 (x) up to the β th order are uniformly bounded for all I. In addition, all the β th order partial derivatives are Hölder continuous with exponent β − β . This implies that q I,0 ∈ Φ(β, C) for some constant C > 0 and any I.
It is shown in Lemma 7 of Farrell et al. (2021) that for any > 0, there exists a deep neural network architecture that approximates q I,0 with the uniform approximation error upper bounded by , and satisfies W I ≤C −p/β (log −1 + 1) and L I ≤C(log −1 + 1) for some constantC > 0. These upper bounds will be used later in Part 2. The detailed value of will be specified below. It follows that for any I, the bias term can be upper bounded by i∈L c We next provide an upper bound for the right-hand-side. Since A has a bounded probability density function, the variance Var{I(A i ∈ I)} is upper bounded by EI(A i ∈ I) ≤c |I| for some universal constantc > 0. It follows from Bernstein's inequality that for any t and I. Set t I = 6 max(c n|I| log n, |I| log n), the right-hand-side is upper bounded by n −4 . Since m n and the number of intervals I in I(m) is upper bounded by m 2 , it follows from Bonferroni's inequality that As such, with probability at least 1 − O(n −2 ), we have that i∈L c I(A i ∈ I) − |L c |EI(A ∈ I) ≤ t I uniformly for all I, or equivalently, i∈L c I(A i ∈ I) ≤ |L c |c|I| + t I . Consider a subset of intervals I with |I| ≥ cγ n for any constant c > 0. Under the given conditions on γ n , we have i∈L c I(A i ∈ I) ≤ nc * |I|, for any I such that |I| ≥ cγ n , for some constantc * > 0. It follows from equation 12 that the following holds with probability at least 1 − O(n −2 ): for any I ∈ I(m) such that |I| ≥ cγ n , we have i∈L c I(A i ∈ I)|q I,0 (X i ) − q(X i )| 2 ≤c * 2 n|I|.
Part 2. For the function class of deep neural networks Q I , we use θ I to denote the parameters in deep neural networks. This allows us to represent Q I as {q I (•, θ I ) : θ I } We will apply the empirical process theory (see e.g., Van Der Vaart and Wellner, 1996) to bound the stochastic error. Let θ I be the estimated parameter in q ( ) I . Define σ 2 (I, θ) = E I(A ∈ I)|q I,0 (X) − q I,0 (X, θ)| 2 , for any θ and I. Consider two separate cases, corresponding to σ(I, θ I ) ≤ |I| 1/2 (n|I|) −β/(2β+p) and σ(I, θ I ) > |I| 1/2 (n|I|) −β/(2β+p) , respectively. We focus our attentions on the latter class of intervals. In Part 3, we will show that for those intervals, for some universal constant O(1). This implies that for any I, we have We consider bounding a scaled version of the stochastic error, Its absolute value can be upper bounded by where the supremum is taken over all θ such that σ(I, θ) > |I| 1/2 (n|I|) −β/(2β+p) .
For a given θ, the empirical sum has zero mean. Under the boundedness assumption on Y , its variance is upper bounded by some universal constant. In addition, each quantity σ −1 (I, θ)I(A i ∈ I){Y i − q I,0 (X i )}{q I,0 (X i , θ) − q I,0 (X i )} is upper bounded by O(1)|I| −1/2 (n|I|) β/(2β+p) for some universal constant O(1). This allows us to apply the tail inequality developed by Massart et al. (2000) to bounded the empirical process. See also Theorem 2 of Adamczak et al. (2008). Specifically, for all t > 0 and I that satisfies σ(I, θ I ) > |I| 1/2 (n|I|) −β/(2β+p) , we obtain with probability at least 1 − exp(t) that for some constantc > 0. By setting t = 3 log n, the probability 1 − exp(t) = 1 − n −3 . Notice that the number of intervals I is upper bounded by O(n 2 ), under the condition that m is proportional to n. By Bonferroni's inequality, we obtain that equation 16 holds with probability at least 1 − O(n −2 ) for any I. Under the given condition on γ n , for any interval I such that |I| ≥ cγ n , the last term on the right-hand-side of equation 16 is o( √ n). It follows that the following occurs with probability 1 − O(n −2 ), for all I such that |I| ≥ cγ n and σ(I, θ I ) > |I| 1/2 (n|I|) −β/(2β+p) .
We next provide an upper bound for EZ(I). Toward that end, we will apply the maximal inequality developed in Corollary 5.1 of Chernozhukov et al. (2014). We first observe that the class of empirical sum indexed by θ belongs to the VC subgraph class with VC-index upper bounded by O(W I L I log(W I )). It follows that for any I such that |I| ≥ cγ n , σ(I, θ I ) > |I| 1/2 (n|I|) −β/(2β+p) , EZ(I) ∝ nW I L I log(W I ) log n + W I L I log(W I ) log n.
Based on the upper bounds on W I and L I developed in Part 1, the right-hand-side is upper bounded by where O(1) denotes some universal constant. It is of the order O{n 1/2 (n|I|) p/(4β+2p) log 4 n}.
This together with equation 16 and equation 17 yields that with probability at least 1 − O(n −2 ), the scaled stochastic error is upper bounded by n 1/2 (n|I|) p/(4β+2p) log 4 n. As such, with probability at least 1 − O(n −2 ), we obtain that for any I such that |I| ≥ cγ n , σ(I, θ I ) > |I| 1/2 (n|I|) −β/(2β+p) . By Cauchy-Schwarz inequality, the left-hand-side can be further upper bounded by where O(1) denotes some universal positive constant. This completes the proof for Part 2.
Part 3. Combining the results in Part 1 and Part 2, we obtain that for any I such that with probability at least 1 − O(n −2 ). As for the left-hand-side, we notice that i∈L c Using similar arguments in Part 2, we can show that the second line is upper bounded by nσ 2 (I, θ I )/8 + O(1)(n|I|) p/(2β+p) log 8 n, with probability at least 1 − O(n −2 ), for any I such that |I| ≥ cγ n , σ(I, θ I ) > |I| 1/2 (n|I|) −β/(2β+p) . Since L c ≥ n/2, we obtain This yields the desired uniform upper bound for σ 2 (I, θ I ). We thus obtain equation 15 holds with probability at least 1 − O(n −2 ).
Under the assumption that the density function b(a|x) is uniformly bounded away from zero, we obtain for some constant c > 0. This assertion thus follows.

D.2 Proof of Lemma D.2
The assertion can be proven in a similar manner as Part 2 of the proof of Lemma D.1. We omit the details to save space.

D.3 Proof of Lemma D.3
Consider a given interval I ∈ D ( ) . Suppose |I| < cγ n . The value of the constant c will be determined later. Then, for sufficiently large n, we can find some interval I ∈ I(m) ∩ D ( ) that is adjacent to I. Thus, we have I ∪ I ∈ I(m), and hence Notice that the left-hand-side of the above expression is nonnegative. It follows that By definition, we have q ( ) It follows that By equation 18, this further implies that Under (A2), the function q I is uniformly upper bounded from above. It thus follows from Cauchy-Schwarz inequality that for some constant c 0 > 0. Using similar arguments in showing equation 13, we can show that with probability at least 1 − O(n −2 ), the following evens hold for all I ∈ I(m), for some constant c 1 > 0. The right-hand-side shall be larger than or equal to γ n . Consequently, we have either |I| ≥ c 2 γ n or |I| ≥ c 2 nγ 2 n / log n for some constant c 2 > 0. Under the given condition on γ n , we obtain that |I| ≥ c 2 γ n for sufficiently large n. The proof is hence completed.

D.4 Proof of Theorem 1
Since the number of folds L is a fixed integer. We will show the assertions in (i) and (ii) holds for each , with probability at least 1 − O(n −2 ). The proof is divided into three parts. In Part 1, we show the consistency of the estimated change point locations and that | D ( ) | ≥ |D 0 | with probability at least 1 − O(n −2 ). In Part 2, we prove that | D ( ) | = |D 0 | with probability at least 1 − O(n −2 ) and derive the rate of convergence of the estimated change point locations and the estimated function Q. In Part 3, we derive the rate of convergence for the value estimator.
Part 1. We first show the consistency of the estimated change-point locations. Assume |D 0 | > 1. Otherwise, the assertion | D ( ) | ≥ |D 0 | trivially hold. Consider the partition D = {[0, 1]} which consists of a single interval and a zero function Q. By definition, we have Under the boundedness assumption on Y , we obtain that |L c |γ n | D ( ) | ≤ C 0 (|L c | + γ n ) for some constant C 0 > 0 and hence for sufficiently large n, as γ n → 0.

Notice that
The second line is non-negative. Under Lemmas D.2 and D.3, the third line is lower bounded by −C 1 I∈ D ( ) (n|I|) p/(p+2β) log 8 n for some constant C 1 > 0 with probability at least 1 − O(n −2 ). In view of equation 19, it can be further lower bounded by −2C 0 C 1 γ −1 n n p/(p+2β) log 8 n. By equation 19 and the given condition on γ n , the third line is o(n). It follows that with probability at least 1 − O(n −2 ).

It follows that
with probability at least 1 − O(n −2 ). It follows from equation 20 that with probability at least 1 − O(n −2 ).
Let us consider η * 2 . We observe that By the uniform approximation property of deep neural networks, there exists some q * I ∈ Q I such that i∈L c |q I,0 (X i ) − q * I (X i )| 2 ∝ n(n|I|) −2β/(2β+p) .
See Part 1 of the proof of Lemma D.1 for details. Similar to equation 13, we can show that the following events occur with probability at least 1 − O(n −2 ), for some constant c 0 > 0 and any I ∈ D 0 . It follows that for some constantc > 0. This together with equation 23 yields that with probability at least 1 − O(n −2 ).
Let K = |D 0 |. For any integer k such that 1 ≤ k ≤ K − 1, let τ * 0,k be the change point location that satisfies τ * 0,k = i/m for some integer i and that |τ 0,k − τ * 0,k | < m −1 . Denoted by D * the oracle partition formed by the change point locations {τ * 0,k } K−1 k=1 . Set τ * 0,0 = 0, τ * 0,K = 1 and q * * Since Y is a bounded variable, q * I is uniformly bounded for any I. The right-hand-side is upper bounded by K k=1 i∈L c I(A i ∈ ∆ k ). Similar to equation 13, The later is upper bounded by O(log n), with probability at least 1 − O(n −2 ).It follows that with probability at least 1 − O(n −2 ). By definition, Combining this together with equation 25 yields that It follows from equation 24 and the condition γ n → 0 that with probability at least 1 − O(n −2 ). Under the event defined above, we show that max τ ∈J(D 0 ) minτ ∈J( D ( ) ) |τ − τ | ≤ δ for any constant δ > 0. This yields the consistency of our estimated change point locations.
Specifically, under the condition that q I 1 ,0 = q I 2 ,0 for any adjacent I 1 , I 2 ∈ D 0 , we have E|q I 1 ,0 (X) − q I 2 ,0 (X)| 2 > 0. Let δ 0 denote the minimum distance between two change point locations. Since the change points are fixed, δ 0 is a fixed positive value. For a given Since δ 0 > δ, the conditional mean function Q is a piecewise function of A in the intervals [a, τ 0 ] and [τ 0 , b]. The left-hand-side thus equals The function q I,0 that minimizes the above objective is given by Consequently, the left-hand-side of equation 28 is greater than or equal to with | D ( ) | ≥ |D 0 |, we obtain that | D ( ) | = |D 0 |. As such, we have that I∈ D ( ) EI(A ∈ I)|Q(X, A) − q I,0 (X)| 2 ∝ n −2β/(p+2β) log 8 n Using similar arguments in establishing the consistency of the estimated change point locations, we can show that under the above event, we have that max τ ∈J(D 0 ) minτ ∈J( D ( ) ) |τ − τ | ∝ n −2β/(p+2β) log 8 n. This completes the proof of this part.
Part 3. For any target policy π, we define a random policy π D ( ) according to the partition D ( ) as follows: where b(I|x) denotes the propensity score function Pr(A ∈ I|X = x). Note that 1 0 π D ( ) (a|x)da = I⊆ D ( ) I{π(x) ∈ I} = 1 for any x. Consequently, π D ( ) is a valid random policy. Since the behavior policy is known, the proposed doubly-robust estimator corresponds to an unbiased estimator for L −1 L =1 V (π D ( ) ). Using similar arguments in the causal inference literature on deriving the asymptotic property of doubly-robust estimators (Chernozhukov et al., 2017), we can show that V (π D ( ) ) = O p (n −1/2 ).
Based on the results obtained in the first two parts, it follows from Cauchy-Schwarz inequality that Note that V (π D ( ) ) = E Similarly, we can show V (π) = I 0 ∈D 0 Eq I 0 (X)I{π(X) ∈ I 0 }.

It follows that
As q I 0 is uniformly bounded, the left-hand-side is upper bounded by Based on the results obtained in Part 2, for each I 0 ∈ D 0 , there exists some I

D.5 Proof of Theorem 2
We break the proof into two parts. In Part 1, we introduce an auxiliary lemma and present its proof. In Part 2, we derive the convergence rate of the proposed value estimator.
Part 1. We first introduce the following lemma.
Lemma D.4 For any interval I ∈ I(m) with |I| γ n and any interval I ∈ D ( ) with I ⊆ I , we have with probability approaching 1 that E|q I,0 (X) − q I ,0 (X)| 2 ≤C|I| −1 γ n , for some constantC > 0.
We next prove Lemma D.4. For a given interval I ∈ D ( ) , the set of intervals I considered in Lemma D.4 can be classified into the following three categories.
Notice that for any interval I 0 , Using similar arguments in bounding the stochastic error term in Part 2 of the proof of Lemma D.1, we can show with probability approaching 1 that the right-hand-side is of the order O{n −2β/(2β+p) log 8 n}, for any I 0 ∈ I(m). As such, we obtain with probability approaching 1 that By the positivity assumption, we have with probability approaching 1 that E[{q I (X) − q I (X)} 2 |{O i } i∈L c ] = O(γ n |I| −1 ) + O{(n|I|) −2β/(2β+p) log 8 n}, uniformly for any I and I . The proof is hence completed by noting that γ n is at least of the order O(n −2β /(2β + p)) log 8 n.
In Case 1, it follows from the Lipschitz property of the conditional mean function Q that |Q(x, a 1 ) − Q(x, a 2 )| ≤ Lγ 1/3 n , for any a 1 , a 2 ∈ I and x. By definition, the function q I can be represented as q I (x) = I Q(x, a)ω(a, x)da for some weight function ω such that I ω(a, x)da = 1. It follows that the right-hand-side of equation 33 is upper bounded by Lγ 1/3 n . In Case 2, for any a ∈ I , we can find an interval I ⊆ I , a ∈ I with length proportional to γ 1/3 n . Using similar arguments in Case 1, we can show that |Q(x, a) − q I,0 (x)| ≤ Lγ 1/3 n . By Lemma D.4 and the Cauchy-Schwarz inequality, we have E|q I,0 (X) − q I ,0 (X)| ≤ C γ 2/3 n =C 1/2 γ 1/3 n , with probability approaching 1. It follows that the right-hand-side of equation 33 is upper bounded by (L + √C )γ 1/3 n , with probability approaching 1. As such, the bias of the proposed estimator is upper bounded by (L + C )γ 1/3 n , with probability approaching 1.
We next consider the standard deviation of our estimator. The proposed estimator is can be represented by L −1 L =1 V (π) where V (π) is the value estimator constructed based on the samples in {O i } i∈L . Since the propensity score function is known to us, each V (π) is unbiased to V (π D ( ) ). Under the positivity assumption and the boundedness assumption on Y and q I , the variance of V (π) is upper bounded by |L | −1 inf I∈ D ( ) |I| −1 . By Lemma D.3, it is upper bounded by O(n −1 γ −1 n ). As such, the standard deviation of our estimator is upper bounded by O(n −1 γ −1 n ). As such, the convergence rate is given by O p (γ 1/3 n + n −1/2 γ −1/2 n ). By setting γ n = n −3/5 , the rate is given by O p (n −1/5 ). The proof is hence completed.