Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
93,112 characters · 18 sections · 85 citation commands
ARMA-Design: Optimal Treatment Allocation Strategies for A/B Testing in Partially Observable Experiments
\def\spacingset#1{ {#1}} \spacingset{1.4}
\if11 \fi
\if01 {
} \fi
{\it Keywords:} ARMA Model; A/B Testing; Experimental Design; Partially Observability; Policy Evaluation; Reinforcement Learning.
\spacingset{1.7}
Background. A growing number of companies, particularly multi-sided platforms like Airbnb, DoorDash, Uber, and retail marketplaces such as Amazon and Zara are increasingly harnessing data-driven approaches to evaluate and refine their policies and products. In particular, A/B testing, which conducts online experiments to compare a standard control policy “A” to an alternate version “B”, plays a crucial role in informing business decisions within these companies and has proven invaluable for their growth and development koning2022experimentation. For instance, ride-sharing platforms, including Uber, Lyft, and DiDi Chuxing, constantly develop new order dispatching, driver repositioning, pricing policies and assess their improvements through A/B testing qin2024reinforcement. Accurate A/B testing enables decision-makers to choose better policies that meet more ride requests, enhance passenger satisfaction, increase driver income, and thus benefit the entire transportation ecosystem xu2018large.
Challenges. In many applications, the experimental units receive treatments sequentially over time. A/B testing in these experiments poses four major challenges:
Contributions. Our primary objective is to develop a statistical framework for A/B testing that addresses the above challenges. Our contributions include the development of optimal designs, efficiency indicators, statistical modeling, estimation methods, and theoretical frameworks. We detail them as follows.
To summarize, our proposal integrates cutting-edge machine learning algorithms, such as RL, with asymptotic theories derived from classical time series models in econometrics to offer guidance for policy deployment in real-world applications.
Outline. We discuss the related literature in Section (ref). In Section (ref), we introduce the controlled ARMA model, elaborate its connection with POMDPs, and develop the associated estimating procedure for ATE. We further derive the asymptotic MSEs of different designs under the small signal assumption and propose two efficiency indicators to assess their effectiveness. In Section (ref), we present the proposed algorithms to estimate the optimal design. In Section (ref), we demonstrate the efficacy of the proposed designs and efficiency indicators through two dispatch simulators and two real datasets from a ride-sharing platform. Finally, we conclude our paper in Section (ref).
Our proposal intersects with a wide range of research fields, including econometrics, statistics, management science, operational research, and machine learning. It particularly engages with three main research branches: experimental designs, POMDPs, ARMA and state space models.
Experimental Designs. The design of experiments, also known as experimental design, is a classical problem in statistics, driven by diverse applications in biology, psychology, agriculture, and engineering fisher1966design. Within the statistics literature, our proposal specifically relates to works that focus on identifying treatment allocation strategies tailored for clinical trials robbins1952some,pocock1975sequential,begg1980treatment,atkinson2007optimum,jones2009d,rosenblum2020optimal,liu2022balancing,ma2024new. These studies typically focus on non-dynamic settings (referred to as contextual bandit settings in the machine learning literature) where observations are assumed to be independent, excluding any carryover effects. In contrast, our research accommodates carryover effects and addresses the more complex challenge of temporal dynamics. While traditional crossover designs laird1992analysis,jones2003design can deal with long-lasting carryover effects, they often require extended washout periods, making them less practical for modern A/B testing with short durations.
More recently, there has been a growing body of literature in management science, econometrics, and machine learning that explores experimental designs for A/B testing in technological companies. Our work differs from them in several aspects: (i) Many papers consider settings without carryover effects over time bajari2021multiple,wan2022safe,viviano2023causal,wang2023b,basse2024randomization. (ii) Some existing works adopt an RL framework to model the experimental data glynn2020adaptive,li2023optimal, wen2024analysis, where the data follows a fully observable MDP. In contrast, our framework is more general and accommodates partial observability, which is a more typical scenario in real applications. (iii) Several recent studies focus on switchback designs where policies alternate at specified intervals under various optimality conditions hu2022switchback,bojinov2023design,xiong2023data,wen2024analysis. In contrast, our approach considers a broader class of designs that allow each treatment assignment to be influenced by the entire treatment history (see Section (ref)).
In the RL literature, the design of the experiment is also referred to as the behavior policy search problem, in which mukherjee2022revar, hanna2017data explored the optimal behavior policy by minimizing the MSE of the policy value estimator in MDPs. Meanwhile, agarwal2022reinforcement employed the D-optimal design for policy learning in MDPs. In contrast to these works, we focus on the evaluation of ATE -- the difference between two policy value estimators — and allow partial observability, offering a more realistic scenario in practice.
Finally, recent works have approached the design problem from an optimization perspective zhao2024experimental. In particular, works in the machine learning literature have proposed the use of deep learning or RL to numerically compute a Bayesian version of the optimal design foster2021deep,blau2022optimizing. In contrast to these methods, we employ a frequentist approach and focus specifically on the evaluation of ATE.
POMDPs. Partial observability often arises in real applications, including autonomous driving levinson2011towards, resource allocation bower2005resource, recommendation li2010contextual, and medical management systems hauskrecht2000planning. POMDP is the most commonly used model to characterize the partial observability of a stochastic dynamics system. Learning the optimal policy in general POMDPs requires the agent to infer the latent belief state krishnamurthy2016partially, which is both statistically and computationally intractable in general papadimitriou1987complexity,vlassis2012computational. Despite these challenges, it is possible to focus on a sub-class of POMDPs to make the estimation tractable kwon2021rl, liu2022partially. Our proposal follows this principle by introducing a controlled (V)ARMA model under a weak signal condition to streamline estimation and design. Different from existing works that proposed partial history importance weighting hu2023off or value-function-based methods uehara2022future to construct policy value estimators, we focus on the experimental design, aiming to optimize the data collection process to enhance policy evaluation.
ARMA and State Space Models. The ARMA model, a cornerstone in time series analysis, has been widely employed in various domains, particularly in econometrics brockwell1991time,hendry1995dynamic,fan2003nonlinear,box2015time,hamilton2020time. Additionally, it is closely related to state space models, which plays a vital role in analyzing continuous dynamic systems harvey1990forecasting, durbin2012time, aoki2013state, kim2017state,komunjer2020likelihood. The ARMA and state space models are also related to POMDPs, which can be seen as controlled state space models with an added dimension of the action or treatment space, allowing state transitions to be influenced by treatments krishnamurthy2016partially; see Section (ref) for detailed discussions about their connections.
In the causal inference literature, menchetti2021estimating proposed causal versions of ARIMA models. More recently, liang2023randomization proposed linear state space models for causal inference. Despite the similarity in the models, these works differ from ours primarily in their focus: they focused on estimating and inferring causal effects, whereas we concentrate on the experimental design to “optimize” the estimated causal effect. Consequently, these works did not utilize the small signal framework we propose to derive closed-form solutions for the asymptotic MSEs under different designs, which we use as criteria to optimize the design. Nor did they explore RL approaches to finding the optimal design. This difference in objectives also influences the choice of models. For instance, the causal ARIMA model is designed for settings with a single persistent treatment over time menchetti2021estimating, making it unsuitable for studying general designs.
{\singlespacing
} This section presents the proposed controlled (V)ARMA model and demonstrates its usefulness in estimating the ATE and comparing different treatment allocation strategies. We first describe the data collected from time series experiments, define the ATE for A/B testing, and introduce three commonly used designs in Section (ref). We further introduce the proposed controlled ARMA model, discuss its connections to POMDPs, and present the estimation procedure for ATE in Section (ref). Next, we propose the small signal asymptotic framework, establish the asymptotic MSE of the estimated ATE, and then derive two efficiency indicators to compare the estimation efficiency under the three designs in Section (ref). Finally, we generalize these results to accommodate multivariate observations and exogenous variables based on the proposed controlled VARMA model in Section (ref).
Data. We divide the experimental period into a series of non-overlapping time intervals, and during each of the time intervals, a specific policy or treatment is implemented. In our collaboration with a ride-sharing company, time intervals are typically set to 30 minutes or 1 hour. The data gathered from the online experiments can be summarized as a sequence of observation-treatment pairs, denoted by $\{(\mathbf{Y}_t, U_t): 1\le t\le T\}$, where $T$ represents the termination time of the experiment. Here, the notations are consistent with those used in control engineering aastrom2012introduction: $\mathbf{Y}_t$ denotes a potentially multivariate observation collected at time $t$, and $U_t$ represents a scalar treatment applied at time $t$. In detail:
ATE. Our ultimate goal lies in estimating the ATE, defined as the difference in the cumulative outcome between the treatment and the control,
provided the limit exists. Here, $\mathbb{E}_{1}$ and $\mathbb{E}_{-1}$ denote expectations under which the treatment $U_t$ is consistently set to $1$ and $-1$ at every time $t$, respectively. This objective is a central focus in A/B testing with carryover effects hu2022switchback,li2023optimal,xiong2023data,wen2024analysis. Both terms on the right-hand-side (RHS) of (ref) should be understood as potential outcomes imbens2015causal, representing the average outcome that would have been observed if either the new treatment or the control had been assigned at all times. Nonetheless, as we focus on experimental design, it eliminates concerns about unmeasured confounders. To simplify the presentation, we choose not to use potential outcome notations. Interested readers may refer to ertefaie2014constructing, luckett2019estimating and viviano2023synthetic for detailed discussions on potential outcomes in dynamic settings.
Design. In our context, each design corresponds to a sequence of treatment allocation strategies $\pi=\{\pi_t\}_{t\ge 1}$ where each $\pi_t$ specifies the conditional distribution of $U_t$ given the past data history up to time $t-1$, denoted by $\mathbf{H}_{t-1}=\{\mathbf{Y}_1,U_1,\ldots,\mathbf{Y}_{t-1}, U_{t-1}\}$. Informally speaking, each design determines the probabilities of applying the treatment and control at each time, given the past history. Our focus is on observation-agnostic designs, where each $\pi_t$ depends on $\mathbf{H}_{t-1}$ only through $\{U_1,U_2,\ldots, U_{t-1}\}$, independent of past observations. This class covers the following three special examples:
Example I: Alternating-time (AT) design. This design alternates between treatment and control at adjacent time intervals and is frequently employed in many ride-sharing companies, such as Lyft and DiDi Chuxing, to compare different order dispatching policies Chamandy2016,luo2024policy. To implement the AT design, the initial treatment $U_1$ is randomly generated with equal probabilities: $\pi_1(1)=\pi_1(-1)=0.5$. For subsequent times, we set $\pi_t(-U_{t-1}|\mathbf{H}_{t-1})=1$ and $\pi_t(U_{t-1}|\mathbf{H}_{t-1})=0$ so that $U_t=-U_{t-1}$ almost surely.
Example II: Alternating-day (AD) design. This design assigns the same treatment throughout each day and switches to the opposite treatment on the following day. Similar to the AT design, the initial treatment $U_1$ in the AD design is also uniformly randomly determined. Let $\tau$ represent the number of time intervals per day. The treatment assignment ensures that $U_1 = U_2 = \cdots = U_\tau = -U_{\tau+1} = -U_{\tau+2} = \cdots = -U_{2\tau} = U_{2\tau+1} = \cdots$, maintaining consistency within each day and alternating on a daily basis.
Example III: Uniform random (UR) design. This design independently assigns treatment and control randomly with equal probabilities each time. Specifically, $\pi_t$ remains a constant function with a value of 0.5, regardless of $t$ and $\mathbf{H}_{t-1}$. Despite its simplicity, designs of this type have been widely adopted in clinical trials.
To conclude this section, we make two remarks here. First, both AT and AD fall under the category of switchback designs, where the duration of each treatment varies from a single time interval to an entire day. Second, while many studies have explored these designs in fully observable Markovian environments, less is known about their efficacy in more realistic, partially observable environments. Addressing this gap is one of our main objectives.
{\singlespacing
}
Controlled ARMA($p, q$). We first introduce the controlled ARMA model, a sub-class of POMDPs, designed to capture carryover effects and partial observability in online experiments; see Figure (ref) for a graphical visualization. The one-dimensional controlled ARMA($p, q$) model is formulated as:
where $\mu$ denotes the intercept, $b,a_1,\ldots,a_p,\theta_1,\ldots,\theta_q$ are parameters, and by convention, $\theta_0=1$. Model ((ref)) consists of three main components:
We next illustrate how model (ref) allows carryover effects and partial observability. First, when $p>0$, the autoregressive structure and the control component allow $U_{t-1}$ to have an indirect effect on subsequent observations (e.g., $Y_t$ and $Y_{t+1}$) through its impact on $Y_{t-1}$, effectively capturing the carryover effects; see the pathway $U_{t-1}\to Y_{t-1}\to Y_t \to Y_{t+1}$ in Figure (ref). Second, when $q>0$, the inclusion of the moving average process renders the time series non-Markovian. For instance, consider the pathway $Y_{t-1}\leftarrow Z_{t-1} \leftarrow \epsilon_{t-1} \to Z_{t+1}\to Y_{t+1}$ in Figure (ref). This pathway is not blocked by $Y_t$ and $U_t$, thus violating the Markov assumption and resulting in a partially observable environment.
Finally, different sets of parameters play different roles in A/B testing: the autoregressive coefficients ($\{a_j\}_{j=1}^p$) and the control parameter ($b$) determine the ATE, whereas the moving average coefficients ($\{\theta_j\}_{j=1}^q$) influence the residual correlation, which in turn determines the optimal design. Formal statements can be found in Lemma (ref) and Theorem (ref).
Connection to POMDPs. We next show that the proposed controlled ARMA($p, q$) model is in essence a sub-class of POMDPs, which have been widely employed to model partially observable environments. Consider the following POMDP with linear state transition and observation emission functions:
In this model: (i) Both the observation $\boldsymbol{Y}_t$ and the treatment $\boldsymbol{U}_t$ can be multi-dimensional. (ii) $\boldsymbol{X}_t$ denotes a vector-valued latent state such that any dependence between the past and future will “funnel” through this latent state. (iii) $\boldsymbol{V}_t$ and $\boldsymbol{W}_t$ are the measurement errors. $F, B, H,$ and $C$ are the parameter matrices, respectively. This model can also be viewed as a variant of the linear state space or dynamic linear model, which incorporates an extra treatment variable $\boldsymbol{U}_t$.
By setting $\boldsymbol{X}_t$ to linear combinations of current and past treatments and observations, the proposed controlled ARMA($p, q$) model (ref) can be transformed into a linear POMDP. See Appendix (ref) of the Supplementary Material for formal proof. The advantage of utilizing the controlled ARMA model over a linear POMDP lies in its ability to provide concise and closed-form expressions for the asymptotic MSE of the ATE estimator (see, e.g., Corollary (ref)), which is crucial for deriving the optimal design.
According to the Wold decomposition theorem wold1938study, any stationary process can be decomposed into two mutually uncorrelated processes: a linear combination of lags of a white noise process (MA($\infty$) process) and a linear combination of its past values (AR($\infty$) process). The stationarity assumption can typically be satisfied in practice by applying periodic filtering to remove seasonal effects, as detailed in Section (ref) and our data analysis in Section (ref). This underlying principle in time series theory indicates that our model is broadly applicable and can represent a diverse range of linear POMDPs.
Estimation of ATE. We begin by deriving the closed-form expression for the ATE under the proposed controlled ARMA$(p,q)$ model.
We make three remarks. First, Assumption (ref) guarantees the ergodicity of the proposed controlled ARMA model, which in turn validates the limits in the definition of the ATE (see (ref)). Second, the ATE can be decomposed into the sum of $2b+2ab/(1-a)$ where the first term corresponds to the direct effect of $U_t$ on $Y_t$ and the second term represents the indirect effect mediated by the past observations $\{Y_{t-j}\}_{j\ge 1}$. Third, as commented earlier, the ATE is exclusively determined by the autoregressive coefficients and the control parameter, and it remains independent of the moving average coefficients. This motivates us to apply the method of moments yule1927vii,walker1931periodicity to estimate the ATE.
Notably, directly applying the ordinary least square method to minimize $\sum_t (Y_{t}-\mu -\sum_{j=1}^{p} Y_{t-j}-bU_{t})^2$ will fail to produce consistent estimators. This failure is due to the correlation between the residual $Z_{t}$ and predictors $\{Y_{t-j}\}_{j=1}^{p}$ under partial observability, as illustrated by the causal pathway $Y_{t-1}\leftarrow Z_{t-1} \leftarrow \epsilon_{t-1}\to Z_{t}$ in Figure (ref). To deal with such exogenous predictors, we employ historical observations $\{Y_{t-q-j}\}_{j=1}^{p}$ as instrumental variables angrist1996identification, which are uncorrelated with $Z_{t}$ to construct unbiased estimating equations. Specifically, by multiplying these historical observations on both sides of (ref) and taking the expectation, we obtain the following Yule-Walker equations:
It yields $p$ equations, but we have $p+2$ parameters to estimate, including $p$ autoregressive coefficients, a control parameter, and an intercept. In light of our concentration on observation-agnostic designs, under which each treatment is independent of the residual process, we further multiply $U_{t}$ and $1$ on both sides of model (ref) and take the expectation, leading to:
We next replace the expectations in (ref) and (ref) by their sample moments from $t=p+q+1$ to $T$ and construct $p+2$ estimating equations. Subsequently, we solve these equations to obtain the Yule-Walker estimators $\{\widehat{a}_j\}_j$ and $\widehat{b}$ for $\{a_j\}_j$ and $b$, respectively, by which we construct the following estimator for ATE:
By definition, the asymptotic property of (ref) depends on those of $\{\widehat{a}_j\}_{j=1}^p$ and $\widehat{b}$. However, deriving their asymptotic variances is extremely challenging, and no closed-form expressions are available to the best of our knowledge. To establish the ATE estimator's asymptotic MSE, we introduce a small signal asymptotic framework, detailed in the next section.
{\singlespacing
} We propose a small signal asymptotic framework to simplify the theoretical analysis in the ATE estimator with two key conditions:
Next, an application of the Delta method oehlert1992note to (ref) leads to
Under the first large sample condition, the third term in (ref) -- a high-order reminder term -- becomes negligible. As such, the first two terms, which measure the discrepancies between the Yule-Walker estimators and their oracle values, become the leading terms. However, as mentioned earlier, deriving their asymptotic variances remains extremely challenging under partial observability.
The second small signal condition further simplifies the calculation in two ways: (i) First, it is immediate to see that the second term is proportional to ATE. Under this condition, the second term also becomes negligible as the ATE decays to zero. The first term, therefore, becomes the sole leading term, and it suffices to calculate the asymptotic variance of the estimated control parameter. (ii) Under this condition, the influence of the treatment on the observation becomes marginal. Consequently, the sequence of treatments becomes asymptotically independent of the sequence of observations, facilitating the derivation of the asymptotic variance of $\widehat{b}$. The following theorem summarizes our findings.
The proof of Theorem (ref) is provided in Appendix (ref). Theorem (ref) may initially appear complex. To elaborate, we first narrow our analysis to the class of controlled AR models by setting $q=0$. In this simplified scenario, the residuals become uncorrelated, and the data follows a $p$-th order Markov process. Such simplification leads to the following corollary.
According to Corollary (ref), the asymptotic MSE of the ATE estimator is determined by three factors: (i) the variance of the white noise; (ii) the autoregressive coefficients; and (iii) $\xi_{\pi}$, which measures the percentage of time the new treatment is applied. Different designs affect the ATE's asymptotic variance only through $\xi_{\pi}$. In other words, designs with the same $\xi_{\pi}$ achieve the same statistical efficiency in estimating the ATE. This uniformity is due to the uncorrelated residuals in the AR model. Additionally, it turns out that any (asymptotically) balanced design with $\xi_{\pi}=0$ is optimal. This principle holds even when $q>0$, as detailed in Theorem (ref) in Section (ref). These observations align with the findings of xiong2023data, highlighting the importance of balancing periodicity in switchback designs under a different model setup.
We now turn our attention to the general controlled ARMA model with $q > 0$. We focus on the three particular designs—AT, AD, and UR—introduced in Section (ref), denoting their treatment assignment strategies as $\pi_{\text{AD}}$, $\pi_{\text{UR}}$, and $\pi_{\text{AT}}$, respectively. We derive the asymptotic MSEs of ATE estimators under these designs in the following corollary. By definition, it is evident that these three designs are balanced with $\xi_{\pi} = 0$. Specifically for the AD design, we additionally require the number of intervals per day $\tau$ to diverge to infinity as $T$ approaches infinity.
The proof is provided in Appendix (ref). According to Corollary (ref), the statistical efficiency of the three designs is primarily determined by the second term on the RHS of (ref), which depends solely on the moving average coefficients $\{\theta_j\}_{j=1}^q$. As previously noted, these coefficients directly influence the correlation of residuals, which in turn affects the designs' efficiencies. Specifically: (i) When all $\theta_j$s are non-negative, it results in non-negatively correlated residuals, and thus AT typically outperforms AD. (ii) Conversely, when the majority of residuals are non-positively correlated, AD tends to outperform AT. These observations align with the findings in xiong2023data and wen2024analysis.
Finally, Corollary (ref) motivates us to define two efficiency indicators $\textrm{EI}_{\textrm{AD}}=\sum_{k=1}^q \sum_{j=k}^q \theta_j\theta_{j-k}$ and $\textrm{EI}_{\textrm{AT}}=\sum_{k=1}^q (-1)^k \sum_{j=k}^q \theta_j\theta_{j-k}$. By (ref), it is immediate to see that
These indicators are useful for comparing the three designs. In practice, one can estimate the moving average coefficients from historical or initial experimental data and plug these estimators into the indicators to determine the most effective design among the three. In Section (ref), we discuss methodologies to search the optimal design within the broader class of observation-agnostic designs, beyond just these three.
To conclude this section, we note that although our asymptotic derivation and the subsequently proposed designs rely on the small signal condition, our proposal remains effective even with large treatment effects. This is because in these scenarios, the design problem itself might not be that critical, and any reasonable design should be able to detect these effects. Therefore, our proposal remains a safe option to use, regardless of whether this assumption holds or not.
In this section, we extend the univariate controlled ARMA model by accommodating multivariate observations and exogenous variables, derive asymptotic MSEs of the estimated ATEs, and propose the resulting efficiency indicators.
Controlled VARMA($p, q$). We define the controlled VARMA($p, q$) model with an additional exogenous variable $\bm{E}_t$ as:
where the bold vectors $\boldsymbol{\mu}, \mathbf{Y}_{t}$, $\mathbf{Z}_{t}$ and $\bm{\epsilon}_{t}$ denote the $d$-dimensional intercept, observation, residual and the white noise, respectively. The treatment $U_t$ remains binary, taking values in $\{-1, 1\}$. The purpose of introducing the extra exogenous variable $\bm{E}_t$ is to further enhance model flexibility. This variable remains unaffected by the treatments and can be regarded as the “non-stationary" components of the model, accounting for a broad range of temporal factors, such as the daily seasonal trends (see Section (ref) for the construction of $\bm{E}_t$).
Model (ref) contains four sets of parameters: (i) the autoregressive coefficient matrices $\mathbf{A}_1,\ldots,\mathbf{A}_p \in \mathbb{R}^{d \times d}$; (ii) the control coefficient vector $\mathbf{b} \in \mathbb{R}^{d}$; (iii) the moving average coefficient matrices $\mathbf{M}_1,\ldots,\mathbf{M}_q \in \mathbb{R}^{d \times d}$ and $\mathbf{M}_0 = \mathbb{I} \in \mathbb{R}^{d \times d} $ as an identity matrix; (iv) the coefficient matrix $\bm{C}$ for the extra exogenous variable.
We next introduce the no unit root assumption for the VARMA model and derive the closed-form expression for the ATE using different treatment allocation strategies.
Motivated by Lemma (ref), we similarly employ the method of moments to estimate $\{\mathbf{A}_j\}_j$ and $\mathbf{b}$ and plug in these estimators to construct the ATE estimator $\widehat{\textrm{ATE}}$. To save space, we relegate the details to Appendix (ref) in the Supplementary Material.
Asymptotic MSEs and Efficiency Indicators. Next, we analyze the asymptotic MSE of the ATE estimator in controlled VAMRA($p, q$). The following theorem extends Theorem (ref) to accommodate multivariate observations.
Similar to Corollary (ref), we next present the asymptotic MSEs of $\widehat{\textrm{ATE}}$ under AD, AT, and UR designs in the following corollary to elaborate Theorem (ref).
The proof of Theorem (ref) and Corollary (ref) are provided in Appendix (ref) in the Supplementary Material. Under the multivariate setting, we define the efficiency indicators as
According to Corollary (ref), they enable us to compare the statistical efficiency of the three designs in estimating the ATE in the controlled VARMA model.
This section focuses on the optimal observation-agnostic design, where the ATE estimator derived from the experimental data achieves the smallest asymptotic MSE. Identifying the optimal design is computationally intractable. To elaborate, each observation-agnostic design is determined by a sequence of treatment allocation strategies $\pi = \{\pi_t\}_{t=1}^T$, where each $\pi_t$ specifies the conditional distribution of $U_t$ given $U_1, \ldots, U_{t-1}$. Consider the class of deterministic treatment allocation strategies where each $\pi_t$ is a degenerate distribution. Since $U_t$s are binary, there are $2^t$ possible $\pi_t$ at each time point. Optimizing over such an exponentially growing number of strategies makes the problem NP-hard.
To address this challenge, we propose two solutions, detailed in Sections (ref) and (ref), respectively. Specifically, in Section (ref), we restrict our attention to Markov and stationary treatment allocation strategies and propose a constrained optimization algorithm to learn the resulting in-class optimal strategy. In Section (ref), we expand the search space to include general history-dependent policies and propose several optimality conditions to characterize the optimal treatment allocation strategy. These conditions significantly reduce the search space, making the computation feasible. We then develop an RL algorithm based on dynamic programming to learn the optimal treatment allocation strategy.
{\singlespacing
} To simplify the computation, we restrict attention to the class of Markov and stationary treatment allocation strategies in our first approach, where each $\pi_t$ is a function of the most recently assigned treatment $U_{t-1}$ only and remains constant with respect to $t$. In A/B testing, this policy class can be parameterized using two parameters $0 \le \alpha, \beta \le 1$, such that: \[ \mathbb{P}(U_{t+1} = 1 | U_{t} = 1) = \alpha \quad \text{and} \quad \mathbb{P}(U_{t+1} = 1 | U_{t} = -1) = \beta. \] By definition, both AT and UR are induced by policies within this class. Specifically, setting $\alpha = 0$ and $\beta = 1$ results in the AT design, whereas $\alpha = \beta = 1/2$ yields the UR design. When $\alpha = 1$, $\beta = 0$, and we alternate the initial treatment on a daily basis, it yields the AD design. This indicates the generality of the considered Markov and stationary policy class, which unifies the AD, UR and AT designs.
Additionally, the sequence $\{U_t\}_t$ forms a Markov chain with binary states. With some calculations, it can be shown that $\xi_{\pi} = (\alpha + \beta - 1) / (\beta + 1 - \alpha)$ in general. To obtain a balanced design, we set $\beta = 1 - \alpha$, leading to $\xi_{\pi} = 0$. It remains to identify the optimal $\alpha$ to minimize the asymptotic MSE of resulting the ATE estimator, which --- under the small signal asymptotic framework --- can be derived as
where $c_0 = {\sum_{j=0}^q \theta^2_j}/{(1-a)^2}$ and $c_k={\sum_{j=k}^q \theta_j \theta_{j-k} }/{(1-a)^2}$ under the controlled ARMA($p, q$) model, while under the controlled VARMA($p, q$) model we have
See Appendix (ref) for more details. This asymptotic MSE formula motivates us to compute $\alpha$ by solving the following constrained $q$-order polynomial optimization:
The above optimization can be efficiently solved using existing convex optimization techniques, such as the limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) algorithm liu1989limited. Notice that $c_k$ and the optimal number of AR and MA lags, $p$ and $q$, depend on the true model, which are typically unknown. However, as discussed in Section (ref), they can be effectively estimated or evaluated using historical data in practice. For instance, the optimal $p$ and $q$ can be selected based on the Akaike information criterion akaike1974new or the Bayesian information criterion schwarz1978estimating. Alternatively, this procedure can be applied sequentially: use current experimental data up to a specific day to learn $\{c_k\}_k$, $p$ and $q$ to estimate the optimal design. Then, this design will be applied on the following day, and the estimating procedure will continue by incorporating data from the subsequent day.
{\singlespacing
} In this section, we consider the more general history-dependent policy class and propose an RL algorithm to identify the {unrestricted} optimal treatment allocation strategy $\pi^*$. The primary objective of RL is to learn an optimal policy, a mapping from time-varying environmental features (referred to as state) to decision rules about which treatment to administer (referred to as action), in order to maximize the expected cumulative outcome (where each intermediate outcome is referred to as a reward). Most existing RL algorithms estimate the optimal policy by modeling these state-action-reward triplets over time as an MDP, wherein each reward and future state are independent of the past history given the current state-action pair.
We begin by providing an optimality condition in Theorem (ref) below to characterize $\pi^*$.
We defer the proof of Theorem (ref) to Appendix (ref) and make a few remarks: (i) Corollary (ref) in Section (ref) proves the optimality of balanced designs for AR processes. Theorem (ref) extends this to (V)ARMA processes, allowing residuals to be correlated over time. (ii) The determinism, stationarity, and $q$-dependency conditions significantly reduce the search space from over $2^{T}$ to less than $2^{q+1}$, simplifying the learning of $\pi^*$. These conditions enable us to focus on this restricted class to find $\pi^*$ by minimizing (ref). (iii) The proof of Theorem (ref) draws from existing proofs establishing the Markov and stationarity properties of the optimal policy in RL puterman2014markov, ljungqvist2018recursive. A crucial step in our proof is to construct an MDP and establish the equivalence between learning the optimal policy that maximizes the average reward in this MDP and identifying the optimal treatment allocation strategies that minimize (ref). To elaborate, we introduce the following sequence of state-action-reward triplets $(\mathbf{S}_t, A_t, R_t)_{t > q}$:
Both the future state $\mathbf{S}_{t+1}$ and the immediate reward $R_t$ are functions of $\mathbf{S}_t$ and $A_t$ only, satisfying the MDP assumption. The expected average reward in this MDP aligns with the objective function in (ref). Consequently, the optimal treatment allocation strategies satisfying (ref) are equivalent to the optimal policies under this MDP. In RL, the optimal policy is a fixed function of the current state-action pair, proving that the optimal treatment allocation strategy is deterministic, $q$-dependent, and stationary over time.
To identify $\pi^*$ that satisfies the conditions in Theorem (ref), we utilize RL as a computational tool to optimize (ref). Specifically, we construct the MDP above and apply dynamic programming to derive the optimal treatment allocation strategy. While an exhaustive policy search might be feasible when \(q\) is small, our RL approach is more computationally efficient in settings with a large \(q\).
We apply the value iteration algorithm sutton2018reinforcement for policy learning; refer to Algorithm (ref) for its pseudocode. The main idea is first to learn an optimal value function $V(s)$, which represents the maximum expected return starting from a given state $s$, and then derive the optimal policy as the greedy policy with respect to this value function (see Line 12 of Algorithm (ref)). Value iteration updates the value function iteratively using the Bellman optimality equation (see Line 8 of Algorithm (ref)) until the changes in the estimated value function are below a predefined small threshold (see Line 9 of Algorithm (ref)), indicating convergence.
We demonstrate the finite sample performance of our proposed methods using two dispatch simulators xu2018large, tang2019deep and two real datasets from a ride-sharing company. Importantly, the two simulators used in Sections (ref) and (ref) are based on physical models that simulate the behaviors of drivers and passengers. Therefore, the data generated from these environments does not necessarily follow the proposed controlled (V)ARMA model, providing a robust evaluation of our proposal under model misspecification.
Our objectives are to (i) validate the effectiveness of the proposed efficiency indicators in comparing AD, UR and AT; (ii) conduct comparisons among the following designs:
We note that Greedy is commonly used in online RL for regret minimization. TMDP and NMDP are variants of AD designs that are proven to be optimal under their model assumptions. Finally, Switch is a variant of AT design that switches back and forth over a fixed period rather than at every decision point. The optimal duration of each switch is determined by the order of the carryover effect, and we select the best duration from $\{2,5,10\}$ to report.
Environment. We simulate a synthetic ride-sharing environment as in xu2018large and li2023optimal, where drivers and customers interact within a 9$\times$9 spatial grid over 20 time steps per day:
Implementation. The outcome of interest is set to the driver's income earned at each time step. In addition to this outcome, we include two other variables in the observation: the number of unassigned orders and the number of idle drivers each time. Implementing both the proposed efficiency indicators and designs requires estimating the AR and MA parameters. To this end, we first generate a historical dataset that lasts for 50 days. Next, we apply the VARMA model to fit this dataset to estimate the AR and MA parameters. The optimal AR and MA orders, $p^*$ and $q^*$, are selected using AIC, resulting in $p^*=q^*=2$. Using these estimators, we compute the proposed efficiency indicators and proceed to implement the proposed designs, comparing them against other previously mentioned designs. Specifically, for each design, we generate 50 days of experimental data to estimate the ATE. Finally, we repeat the entire procedure 30 times to compute the MSE of the ATE estimator under each design. The oracle ATE is evaluated via the Monte Carlo method, resulting in a value of 2.24, leading to a 6% improvement.
Results. We visualize the efficiency indicators and the MSEs of ATE estimators under different designs in Figure (ref). The values of these MSEs are detailed in Table (ref). The results are summarized as follows:
Environment. We further conduct A/B testing by using a more complicated and realistic city-level order dispatching simulator tang2019deep. To mimic real-world ride-sharing markets, this simulator is trained based on a historical dataset collected from a world-leading ride-sharing company in a particular city. We do not disclose the names of the cities or the company for privacy concerns. Compared with the $9\times 9$ synthetic simulator in Section (ref), this dispatch simulator is more realistic in the following ways:
Similar to Section (ref), the observation in this city-scale simulator is also three-dimensional, including the number of orders, the number of drivers, and the driver income, which is the outcome of interest. For each design, we conduct the online experiments over four days to estimate the ATE and replicate the experiment 30 times to calculate its root MSE.
Results. The empirical distribution of efficiency indicators and the root MSE (RMSE) of ATE estimators under different designs are visualized in Figure (ref)(a) and (b), respectively. Additionally, the values of these RMSEs are reported in Table (ref) as well. We summarize the results as follows.
Data. We use two real datasets from two different cities, provided by the ride-sharing company, to create simulation environments for investigating the finite sample performance of the proposed efficiency indicators and designs. Both datasets are generated under A/A experiments, where a single order dispatching strategy is consistently deployed over time. Each dataset contains 40 days of data and is summarized as a three-dimensional time series. The first dimension records the drivers' total income at each time interval, serving as the outcome. The last two elements are the number of order requests and drivers' online time at each time interval, respectively, measuring the demand and supply of the market. The time units in the datasets differ, with the first being 30 minutes and the second being one hour. See Figure (ref) for visualizations of these three-dimensional time series.
Bootstrap-based Simulation. Figure (ref) reveals clear daily trends in both time series, with a significant rise and a subsequent decline in driver income and the number of call orders during the morning and evening peak hours. To effectively capture these seasonal patterns, we incorporate a dummy variable, $D_t$, as an exogenous variable in our controlled VARMA model to fit the three-dimensional observation. This variable is set to one during peak hours between 8 am to 8 pm and zero otherwise.
Next, we employ the parametric bootstrap to create simulated data. Specifically, we first fit the following VARMA model
and record all estimated parameters, i.e., $\widehat{\bm{\mu}}$, $\{\widehat{\mathbf{A}}_i\}_{i=1}^p$, $\widehat{\boldsymbol{\eta}}$, $\widehat{\bm{\Sigma}}$, $\{\widehat{\mathbf{M}}_{i=1}^q \}$. Next, we simulate time series $\{ \widehat{\mathbf{Y}}_{t} \}$ according to the following equation:
where $\mathbf{1}$ denotes a vector of ones, $b$ is some pre-specified parameter that determines the size of the ATE, $\{U_t\}_t$ are determined by different designs, and $\{\widehat{\mathbf{Z}}_{t}\}$ follow the estimated MA process and are generated prior to $\{\widehat{\mathbf{Y}}_{t}\}$.
Evaluation and Results. For each design and each choice of $b$, we apply the bootstrap-based simulation to generate an experimental dataset. We next apply the controlled VARMA model to this experimental dataset to estimating the ATE and evaluating its MSE. We choose an appropriate range of $b$ for each city to ensure that the resulting ATE falls between 0.5$\%$ and 2$\%$, a range that aligns with our empirical observations tang2019deep.
Given that the magnitude of the estimated ATE and the associated MSE vary with $b$, averaging all MSEs across different values of $b$ may not accurately evaluate each design. To address this, we report a performance ranking metric across the eight considered designs, which serves as a more robust measure alongside the average MSE. All results are summarized in Figure (ref) and Table (ref).
In this paper, we study experimental designs for A/B testing in partially observable online experiments where the data does not satisfy the Markov assumption. Specifically, we propose the controlled (V)ARMA model---a rich subclass of POMDPs---for fitting experimental data, establish asymptotic MSEs of ATE estimators, derive two efficiency indicators to assess the statistical efficiency of three commonly used designs, and develop two data-driven algorithms to learn the optimal observation-agnostic design. Our work bridges several vital research areas, including time series analysis, experimental design, causal inference, RL, and A/B testing, opening numerous exciting avenues for future research across these fields, covering theory, methodology, and applications. {\singlespacing
} In our ride-sharing example, we consider evaluating order dispatching policies, which are typically randomized over time and implemented in the whole city at each time, leading to temporally dependent experiments where the data generated is summarized into a single time series. However, the company is equally keen on applying different subsidizing policies in different spatial locations of the city to balance the driver supply and customer demand across the city shi2023multiagent. These experiments are inherently spatially dependent.
Spatially dependent experiments are common in many applications that involve a group of experimental units that receive (sequential) treatments across different locations ugander2013graph, baird2018optimal, johari2022experimental, leung2022rate, jia2023clustered, viviano2023causal, liu2024cluster, zhan2024estimating; see also bajari2023experimental for a recent overview. In these experiments, in addition to the carryover effect over time, the spatial spillover effect also exists where the treatment of one experimental unit can affect the outcome of others. There is growing interest in developing causal inference methods that account for spatial interference reich2021review. It would be practically interesting to integrate the proposed designs with these methodologies to adapt them to such settings.
Our proposal is model-based in that it employs classical time series models for estimating the ATE and designing the experiment. Alternatively, model-free methods that do not directly model the time series are equally applicable. Depending on how they estimate the ATE, these model-free methods can be roughly categorized into two types:
Both model-based and model-free methods have their own merits. Model-based approaches are often more data efficient, leading to less variable ATE estimators. However, they can be vulnerable to model misspecification. The first type of model-free method does not rely on a specific model, but it may fail under partial observability. The second type of method allows partial observability, but it leads to more variable ATE estimators. This increased variability is undesirable for A/B testing, particularly in settings with small sample sizes and weak signals.
Our choice of the model-based approach is guided by the principle that “all models are wrong, but some are useful” box1979all. Unlike the aforementioned model-free methods, our approach not only addresses the four practical challenges mentioned in the introduction but also demonstrates its usefulness in our numerical experiments under model misspecification. Additionally, our collaborators in the ride-sharing company prefer the model-based approach for its interpretability.
Meanwhile, the proposed controlled (V)ARMA model serves as a stepping stone. It can be extended to a variety of models for future research. For instance, autoregressive fractionally integrated moving-average (ARFIMA) models could be explored to handle long-term dependencies, which are common in practice. Similarly, more general linear state space models such as those in liang2023randomization could be considered. Despite these changes in modeling, our paper establishes a foundational framework for analysis and design. It covers both theoretical techniques, such as the small signal asymptotics, and methodological developments, including the RL algorithms, which can be adapted to these new models.
At the core of our asymptotic theories is the proposed small signal asymptotic framework, which substantially simplifies the asymptotic calculations in time dependent experiments. As mentioned earlier, it aligns with our empirical observations where most improvements from new strategies are not substantial. When this assumption is violated, our designs are not guaranteed to be optimal. However, in such cases, the treatment effect becomes non-negligible. Our approach, although potentially sub-optimal, remains consistent in detecting this effect. This ensures that our design remains safe to use, regardless of whether the signal is small or large.
Meanwhile, similar assumptions have been proposed in the literature on either A/B testing or other fields to simplify theoretical or methodological development. For instance, kuang2023weak introduced a weak signal asymptotic framework in a different context for solving multi-armed bandit problems. The main differences include: (i) Our small signal condition requires the ATE to decay to zero at an arbitrary rate, whereas kuang2023weak requires the difference in mean outcomes between different arms to decay to zero at a more restrictive parametric rate. (ii) Unlike our framework, which is designed to simplify the asymptotic analysis, their theoretical framework is developed to derive a diffusion convergence limit theorem for sequentially randomized Markov experiments.
Additionally, farias2022markovian and wen2024analysis imposed similar small signal conditions in the same context as ours for A/B testing in time dependent experiments, aiming either to derive more efficient ATE estimators methodologically or to analyze these estimators theoretically. However, their focus on Markovian environments is more restrictive than ours. They also required the difference of Markov state transition probabilities under the two treatments to be small, a condition less interpretable than our requirement for the ATE to approach zero. viviano2023causal assumed a small peer effect condition (Assumption 3) to handle the spatial spillover effect. Such an assumption shares similar spirits to ours, but is designed to develop the optimal design in spatially dependent experiments.
\spacingset{1.4}