Appendix

#  Off-policy Evaluation in Doubly Inhomogeneous Environments

Zeyu Bian1, Chengchun Shi2, Zhengling Qi3 and Lan Wang1

1Department of Management Science, University of Miami
2Department of Statistics, London School of Economics and Political Science
3 Department of Decision Sciences, The George Washington University

###### Abstract

This work aims to study off-policy evaluation (OPE) under scenarios where two
key reinforcement learning (RL) assumptions – temporal stationarity and
individual homogeneity are both violated. To handle the “double
inhomogeneities”, we propose a class of latent factor models for the reward
and observation transition functions, under which we develop a general OPE
framework that consists of both model-based and model-free approaches. To our
knowledge, this is the first paper that develops statistically sound OPE
methods in offline RL with double inhomogeneities. It contributes to a deeper
understanding of OPE in environments, where standard RL assumptions are not
met, and provides several practical approaches in these settings. We establish
the theoretical properties of the proposed value estimators and empirically
show that our approach outperforms competing methods that ignore either
temporal nonstationarity or individual heterogeneity. Finally, we illustrate
our method on a data set from the Medical Information Mart for Intensive Care.

Keywords: Double Inhomogeneities; Off-policy Evaluation; Reinforcement
Learning; Two-way Fixed Effects Model.

##  1 Introduction

Reinforcement learning (RL, Sutton and Barto,, 2018) aims to optimize an
agent’s long-term reward by learning an optimal policy that determines the
best action to take under every circumstance. RL is closely related to the
dynamic treatment regimens (DTR) or adaptive treatment strategies in
statistics research for precision medicine (Murphy,, 2003; Robins,, 2004; Qian
and Murphy,, 2011; Zhang et al.,, 2012; Chakraborty and Moodie,, 2013; Kosorok
and Moodie,, 2015; Zhu et al.,, 2017; Tsiatis et al.,, 2019; Kosorok and
Laber,, 2019; Qi et al.,, 2020; Zhou et al., 2022a, ), which seeks to obtain
the optimal treatment policy in finite horizon settings with a few treatment
stages that maximizes patient’s expected outcome. Nevertheless, statistical
methods for DTR mentioned above normally cannot handle large or infinite
horizon settings. They require the number of trajectories to tend to infinity
to achieve a consistent estimation, unlike RL, which can work even with finite
sample size under certain conditions. In addition to precision medicine, RL
has been applied to various fields, such as games (Silver et al.,, 2016),
ridesharing (Xu et al.,, 2018), mobile health (Liao et al.,, 2021) and
robotics (Levine et al.,, 2020).

In this article, we focus on the topic of off-policy evaluation (OPE) in RL,
whose objective is to evaluate the value function of a given target policy
using data collected from a potentially different policy, known as the
behavior policy. OPE is important in applications in which directly
implementing the target policy involves potential risks and high costs. For
instance, in healthcare, it would be expensive (in terms of both time and
budget) to conduct a randomized experiment to recruit a large number of
individuals and follow them up for the duration of the entire experiment.
Meanwhile, it might be unethical to directly apply a new treatment policy to
some individuals without offline validation. It is therefore important to
study and develop RL methods only using the observed and historical data
(offline RL), and OPE appears to be particularly vital in offline RL.
Generally speaking, the existing OPE methods can be divided into four
categories: model-based methods (Yin and Wang,, 2020), importance sampling-
based methods (Liu et al.,, 2018; Wang et al.,, 2021), value function-based
methods (Luckett et al.,, 2020; Hao et al.,, 2021; Liao et al.,, 2021; Shi et
al.,, 2022), and doubly robust estimation methods (Jiang and Li,, 2016; Thomas
and Brunskill,, 2016; Bibaut et al.,, 2021; Kallus and Uehara,, 2022; Liao et
al.,, 2022), which typically are a fusion of the former two approaches. See
Uehara et al., (2022) and the references therein for an overview.

Motivation. Most methods in the RL literature would require the following two
critical assumptions in order to obtain a valid estimator: temporal
stationarity and individual homogeneity. The temporal stationarity assumption
requires that the system dynamics for each subject do not depend on the time
whereas individual homogeneity requires the system dynamics at each time to be
identical across all subjects. Nonetheless, both conditions are likely to be
violated in many of RL applications, e.g., mobile health and infectious
disease control; see Hu et al., (2022) for a thorough discussion and examples.
This work draws partial motivation from the longitudinal data of septic
patients obtained from the Medical Information Mart for Intensive Care (MIMIC-
III, Johnson et al.,, 2016), a database that contains information on critical
care patients. Sepsis is a severe and potentially fatal condition that occurs
when the human body’s response to an infection injures its own tissues and
organs (Singer et al.,, 2016). It can progress rapidly and cause multiple
organ failures, resulting in a decline in a patient’s health and an increased
risk of death. Prompt treatment of sepsis is thus essential for improving
patient outcomes and reducing mortality rates. However, the heterogeneity in
patients’ response to sepsis treatment (Evans et al.,, 2021), as well as a
potential non-stationary environment (the data includes patients’ medical
information over 10 years) make it extremely challenging to effectively manage
the illness using existing statistical methods. Our analysis provides insights
into the impact of different treatment options on patient outcomes which helps
develop more effective and personalized approaches to sepsis care.

In the statistical literature, Li et al., (2022) developed a hypothesis
testing procedure to assess the stationarity assumption in RL, based on which
a policy learning procedure is proposed to handle possibly nonstationary
environments. Chen et al., 2022a ; Chen et al., 2022b  developed a transferred
Q-learning algorithm and an auto-clustered policy iteration algorithm to
handle heterogeneous data. However, these methods require either temporal
stationarity or individual homogeneity, and would fail in doubly inhomogeneous
environments when both assumptions are violated. Hu et al., (2022) proposed an
algorithm to adaptively split the data into rectangles in which the system
dynamics are the same over time and across individuals. They studied policy
learning instead of OPE. In addition, they imposed a latent group structure
over time and population. This structural assumption can be violated when the
dynamics vary smoothly over both population and time.

Challenges. OPE is substantially more challenging under the doubly
inhomogeneous environments. First, the evaluation target is different. In
particular, most existing solutions developed in doubly homogeneous
environments have predominantly focused on evaluating the expected long-term
reward following the target policy aggregated over time and population. In
contrast, the following four time- and/or individual-specific values are of
particular interest in the presence of double inhomogeneities:

  1. 1.

The expected long-term reward aggregated over both time and population;

  2. 2.

The expected long-term reward aggregated over time for a given subject;

  3. 3.

The expected reward at a given time aggregated over population;

  4. 4.

The expected reward at a given time for a given subject.

Second, an unresolved challenge is how to efficiently borrow information over
time and population for OPE. On one hand, to account for the subject
heterogeneity or temporal nonstationarity, one could conduct OPE based on the
data within each individual trajectory or at a given time. However, this
approach may result in estimators with high variance. On the other hand,
naively pulling data over population and time without careful considerations
would lead to biased estimators.

Contributions. This work makes the following contributions. First, to the best
of our knowledge, it is the first study to investigate OPE in doubly
inhomogeneous RL domains. Unlike prior works that primarily focused on
evaluating the average effect over time and population, we provide a
systematic approach by examining values that are specific to time and/or
individuals. These values hold particular importance in the context of double
inhomogeneities.

Second, we present a comprehensive framework for doubly inhomogeneous OPE
which comprises both model-free and model-based methods. To effectively
utilize information in the presence of temporal nonstationarity and individual
heterogeneity, we introduce a class of two-way doubly inhomogeneous decision
process (TWDIDP) models and develop policy value estimators based on these
models. Our proposal shares similar spirits with the two-way fixed effects
model that is widely studied in the economics and social science literature
(Angrist and Pischke,, 2009; Wooldridge,, 2010; De Chaisemartin and
d’Haultfoeuille,, 2020; Arkhangelsky et al.,, 2021; Imai and Kim,, 2021).
Nonetheless, our model is more complicated in that the current treatment not
only affects its immediate outcome, but impacts the future outcomes as well
through its effect on the future observation via the transition function. In
contrast, the fixed effects models commonly employed in the panel data
literature tend to exclude the possibility of carryover effects (Imai and
Kim,, 2019; Arkhangelsky et al.,, 2021).

Finally, we systematically explore the theoretical properties of the proposed
model-free method. In particular, we derive the convergence rates of various
proposed value estimators, showing that the estimated average effect,
individual-specific effect, time-specific effect and individual- and time-
specific effect converge at a rate of $(NT)^{-1/2}$, $T^{-1/2}$, $N^{-1/2}$
and $\min^{-1/2}(N,T)$ respectively, up to some logarithmic factors, where $N$
is the number of trajectories and $T$ is the number of time points. We further
establish the limiting distributions of these estimators.

Organization. The rest of this paper is organized as follows. In Section 2, we
introduce the proposed doubly inhomogeneous decision process model to
incorporate temporal nonstationarity and individual heterogeneity. In Sections
3 and 4, we present our proposed model-free and model-based methods. We
analyze their statistical properties in Section 5. A series of comprehensive
simulation studies are conducted in Section 6. Finally, in Section 7, we
illustrate the proposed approach using the MIMIC-III dataset.

##  2 Two-way Doubly Inhomogeneous Decision Processes

Data. We first describe the dataset. We assume the offline data consist of $N$
independent trajectories, each with $T$ many time points, and can be
summarized as the following observation-action-reward triplets
$\\{(O_{i,t},A_{i,t},R_{i,t}):1\leq i\leq N,1\leq t\leq T\\}$ where $i$
indexes the $i$th individual and $t$ indexes the $t$th time point. For
example, in mobile health applications, $O_{i,t}\in\mathbb{R}^{d}$ denotes the
vector of covariates measured from the $i$th individual at time $t$, $A_{i,t}$
denotes the treatment assigned to the $i$th individual at time $t$, and
$R_{i,t}\in\mathbb{R}$ denotes the $i$th individual’s clinical outcome at time
$t$. Let $\mathcal{O}$ and $\mathcal{A}$ denote the observation and action
space, respectively. We assume $\mathcal{A}$ is a discrete space whereas,
$\mathcal{O}$ is a compact subspace of $\mathbb{R}^{d}$ where $d$ denotes the
dimension of the observation and that the reward is uniformly bounded. The
bounded rewards assumption is commonly imposed in the RL literature (see e.g.,
Fan et al.,, 2020; Li et al.,, 2023).

Model. We next present the proposed two-way doubly inhomogeneous decision
process model. In the RL literature, a common practice is to employ the Markov
decision process (MDP, Puterman,, 2014) to model the data generating process.
Assuming both the observation and reward spaces are discrete, the MDP model
essentially requires the reward and future observation to be conditionally
independent of the past data history given the current observation-action pair
so that the system dynamics are uniquely determined by the following Markov
transition function $p$,

$\displaystyle{\mathbb{P}}(O_{i,t+1}=o^{\prime},R_{i,t}=r|A_{i,t}=a,O_{i,t}=o,\\{O_{i,j},A_{i,j},R_{i,j}\\}_{1\leq
j<t})=p(o^{\prime},r|a,o),$ (2.1)

which is assumed to be doubly homogeneous, i.e., constant over time and
population.

The proposed model relies on two key assumptions. First, to model double
inhomogeneities, we assume the existence of a set of individual- and time-
specific latent factors $\\{U_{i}\\}_{i=1}^{N}$ and $\\{V_{t}\\}_{t=1}^{T}$
conditional on which the Markov assumption holds. More specifically, for any
$i$ and $t$, we assume

$\displaystyle{\mathbb{P}}(O_{i,t+1}=o^{\prime},R_{i,t}=r|U_{i}=u_{i},V_{t}=v_{t},A_{i,t}=a,O_{i,t}=o,\\{O_{i,j},A_{i,j},R_{i,j},V_{j}\\}_{1\leq
j<t})$ (2.2) $\displaystyle=$ $\displaystyle p(o^{\prime},r|u_{i},v_{t},a,o).$

######  Remark 1.

Unlike (2.1), the transition function in (2.2) is both individual- and time-
dependent due to the inclusion of $U_{i}$ and $V_{t}$. The individual-specific
factors can be viewed as certain individual baseline information (e.g.,
educational background, salary) that does not vary over time whereas the time-
specific factors correspond to certain external factors (e.g., holidays) that
have common effects on all individuals.

######  Remark 2.

Both $\\{U_{i}\\}_{i=1}^{N}$ and $\\{V_{t}\\}_{t=1}^{T}$ are unobserved in
practice, leading to the violation of the Markov assumption. Indeed, the
proposed data generating process can be viewed as a special class of partially
observable MDPs (POMDPs, Sutton and Barto,, 2018) where the unobserved factors
either do not evolve over time (e.g., $\\{U_{i}\\}_{i=1}^{N}$) or do not vary
across individuals (e.g., $\\{V_{t}\\}_{t=1}^{T}$). More generally, one may
allow the latent factors to evolve over both time and population. However,
this makes the subsequent policy evaluation extremely challenging. In
contrast, our proposal decomposes these factors into individual-only and time-
only effects, which can be consistently estimated when both $N$ and $T$
diverge to infinity. We also remark that latent factor models are widely used
in finance (Ross,, 1976), economics (Bai and Ng,, 2002), psychology (Bollen,,
2002) and machine learning (Nickel et al.,, 2011).

######  Remark 3.

When the observation and reward spaces are continuous, we use
$p(\bullet,\bullet|U_{i},V_{t},A_{i,t},O_{i,t})$ to denote the conditional
density function of $(O_{i,t+1},R_{i,t})$ given $U_{i}$, $V_{t}$, $A_{i,t}$
and $O_{i,t}$. Our proposal is applicable to both discrete and continuous
observation/reward spaces.

Second, we impose an additivity assumption, which requires the transition
function $p$ to be additive in $u$, $v$ and $(a,o)$, i.e.,

$\displaystyle
p(o^{\prime},r|u_{i},v_{t},a,o)=\pi_{u}p_{u_{i}}(o^{\prime},r|u_{i})+\pi_{v}p_{v_{t}}(o^{\prime},r|v_{t})+\pi_{0}p_{0}(o^{\prime},r|a,o)$
(2.3)

for some non-negative constants $\pi_{u}$, $\pi_{v}$, and $\pi_{0}$ that
satisfy $\pi_{u}+\pi_{v}+\pi_{0}=1$ and some unknown conditional probability
density (mass) functions $p_{u}$, $p_{v}$ and $p_{0}$.

The additivity assumption in (2.3) essentially assumes that the transition
function corresponds to a mixture of $p_{u}$, $p_{v}$ and $p_{0}$, with the
mixing weights given by $\pi_{u}$, $\pi_{v}$ and $\pi_{0}$, respectively.
Under the additivity assumption, $p_{u}$ and $p_{v}$ correspond to the
individual- and time-specific effects, and are independent of the current
observation-action pair. The function $p_{0}$ corresponds to the main effect
shared over time and subjects.

Multiplying both sides of (2.3) by $r$ and integrating with respect to $r$ and
$o^{\prime}$, we obtain that

$\displaystyle R_{i,t}=\theta_{i}+\lambda_{t}+r_{1}(a,o)+\varepsilon_{i,t},$
(2.4)

where $\theta_{i}=\pi_{u}\int rp_{u_{i}}(o^{\prime},r|u_{i})drdo^{\prime}$,
$\lambda_{t}=\pi_{v}\int rp_{v_{t}}(o^{\prime},r|v_{t})drdo^{\prime}$,
$r_{1}(a,o)=\pi_{0}\int rp_{0}(o^{\prime},r|a,o)drdo^{\prime}$, and
$\varepsilon_{i,t}=R_{i,t}-{\mathbb{E}}(R_{i,t}|A_{i,t}=a,O_{i,t}=o)$ has
conditional mean zero. To ease notation, in what follows, we assume the latent
factors $\\{U_{i}\\}_{i}$ and $\\{V_{t}\\}_{t}$ are fixed, and use
$\\{u_{i}\\}_{i}$ and $\\{v_{t}\\}_{t}$ to denote their realizations. All the
expectations mentioned below are assumed to be implicitly conditional on
$\\{U_{i}\\}_{i}$ and $\\{V_{t}\\}_{t}$. Notably, the policy value of interest
is also defined as the conditional expectation of the reward given
$\\{U_{i}\\}_{i}$ and $\\{V_{t}\\}_{t}$. Models of this type (Equation 2.4)
are referred to as the two-way fixed-effects model in the panel data
literature (see e.g., Angrist and Pischke,, 2009; Imai and Kim,, 2021).
Nonetheless, our model has the distinctive feature that it allows the current
treatment to not only affect the immediate outcome, but also to impact the
future outcomes through its effect on the future observation via the
transition function in (2.3).

######  Remark 4.

Our additivity assumption (2.3) is motivated by the increased popularity of
the fixed-effect models in the panel data literature, due to its ability to
account for unobserved variables. As commented by Green et al., (2001), fixed
effects regression can scarcely be faulted for being the bearer of bad
tidings. On the other hand, when $\pi_{0}=1$, the proposed model reduces to
the standard MDP studied in doubly homogeneous environments.

Estimand. Finally, we define our target estimand of interest. A policy
prescribes how an agent should act and make decisions. Mathematically, it maps
the space of observed data history to a probability mass function on the
action space, determining the probability that a given individual receives a
given treatment at each time point. Throughout this paper, we focus on
evaluating stationary policies where the action selection probability depends
on history only through the current observation and this dependence is
stationary over time. More specifically, following a given stationary policy
$\pi$, the $i$th individual will receive treatment $a$ with probability
$\pi(a|O_{i,t})$. This policy class includes observation-agnostic policies,
i.e., $\pi(a|O_{i,t})=\pi(a)$ for any $a$. Meanwhile, the proposed method can
be extended to evaluate possibly nonstationary policies.

For a given target policy $\pi$, we define the following four estimands of
interest: (i) the average effect
$\eta^{\pi}\equiv(NT)^{-1}\sum_{i=1}^{N}\sum_{t=1}^{T}{\mathbb{E}}^{\pi}(R_{i,t})$;
(ii) the individual-specific effect given the observed initial observation
$\eta^{\pi}_{i}\equiv T^{-1}\sum_{t=1}^{T}{\mathbb{E}}^{\pi}(R_{i,t}|O_{i,1})$
(which is a scalar instead of a function of the initial observation); (iii)
the time-specific effect $\eta^{\pi}_{t}\equiv
N^{-1}\sum_{i=1}^{N}{\mathbb{E}}^{\pi}(R_{i,t})$ and (iv) the individual- and
time-specific effect
$\eta^{\pi}_{i,t}\equiv{\mathbb{E}}^{\pi}(R_{i,t}|O_{i,1})$. Here, the
notation $\mathbb{E}^{\pi}$ means that the expectation is taken by assuming
the system dynamics follow the target policy $\pi$. In defining
$\eta_{i}^{\pi}$ and $\eta_{i,t}^{\pi}$, we include $O_{i,1}$ in the
conditioning set to eliminate their variability resulting from marginalizing
over the initial observation distribution. This is reasonable, as the initial
observation distribution may no longer be identical across different subjects
due to individual heterogeneity, making it impossible to infer consistently
from the data.

We focus on estimating (iv) $\eta^{\pi}_{i,t}$ in the next two sections, based
on which estimators for (i)–(iii) can be easily derived by taking the average
over time and/or population.

######  Remark 5.

In stationary environments with certain mixing conditions (see e.g., Bradley,,
2005), the system will reach its stationary distribution. In that case, it is
immediate to see that $\lim_{T}\eta_{i}^{\pi}=\lim_{t}\eta_{i,t}^{\pi}$ and
$\lim_{T}\eta^{\pi}=\lim_{t}\eta_{t}^{\pi}$. Meanwhile, under individual
homogeneity, we obtain $\lim_{T}\eta^{\pi}=\lim_{T}\eta_{i}^{\pi}$. As such,
the aforementioned four targets are asymptotically the same in the classical
doubly homogeneous environments.

##  3 Model-free OPE

In this section, we develop model-free methodologies to learn
$\eta^{\pi}_{i,t}$: the $i$th subject’s average reward at time $t$ given
$O_{i,1}$. Model-free methods construct the policy value estimator without
directly learning the transition function. Compared to model-based methods
which directly learn the transition function to derive the estimator, they are
preferred in settings with large observation space, or where the transition
function is highly complicated and can be easily misspecified.

Challenge. Before presenting our proposal, we outline the challenges in
consistently estimating the policy value. First, existing model-free methods
developed in the RL literature (see e.g., Luckett et al.,, 2020; Liao et al.,,
2021; Kallus and Uehara,, 2022; Liao et al.,, 2022; Shi et al.,, 2022) focused
on learning the long-term reward in a stationary environment. These methods
are not applicable to learn the expected reward at a given time with
nonstationary transition functions. Second, in the DTR literature, backward
induction or dynamic programming (Bellman,, 1957; Bather,, 2000) is widely
employed to evaluate the value function in the sparse reward setting where the
reward is obtained at the last stage and all the immediate rewards equal zero.
It is applicable to evaluate ${\mathbb{E}}^{\pi}(R_{i,t})$ in nonstationary
environments. Nonetheless, it requires all individual trajectories to follow
the same distribution and is thus inapplicable to our setting.

Q-function. Our proposal extends the classical backward induction method to
the doubly inhomogeneous environments. To begin with, we define the following
individual- and time-specific Q-function

$\displaystyle
Q^{\pi}_{i,t_{1},t_{2}}(o,a)={\mathbb{E}}^{\pi}(R_{i,t_{2}}|A_{i,t_{1}}=a,O_{i,t_{1}}=o),$
(3.1)

for any $1\leq i\leq N$ and $1\leq t_{1}\leq t_{2}\leq T$. To elaborate on
this definition, we consider two particular choices of $t_{1}$. First, when
$t_{1}=t_{2}$, (3.1) reduces to the conditional mean of $R_{i,t_{2}}$ given
($A_{i,t_{2}}$, $O_{i,t_{2}}$), which equals
$\theta_{i}+\lambda_{t}+r_{1}(A_{i,t},O_{i,t})$ (see Equation (2.4)) under
additivity. Second, when $t_{1}=1$, it is immediate to see that

$\displaystyle\eta^{\pi}_{i,t_{2}}=\sum_{a}Q^{\pi}_{i,1,t_{2}}(O_{i,1},a)\pi(a|O_{i,1}).$
(3.2)

As such, it suffices to learn $Q^{\pi}_{i,1,t}$ to construct estimators for
$\eta_{i,t}^{\pi}$.

######  Remark 6.

In both the DTR and RL literature, the Q-function is typically defined as the
cumulative reward starting from a given time $t_{1}$. Our Q-function in (3.1)
differs in that: (i) it is individual-specific where the subscript $i$ encodes
its dependence upon the latent factor $u_{i}$; (ii) it is the conditional mean
of the immediate reward at time $t_{2}$ only instead of the cumulative reward
since our objective here lies in evaluating ${\mathbb{E}}^{\pi}(R_{i,t_{2}})$.

Backward induction. We propose to use backward induction to compute an
estimated Q-function $\widehat{Q}_{i,1,t}^{\pi}$ for $Q_{i,1,t}^{\pi}$ and
then plug the estimator into (3.2) to construct the policy value estimator. To
begin with, consider the reward function $\\{Q^{\pi}_{i,t,t}\\}_{i,t}$. As
shown in (2.4), under the two-way fixed-effect model, we have
$Q^{\pi}_{i,t,t}(o,a)=r_{1}(o,a)+\theta_{i}+\lambda_{t}$ for any $i$ and $t$.
This motivates us to consider the following optimization problem:

$\displaystyle(\widehat{\bm{\theta}},\widehat{\bm{\lambda}},\widehat{r}_{1})=\operatorname*{arg\,min}_{\bm{\theta},\bm{\lambda},\,r_{1}}\sum_{i,t}[R_{i,t}-\theta_{i}-\lambda_{t}-r_{1}(O_{i,t},A_{i,t})]^{2},$
(3.3)

where $\bm{\theta}=(\theta_{1},\dots,\theta_{N})^{\top}\in\mathbb{R}^{N}$,
$\bm{\lambda}=(\lambda_{1},\dots,\lambda_{T})^{\top}\in\mathbb{R}^{T}$.

To guarantee the uniqueness of the solution to (3.3), one could impose the
identifiability constraints $\sum_{i}\theta_{i}=\sum_{t}\lambda_{t}=0$. There
are other constraints to consider, but they all lead to the same estimation of
the Q-function. Many supervised learning algorithms can be employed to solve
(3.3). We provide a concrete proposal based on the method of sieves
(Grenander,, 1981) in Section 3.1. Alternatively, one can model $r_{1}$ via
deep neural networks (DNN) and estimate the parameters involved in the DNN as
well as $\\{\theta_{i}\\}_{i}$, $\\{\lambda_{t}\\}_{t}$ via e.g., the Adam
algorithm (Kingma and Ba,, 2015).

We next estimate $\\{Q_{i,t-1,t}\\}_{i,t}$. According to the Bellman equation,
we obtain

$\displaystyle
Q_{i,t-1,t}(O_{i,t-1},A_{i,t-1})={\mathbb{E}}\Big{[}\sum_{a}\pi(a|O_{i,t})Q_{i,t,t}(O_{i,t},a)\big{|}A_{i,t-1},O_{i,t-1}\Big{]}.$

Under the additivity assumption, we can similarly obtain a two-way
decomposition for $Q_{i,t-1,t}$; see Proposition 1 below for a formal
statement. This allows us to solve a similar constrained optimization problem
to (3.3) to estimate $Q_{i,t-1,t}$. We next repeat this procedure to
recursively estimate $\\{Q_{i,t-2,t}\\}_{i,t}$, $\\{Q_{i,t-3,t}\\}_{i,t}$,
$\cdots$, $\\{Q_{i,1,t}\\}_{i,t}$ based on the Bellman equation and finally
construct the policy value estimator via (3.2). We summarize our estimating
procedure in Algorithm 1. The following proposition formally states the two-
way structure of these Q-functions.

######  Proposition 1.

For any integer $k$ such that $1\leq k<t$, the Q-function
$Q^{\pi}_{i,t-k+1,t}(o,a)$ satisfies

$Q^{\pi}_{i,t-k+1,t}(o,a)=r^{\pi}_{k}(o,a)+\theta^{\pi}_{k,i}+\lambda^{\pi}_{k,t},$

where $\theta^{\pi}_{k,i}$ and $\lambda^{\pi}_{k,t}$ depend only on $i,k,\pi$
and $t,k,\pi$ respectively.

In what follows, we will omit the superscript $\pi$ in
$r^{\pi}_{k}(o,a),\mbox{ }\theta^{\pi}_{k,i}$, and $\lambda^{\pi}_{k,t}$ when
there is no confusion.

To conclude this section, we draw a comparison with the classical backward
induction in the DTR literature (see e.g., Murphy,, 2003; Robins,, 2004).
First, the classical backward induction algorithm aims to learn the Q-function
under an optimal policy and derive the optimal policy as the greedy policy
with respect to the estimated Q-function. To the contrary, the proposed
algorithm learns the Q-function under a fixed target policy for the purpose of
policy evaluation.

Second, classical backward induction requires the computation of the
Q-function recursively from time $t$ till the beginning. However, it is worth
mentioning that the Q-function converges exponentially fast to a constant
function with respect to the lag $k$ (see Section C.2 of the supplementary
material). We refer to this phenomenon as Q-function degeneracy. As such,
early stopping can be potentially employed in Algorithm 1 to facilitate the
computation.

Third, the proposed backward induction allows us to to efficiently borrow
information under the additivity assumption. Specifically, during each
iteration, we pull all the relevant data together to estimate the Q-function.
This allows us to consistently estimate the main effect (shared by all
observations) at a rate of $(NT)^{-\alpha}$ which depends on both $N$ and $T$,
and the exponent $0<\alpha\leq 1/2$ depends on the nonparametric methods being
used to solve the constrained optimization problem. Meanwhile, the two-way
fixed effects $\theta$s and $\lambda$s will converge at $T^{-1/2}$ and
$N^{-1/2}$ respectively, up to some logarithmic factors. To the contrary, the
estimator obtained via the classical backward induction algorithm typically
converges at a rate of $N^{-\alpha^{\prime}}$ for some $0<\alpha^{\prime}\leq
1/2$ in individual-homogeneous and history-dependent111The transition function
depends on the entire history instead of the current observation-action pair.
environments.

Algorithm 1 Pseudocodes for Estimating $\eta_{i,t^{*}}^{\pi}$.

1:function ($R_{i,t},A_{i,t},O_{i,t}$ for $1\leq i\leq N,1\leq t\leq T$)

2: Set iteration counter $k\leftarrow 1$.

3: Solve
$(\widehat{\theta}_{i},\widehat{\lambda}_{t},\widehat{r}_{1})=\operatorname*{arg\,min}_{\theta_{i},\lambda_{t},r_{1}}\sum_{i,t}[R_{i,t}-\theta_{i}-\lambda_{t}-r_{1}(O_{i,t},A_{i,t})]^{2}.$

4: Compute
$\widehat{Q}^{\pi}_{i,t,t}(O_{i,t},a)=\widehat{r}_{1}(O_{i,t},a)+\widehat{\theta}_{i}+\widehat{\lambda}_{t}$,
for all $i$ and $t$, and $a\in\mathcal{A}$.

5: repeat

6: $k\leftarrow k+1$

7: For all $i$ and $t\geq k$, solve
$\displaystyle(\widehat{\theta}_{k,i},\widehat{\lambda}_{k,t},\widehat{r}_{k})$
$\displaystyle=\operatorname*{arg\,min}_{\theta_{k,i},\lambda_{k,t},r_{k}}\sum_{i,t}\;\left[\sum_{a\in\mathcal{A}}\pi(a|O_{i,t-k+2})\widehat{Q}^{\pi}_{i,t-k+2,t}(O_{i,t-k+2},a)-\theta_{k,i}-\lambda_{k,t}-r_{k}(O_{i,t-k+1},A_{i,t-k+1})\right]^{2}.$

8: Compute
$\widehat{Q}^{\pi}_{i,t-k+1,t}(O_{i,t-k+1},a)=\widehat{r}_{k}(O_{i,t-k+1},a)+\widehat{\theta}_{k,i}+\widehat{\lambda}_{k,t}$,
for $a\in\mathcal{A}$.

9: until $k=t^{*}$

10:
$\widehat{\eta}^{\pi}_{i,t^{*}}=\sum_{a}\widehat{Q}^{\pi}_{i,1,t^{*}}(O_{i,1},a)\pi(a|O_{i,1})$

###  3.1 A Linear Sieve Estimator for Two-way Fixed Effects Model

Notation. Given arbitrary $\\{x_{i,t}\\}_{1\leq i\leq N,1\leq t\leq T}$, let
$\bm{x}\in\mathbb{R}^{NT}$ denote the vector whose $((t-1)N+i)$-th element
equals $x_{i,t}$. That is, $\bm{x}$ is constructed by stacking the $N$
elements at the first time point, followed by the $N$ elements at the second
time point, and so on, until the $N$ elements at the final time point, e.g.,

$\bm{x}=(x_{1,1},x_{2,1},\dots,x_{N,1},x_{1,2},\dots,x_{N-1,T},x_{N,T})^{\top}.$

Similarly, given a set of vectors $\\{\bm{x}_{i,t}\\}_{i,t}$, let $\bm{X}$
denote the matrix whose $((t-1)N+i)$th row of $\bm{X}$ equals $\bm{x}_{i,t}$.

To implement Algorithm 1, we need to solve two-way fixed effects models
repeatedly for value function estimation. To simplify the presentation, we
focus on the estimation of $Q^{\pi}_{i,t,t}(O_{i,t},a_{i,t})$ (see Equation
(3.3)). We propose to approximate the main effect function $r_{1}(o,a)$ using
linear sieves (Huang,, 1998; Chen and Christensen,, 2015), i.e.,
$r_{1}(o,a)\approx\bm{\Phi}_{L}^{\top}(o)\bm{\beta}_{a}$, where
$\bm{\Phi}_{L}(o)$ is a vector consisting of $L$ sieve basis functions, e.g.,
splines or wavelet bases. Under mild conditions, for spline basis functions,
there exists a vector $\\{\bm{\beta}^{*}_{a}\\}$ such that the approximation
error is negligible, i.e.,
$\sup_{o,a}|r_{1}(o,a)-\bm{\Phi}_{L}(o)^{\top}\bm{\beta}^{*}_{a}|=O(L^{-p/d})$
where $p>0$ measures the smoothness of the system dynamics; see Assumption 1
in Section 5 (Schumaker,, 2007). For simplicity, we now focus on the binary
action space setting, in which $\mathcal{A}=\\{0,1\\}$.

The two-way fixed effects model in (2.3) can be represented in the following
matrix form: $\bm{R}=\bm{B}\bm{\alpha}+\bm{M}+\bm{\varepsilon},$ where
$\bm{R}=(R_{1,1},R_{2,1},\dots,R_{N,1},R_{2,1},\dots,R_{N-1,T},R_{N,T})^{\top}\in\mathbb{R}^{NT}$,
$\bm{\alpha}=(\bm{\theta}^{\top},\bm{\lambda}^{\top})^{\top}$,
$\bm{B}=(\bm{1}_{T}\otimes\bm{I}_{N},\bm{I}_{T}\otimes\bm{1}_{N})\in\mathbb{R}^{NT\times(N+T)}$
is the design matrix, $\bm{I}_{N}$ is a $N\times N$ identity matrix,
$\bm{1}_{T}$ is a vector of length $T$ with all elements one, and $\otimes$ is
the Kronecker product. In what follows, we will omit the indices of these
matrices and vectors when there is no confusion.

Let
$\bm{\Phi}_{i,t}=((1-A_{i,t})\bm{\Phi}_{L}^{\top}(O_{i,t}),A_{i,t}\bm{\Phi}_{L}^{\top}(O_{i,t}))^{\top}$,
and let $\bm{\Phi}$ be the $\mathbb{R}^{NT\times 2L}$ matrix

$(\bm{\Phi}_{1,1}^{\top},\bm{\Phi}_{2,1}^{\top},\dots,\bm{\Phi}_{N,1}^{\top},\bm{\Phi}_{1,2}^{\top},\dots,\bm{\Phi}_{N-1,T}^{\top},\bm{\Phi}_{N,T}^{\top})^{\top}.$

By the Frisch–Waugh–Lovell theorem (Frisch and Waugh,, 1933; Lovell,, 1963), a
closed-form estimator of
$\bm{\beta}=(\bm{\beta}_{0}^{\top},\bm{\beta}_{1}^{\top})^{\top}$ can be
obtained accordingly:

$\displaystyle\widehat{\bm{\beta}}=(\bm{\Phi}^{\top}(\bm{I}-\bm{P})\bm{\Phi})^{-1}\bm{\Phi}^{\top}(\bm{I}-\bm{P})\bm{R},$
(3.4)

where $\bm{P}$ is the projection matrix:
$\bm{P}=\bm{B}(\bm{B}^{\top}\bm{B})^{+}\bm{B}^{\top}$, and
$(\bm{B}^{\top}\bm{B})^{+}$ is the Moore–Penrose inverse of the matrix
$\bm{B}^{\top}\bm{B}$.

Given the estimated $\bm{\beta}$, the estimator of main effect
$r_{1}(O_{i,t},A_{i,t})$ (denoted by $\widehat{r}_{1}$) can also be obtained,
based on which the fixed effects can be calculated. Under the constraints that
$\sum_{i=1}^{N}\theta_{i}=\sum_{t=1}^{T}\lambda_{t}=0$, we have

$\widehat{\theta}_{i}=T^{-1}\sum_{t=1}^{T}(R_{i,t}-\widehat{r}_{1}(O_{i,t},A_{i,t})),\mbox{
and
}\widehat{\lambda}_{t}=N^{-1}\sum_{i=1}^{N}(R_{i,t}-\widehat{r}_{1}(O_{i,t},A_{i,t})).$

The resulting estimated Q-function is given by

$\widehat{Q}^{\pi}_{i,t,t}(O_{i,t},a)=\widehat{\theta}_{i}+\widehat{\lambda}_{t}+\widehat{\bm{\Phi}}_{L}(O_{i,t})^{\top}\widehat{\bm{\beta}}_{a}.$

######  Remark 7.

By design, the matrix $\bm{B}$ is a singular matrix, since the sum of the
first $N$ columns is equal to the sum of the last $T$ columns. Consequently,
the generalized inverse $(\bm{B}^{\top}\bm{B})^{+}$ is used to compute the
matrix $\bm{P}$.

######  Remark 8.

When $N$ or $T$ is large, solving the two-way fixed effects models directly
can be computationally burdensome. An alternative approach is to separate the
estimation procedure into two tasks so that the main effect can be computed
first, followed by the estimation of time and individual-specific effects.
This avoids inverting a large $(N+T+2L)\times(N+T+2L)$ square matrix, making
it computationally efficient (see Algorithm 2 in the Appendix).

##  4 Model-based OPE

In this section, we develop model-based methods that derive the off-policy
value estimator by learning the system dynamics. Recall that under the
additivity assumption,

$\displaystyle
R_{i,t}=\theta_{i}+\lambda_{t}+r_{1}(O_{i,t},A_{i,t})+\varepsilon_{i,t}.$

As we discussed in Section 3, the main effect $r_{1}$, as well as the
individual- and time-specific effects can be estimated by solving the
following optimization,

$\displaystyle\operatorname*{arg\,min}_{\\{\theta_{i}\\}_{i},\\{\lambda_{t}\\}_{t},r_{1}}\sum_{i,t}[R_{i,t}-\theta_{i}-\lambda_{t}-r_{1}(O_{i,t},A_{i,t})]^{2}.$

In addition, we need to estimate the mixing probabilities $\pi_{u}$,
$\pi_{v}$, $\pi_{0}$ as well as the distribution functions
$p_{u_{i}}(o^{\prime}|u_{i})$, $p_{v_{t}}(o^{\prime}|v_{t})$,
$p_{0}(o^{\prime}|a,o)$, obtained by marginalizing over
$p_{u_{i}}(o^{\prime},r|u_{i})$, $p_{v_{t}}(o^{\prime},r|v_{t})$,
$p_{0}(o^{\prime},r|a,o)$ in Equation (2.3).

Given these estimators, we employ simulation-based method to construct the
policy value. To be more specific, based on the estimated transition function,
we can simulate an observation $O_{i,2}^{*}$ based on the observed $O_{i,1}$
under the target policy $\pi$. We next sequentially simulate a sequence of
observations $\\{O_{i,t}^{*}\\}_{t}$ under $\pi$ and compute the estimated
reward
$\pi(a|O_{i,t}^{*})(\widehat{r}_{1}(O_{i,t}^{*},a)+\widehat{\theta}_{i}+\widehat{\lambda}_{t})$.
Finally, we repeat this procedure sufficiently many times and average all the
estimated rewards across different simulations.

Likelihood. It remains to estimate $p_{u}$, $p_{v}$, $p_{0}$ and $\pi_{u}$,
$\pi_{v}$, $\pi_{0}$. Given the latent factors, the likelihood function is
proportional to the following,

$\displaystyle\prod_{i=1}^{N}\prod_{t=2}^{T}p(O_{i,t}|u_{i},v_{t-1},A_{i,t-1},O_{i,t-1};\bm{\Theta})$
$\displaystyle=$
$\displaystyle\prod_{i=1}^{N}\prod_{t=2}^{T}[\pi_{u}p_{u_{i}}(O_{i,t}|u_{i};\bm{\Theta})+\pi_{v}p_{v_{t}}(O_{i,t}|v_{t-1},\bm{\Theta})+\pi_{0}p_{0}(O_{i,t}|A_{i,t-1},O_{i,t-1};\bm{\Theta})],$
(4.1)

where we parameterize the transition model by
$\bm{\Theta}=\\{\pi_{0},\pi_{u},\pi_{v},\bm{\Theta}_{0},\bm{\Theta}_{v},\bm{\Theta}_{u}\\}$,
and $\bm{\Theta}_{0},\bm{\Theta}_{v},\bm{\Theta}_{u}$ are the parameters
corresponding with models $p_{0},p_{v}$ and $p_{u}$ respectively.

We introduce a latent variable $Z_{i,t}\in\\{0,1,2\\}$ such that

$Z_{i,t}=\begin{cases}0,\mbox{ if $O_{i,t}$ is generated by $p_{0}$}.\\\
1,\mbox{ if $O_{i,t}$ is generated by $p_{v}$},\\\ 2,\mbox{
otherwise}.\end{cases}$

However, directly maximizing the likelihood in (4.1) is challenging, since it
requires to marginalize over $Z_{i,t}$. Toward that end, we employ the
expectation-maximization (EM, Dempster et al.,, 1977; Wu,, 1983; Meng and Van
Dyk,, 1997) algorithm for parameter estimation. The EM algorithm is a
recursive algorithm which alternates between an E-step for computing
conditional expectation and an M-step for maximizing the likelihood. It has
been widely used in various statistical problems including missing data,
mixture model and clustering (Little and Rubin,, 2019). We detail the two
steps below.

E-step. Similar to (4.1), the complete log-likelihood involving both
$\\{O_{i,t}\\}_{i,t}$ and $\\{Z_{i,t}\\}_{i,t}$ is given by

$\displaystyle l(\bm{O},\bm{Z}|\bm{U},\bm{V},\bm{A};\bm{\Theta})\propto$
$\displaystyle\sum_{i=1}^{N}\sum_{t=2}^{T}\log
p(O_{i,t}|Z_{i,t},u_{i},v_{t-1},A_{i,t-1},O_{i,t-1};\bm{\Theta})p(Z_{i,t};\bm{\Theta})$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{N}\sum_{t=2}^{T}[\mathbb{I}(Z_{i,t}=0)\log(\pi_{0}p_{0}(O_{i,t}|A_{i,t-1},O_{i,t-1};\bm{\Theta}_{0}))$
$\displaystyle+$
$\displaystyle\mathbb{I}(Z_{i,t}=1)\log(\pi_{v}p_{v_{t}}(O_{i,t}|v_{t};\bm{\Theta}_{v}))+\mathbb{I}(Z_{i,t}=2)\log(\pi_{u}p_{u_{i}}(O_{i,t}|u_{i};\bm{\Theta}_{u}))].$
(4.2)

Given an current guess of $\bm{\Theta}$, say $\widetilde{\bm{\Theta}}$, define
$\Gamma(\bm{\Theta}|\widetilde{\bm{\Theta}})$ as the expected value of the
$l(\bm{O},\bm{Z}|\bm{U},\bm{V},\bm{A};\bm{\Theta})$ with respect to the
currently estimated conditional distribution
$p(\bm{Z}|\bm{U},\bm{V},\bm{O},\bm{A};\widetilde{\bm{\Theta}})$, i.e.,
$\Gamma(\bm{\Theta}|\widetilde{\bm{\Theta}})=\mathbb{E}_{Z\sim
p(\bm{Z}|\bm{U},\bm{V},\bm{O},\bm{A};\widetilde{\bm{\Theta}})}l(\bm{O},\bm{Z}|\bm{U},\bm{V},\bm{A};\bm{\Theta})$.
We aim to calculate $\Gamma$ in this step. It follows from (4) that

$\displaystyle\begin{split}\Gamma(\bm{\Theta}|\widetilde{\bm{\Theta}})=&\sum_{i=1}^{N}\sum_{t=2}^{T}[p(Z_{i,t}=0|O_{i,t},A_{i,t-1},O_{i,t-1};\widetilde{\bm{\Theta}}_{0})\log(\pi_{0}p_{0}(O_{i,t}|A_{i,t-1},O_{i,t-1};\bm{\Theta}_{0}))\\\
&+p(Z_{i,t}=1|O_{i,t},v_{t};\widetilde{\bm{\Theta}}_{v})\log(\pi_{v}p_{v_{t}}(O_{i,t}|v_{t};\bm{\Theta}_{v}))\\\
&+p(Z_{i,t}=2|O_{i,t},u_{i};\widetilde{\bm{\Theta}}_{u})\log(\pi_{u}p_{u_{i}}(O_{i,t}|u_{i};\bm{\Theta}_{u}))].\end{split}$
(4.3)

M-step. We aim to update the model parameter $\bm{\Theta}_{new}$ that
maximizes $\Gamma(\bm{\Theta}|\widetilde{\bm{\Theta}})$ with respect to
$\bm{\Theta}$, i.e.,
$\bm{\Theta}_{new}=\operatorname*{arg\,max}_{\bm{\Theta}}\Gamma(\bm{\Theta}|\widetilde{\bm{\Theta}})$.
It follows from (4.3) that

$\displaystyle\pi_{0,new}$ $\displaystyle=$
$\displaystyle\frac{1}{N(T-1)}\sum_{i=1}^{N}\sum_{t=2}^{T}p(Z_{i,t}=0|O_{i,t},A_{i,t-1},O_{i,t-1};\widetilde{\bm{\Theta}}_{0}),$
$\displaystyle\pi_{v,new}$ $\displaystyle=$
$\displaystyle\frac{1}{N(T-1)}\sum_{i=1}^{N}\sum_{t=2}^{T}p(Z_{i,t}=1|O_{i,t},v_{t};\widetilde{\bm{\Theta}}_{v}),$
$\displaystyle\pi_{u,new}$ $\displaystyle=$
$\displaystyle\frac{1}{N(T-1)}\sum_{i=1}^{N}\sum_{t=2}^{T}p(Z_{i,t}=2|O_{i,t},u_{i};\widetilde{\bm{\Theta}}_{u}).$

The rest of the parameters can be updated using any derivative-based (e.g.,
quasi-Newton) or derivative-free (e.g., Nelder-Mead) algorithm. Our final
estimator is obtained by repeating the E-step and the M-step until
convergence.

Choice of the parametric family. In our implementation, when the observation
is continuous, we posit normal distribution functions for
$p_{u_{i}}(o^{\prime}|u_{i})$, $p_{v_{t}}(o^{\prime}|v_{t})$,
$p_{0}(o^{\prime}|a,o)$, i.e.,
$p_{u_{i}}(o^{\prime}|u_{i})=\phi(o^{\prime};\mu_{u_{i}},\Sigma_{u_{i}})$,
$p_{u_{i}}(o^{\prime}|v_{t})=\phi(o^{\prime};\mu_{v_{t}},\Sigma_{v_{t}})$ and
$p_{0}(o^{\prime}|a,o)=\phi(o^{\prime};\mu_{0}(a,o),\Sigma_{0}(a,o))$ where
$\phi(\bullet;\mu,\Sigma)$ denotes a $d$-dimensional multivariate normal
density function with mean vector $\mu$ and covariance matrix $\Sigma$. We
further use a linear model for the mean function $\mu_{0}$, i.e.,
$\mu_{0}(a,o)=\Lambda o+\psi a$ and a constant model for the covariance
function, i.e., $\Sigma_{0}(a,o)=\Sigma_{0}$ for any $a$ and $o$. As such, the
set of parameters $\bm{\Theta}$ can be summarized by
$\\{\pi_{0},\pi_{v},\pi_{u},\\{\mu_{u_{i}}\\}_{i},\\{\Sigma_{u_{i}}\\}_{i},\\{\mu_{v_{t}}\\}_{t},\\{\Sigma_{v_{t}}\\}_{t},\Lambda,\psi,\Sigma_{0}\\}$.

##  5 Theoretical Results

In this section, we focus on investigating the theoretical properties of our
proposed model-free estimators. Consistencies and convergence rates of the
model-based estimators can be established based on existing analyses of EM
algorithms (see e.g., Wu,, 1983; Balakrishnan et al.,, 2017) and we omit the
details to save space.

We begin with a summary of our theoretical results. Theorem 1 is concerned
with the convergence rates of the proposed value estimators. In particular,
for a sufficiently large $L$, we show that the estimated average effect
$\widehat{\eta}^{\pi}$, individual-specific effect $\widehat{\eta}_{i}^{\pi}$,
time-specific effect $\widehat{\eta}_{t}^{\pi}$ and individual- and time-
specific $\widehat{\eta}_{i,t}^{\pi}$ converge at a rate of $(NT)^{-1/2}$,
$T^{-1/2}$, $N^{-1/2}$ and $\min^{-1/2}(N,T)$, respectively, up to some
logarithmic factors. Theorem 2 establishes the limiting distributions of these
estimators.

We next impose some technical assumptions.

######  Assumption 1 (Hölder smoothness).

Assume that there exists some constants $p$ and $C$, such that for any
$a\in\mathcal{A}$ and $o\in\mathcal{O}$, both the function $r_{1}(\cdot,a)$
and $p_{0}(o^{\prime}|\cdot,a)$ belong to the class of $p$-smooth functions
$\Lambda(p,C)$; see Appendix C.1 for the detailed definition.

######  Assumption 2 (Basis functions).

(i) $\sup_{o}\|\bm{\Phi}_{L}(o)\|_{2}=O(\sqrt{L})$ and
$\lambda_{\max}[\int_{o\in\mathcal{O}}\bm{\Phi}_{L}(o)\bm{\Phi}^{\top}_{L}(o)do]=O(1)$;
(ii) For any $C>0$,
$\sup_{f\in\Lambda(p,C)}\inf_{\beta\in\mathbb{R}^{L}}\sup_{o}|\Phi^{\top}_{L}(o)\beta-f(o)|=O(L^{-p/d})$;
(iii) $L\ll\min(N,T)/\log(NT)$.

######  Assumption 3 (System dynamics).

(i) Assume that there exist random errors $\\{e_{i,t}\\}_{i,t}$ that are i.i.d
copies of $E$ such that the future observation $O_{i,t+1}$ can be represented
as $\kappa(O_{i,t},A_{i,t},u_{i},v_{t},e_{i,t})$ for some function $\kappa$
that satisfies

$\displaystyle\sup_{a,u,v}\mathbb{E}\lVert\kappa(o,a,u,v,E)-\kappa(o^{\prime},a,u,v,E)\rVert_{2}\leq
q\lVert o-o^{\prime}\rVert_{2},$
$\displaystyle\sup_{o,a}\lVert\kappa(o,a,u,v,E)-\kappa(o,a,u,v,E^{\prime})\rVert_{2}=O(\lVert
E-E^{\prime}\rVert_{2}),$

for some $0\leq q<1$; (ii) each element of the error $E$ vector has sub-
exponential tail, i.e., $\max_{j}{\mathbb{E}}\exp(t|E_{j}|)<\infty$ for some
$t>0$, where $E_{j}$ denotes the $j$th element of $E$; (iii) the density
functions $p_{u}$, $p_{v}$ and $N^{-1}\sum_{i=1}^{N}p_{O_{i,1}}$
($p_{O_{i,1}}$ denotes the density function of $O_{i,1}$) are uniformly
bounded.

######  Assumption 4 (Stability).

For any backward step $k$ (the $k$th iteration in Algorithm 1),

$\displaystyle\lambda_{\min}[\mathbb{E}(\bm{\Phi}_{k}^{\top}\bm{S}_{k}\bm{\Phi}_{k-1}^{new})]\geq(NT)\rho_{0}\mbox{
and
}\lVert[\mathbb{E}(\bm{\Phi}_{k}^{\top}\bm{S}_{k}\bm{\Phi}_{k})]^{-1}\mathbb{E}(\bm{\Phi}_{k}^{\top}\bm{S}_{k}\bm{\Phi}_{k-1}^{new})\rVert_{2}\leq\rho_{1},$

for some constants $\rho_{0}>0$ and $0<\rho_{1}<1$, where $\bm{\Phi}_{k}$ are
the matrices consisting of the first $N(T-k+1)$ rows of matrix $\bm{\Phi}$,
$\bm{S}_{k}$ and $\bm{B}_{k}$ is the residual maker matrix, the design matrix
for the fixed effects at step $k$ respectively (similarly, $\bm{S}_{k}$ and
$\bm{B}_{k}$ is the matrix consisting of the first $N(T-k+1)$ rows of matrix
$\bm{S}$ and $\bm{B}$ respectively, see Section C.4.5 of the supplementary
article for the detailed formulation), and $\bm{\Phi}_{k}^{new}$ is the design
matrix such that

$\displaystyle\bm{\Phi}^{new}_{k,i,t}=\frac{1}{T-k+1}\sum_{t=1}^{T-k+1}\bm{\Phi}_{k,i,t}+\frac{1}{N}\sum_{i=1}^{N}\bm{\Phi}_{k,i,t}+\sum_{a\in\mathcal{A}}\pi(a|O_{k,i,t})\bm{\Phi}_{k,i,t}.$
(5.1)

Assumption 1 is frequently imposed in the sieve estimation literature (see
e.g., Huang,, 1998; Chen and Christensen,, 2015). It has been recently
employed in the RL literature to obtain sharp convergence rate for the
estimated Q- and value estimator (Fan et al.,, 2020; Chen and Qi,, 2022; Li et
al.,, 2022; Shi et al.,, 2022). In doubly inhomogeneous environments, it
allows us to show that $Q_{i,t_{1},t_{2}}^{\pi}$ is $p$-smooth for any
$i,t_{1},t_{2}$ and $\pi$; see Lemma 1 in Section C.1 of the supplementary
article.

Assumption 2(i) and (ii) are automatically satisfied when tensor product
B-spline or wavelet bases are employed; see Section 6 of Chen and Christensen,
(2015) for a review of these basis functions. Assumption 2(iii) is to
guarantee the consistency of the sieve estimator.

Assumption 3(i) and (ii) allow us to establish concentration inequalities in
doubly inhomogeneous environments. Under the additivity assumption (2.3),
there exist functions $\kappa_{0}$, $\kappa_{u}$, $\kappa_{v}$ and random
errors $E_{0}$, $E_{u}$ and $E_{v}$ such that
$\kappa(o,a,u,v,E)\stackrel{{\scriptstyle
d}}{{=}}\mathbb{I}(Z=0)\kappa_{0}(o,a,E_{0})+\mathbb{I}(Z=1)\kappa_{v}(v,E_{v})+\mathbb{I}(Z=2)\kappa_{u}(u,E_{u})$
where the latent variable $Z$ is independent of $(E_{0},E_{u},E_{v})$ and
satisfies that ${\mathbb{P}}(Z=0)=\pi_{0}$, ${\mathbb{P}}(Z=1)=\pi_{v}$,
${\mathbb{P}}(Z=2)=\pi_{u}$. As such, Assumption 3 is automatically satisfied
if $E_{0},E_{v}$ and $E_{u}$ have sub-exponential tails, $\kappa_{0}$,
$\kappa_{u}$ and $\kappa_{v}$ are Lipschitz continuous as functions of the
error term and

$\displaystyle\sup_{a}{\mathbb{E}}\|\kappa_{0}(o,a,E_{0})-\kappa_{0}(o^{\prime},a,E_{0})\|_{2}\leq
q\|o-o^{\prime}\|_{2},$ (5.2)

for some $0\leq q<1$. Notice that (5.2) is automatically satisfied for the
auto-regressive model $O^{\prime}=f(O,A)+g(E_{0})$ for any $f$ such that
$\sup_{a}|f(o,a)-f(o^{\prime},a)|\leq q\|o-o^{\prime}\|_{2}$. Other examples
are provided in Diaconis and Freedman, (1999). Notice that Assumption 1
implicitly implies the boundedness of the transition density function $p_{0}$.
This together with Assumption 3(iii) yields the uniform boundedness of
marginal density functions of $\\{O_{i,t}\\}_{i,t}$.

The first part of Assumption 4 essentially requires
$(NT)^{-1}\mathbb{E}(\bm{\Phi}_{k}^{\top}\bm{S}_{k}\bm{\Phi}_{k})$ to be
invertible. The second part is closely related to the irrepresentable or
mutual incoherence condition in the variable selection literature for the
selection consistency of the least absolute shrinkage and selection operator
(Meinshausen and Bühlmann,, 2006; Zhao and Yu,, 2006; Zou,, 2006). This type
of assumption is necessary to ensure the consistency of the subsequent value
estimator (Perdomo et al.,, 2022). Similar assumptions have been imposed in
the statistics literature (Ertefaie and Strawderman,, 2018; Luckett et al.,,
2020; Shi et al.,, 2022).

Finally, we present our theories. Recall that both $\eta_{i}^{\pi}$,
$\eta_{i,t}^{\pi}$ as well as their estimators implicitly depend on the
initial observation $O_{i,1}$. As such, it is more proper to write them as
functions of $O_{i,1}$, e.g.,
$\eta_{i,t}^{\pi}(o)={\mathbb{E}}^{\pi}(R_{i,t}|O_{i,1}=o)$,
$\widehat{\eta}_{i,t}^{\pi}(o)=\sum_{a}\widehat{Q}_{i,1,t}^{\pi}(o,a)\pi(a|o)$
($\eta_{i}^{\pi}(o)$ and $\widehat{\eta}_{i}^{\pi}(o)$ can be similarly
defined). For these effects, instead of considering the differences
$\widehat{\eta}_{i,t}^{\pi}(O_{i,1})-\eta_{i,t}^{\pi}(O_{i,1})$ and
$\widehat{\eta}_{i,t}^{\pi}(O_{i,1})-\eta_{i}^{\pi}(O_{i,1})$, we focus on the
aggregated differences
$\int_{o\in\mathcal{O}}[\widehat{\eta}_{i,t}^{\pi}(o)-\eta_{i,t}^{\pi}(o)]do$
and
$\int_{o\in\mathcal{O}}[\eta_{i}^{\pi}(o)-\widehat{\eta}_{i,t}^{\pi}(o)]do$ to
eliminate the variability due to $O_{i,1}$.

######  Theorem 1 (Rates of Convergence).

Assume Assumptions 1 – 4 hold. Then with probability approaching $1$, we have
for any $1\leq i\leq N$ and $1\leq t\leq T$ that

$\displaystyle\max_{i,t}\left|\int_{o\in\mathcal{O}}[\widehat{\eta}_{i,t}^{\pi}(o)-\eta_{i,t}^{\pi}(o)]do\right|=O(L^{-p/d})+O(N^{-1/2}\sqrt{\log(NT)})+O(T^{-1/2}\sqrt{\log(NT)}),$
$\displaystyle\max_{i}\left|\int_{o\in\mathcal{O}}[\widehat{\eta}_{i}^{\pi}(o)-\eta_{i}^{\pi}(o)]do\right|=O(L^{-p/d})+O(T^{-1/2}\sqrt{\log(NT)}),$
$\displaystyle\max_{t}|\widehat{\eta}_{t}^{\pi}-\eta_{t}^{\pi}|=O(L^{-p/d})+O(N^{-1/2}\sqrt{\log(NT)})\,\,\mbox{and
}\,\,|\widehat{\eta}^{\pi}-\eta^{\pi}|=O(L^{-p/d})+O(\sqrt{\log(NT)/NT}).$

Theorem 1 highlights a noteworthy property of our method: the error bounds of
value estimator depend solely on $p$, $d$, $L$, $N$ and $T$, and are
independent of the number of backward inductions conducted. This is due to an
important feature of our approach: the error term at the $k$th backward stage
is of the order $O(\pi_{0}^{k})$ (as demonstrated by Lemma 2 in the
supplementary article). Specifically, the error bounds for each value
estimator comprise two components: the bias term $O(L^{-p/d})$ and the
variance term $O(N^{-1/2}\sqrt{\log(NT)})$, $O(T^{-1/2}\sqrt{\log(NT)})$ or
$O(\sqrt{\log(NT)/NT})$. Moreover, for sufficiently large $L$, it is evident
that due to aggregation over time and population, the convergence of the
average effect $\widehat{\eta}^{\pi}$ is the fastest, whereas the individual-
and time-specific effect $\widehat{\eta}^{\pi}_{i,t}$ demonstrates the slowest
convergence.

######  Theorem 2 (Asymptotically Normality).

Assume Assumptions 1 – 4 hold, and there exists some constant $c\geq 1$ such
that $\mathbb{E}(\varepsilon_{i,t}^{2}|O_{i,t}=o,A_{i,t}=a)>c^{-1}$, for any
$i$, $t$, $o\in\mathcal{O}$ and $a\in\mathcal{A}$. Furthermore, assume that
the number of basis functions $L$ satisfies $NT\ll L^{2p/d}$. Then when both
$N$ and $T$ goes to infinity,

$\displaystyle\sqrt{\min(N,T)}\sigma_{\eta_{i,t}^{\pi}}^{-1}\int_{o\in\mathcal{O}}\left(\widehat{\eta}_{i,t}-\eta_{i,t}\right)do\xrightarrow{\enskip
d\enskip}\mathcal{N}(0,1)\mbox{,
}\sqrt{T}\sigma_{\eta_{i}^{\pi}}^{-1}\int_{o\in\mathcal{O}}(\widehat{\eta}_{i}^{\pi}-\eta_{i}^{\pi})do\xrightarrow{\enskip
d\enskip}\mathcal{N}(0,1),$
$\displaystyle\sqrt{N}\sigma_{\eta_{t}^{\pi}}^{-1}(\widehat{\eta}_{t}^{\pi}-\eta_{t}^{\pi})\xrightarrow{\enskip
d\enskip}\mathcal{N}(0,1)\mbox{, and
}\sqrt{NT}\sigma_{\eta^{\pi}}^{-1}(\widehat{\eta}^{\pi}-\eta^{\pi})\xrightarrow{\enskip
d\enskip}\mathcal{N}(0,1),$

where
$\sigma_{\eta_{i,t}^{\pi}},\sigma_{\eta_{i}^{\pi}},\sigma_{\eta_{t}^{\pi}}$
and $\sigma_{\eta^{\pi}}$ are some quantities bounded from below and above
(for a detailed formulation, please refer to Section D in the Appendix).

Theorem 2 establishes the asymptotic normality of the value estimators when
both $N$ and $T$ diverges. Specifically, in Section D of the Appendix, we
demonstrate that $\sqrt{\min(N,T)}(\widehat{\eta}_{i,t}-\eta_{i,t})$ can be
expressed as a sum of the martingale difference sequences plus some
asymptotically negligible term, and a direct application of the martingale
difference central limit theorem (McLeish,, 1974) yields the above results.

##  6 Numerical Studies

In this section, we evaluate our proposed model-based and model-free
approaches through extensive simulations. Our focus is on a binary action
space, with a tabular observation space and a continuous observation space. As
discussed earlier, we focus on evaluating the following four targets: the
average reward $\eta^{\pi}$; the $i$th subject’s average reward aggregated
over time $\eta^{\pi}_{i}$; the average reward in the population at time $t$,
denoted by $\eta^{\pi}_{t}$; the $i$th subject’s expected reward at time $t$,
denoted by $\eta^{\pi}_{i,t}$.

Baseline methods. We compare our proposed approach against the following:

  1. 1.

Doubly homogeneous model-free OPE. As commented earlier, the four targets
become asymptotically the same in doubly homogeneous environments. The first
baseline method ignores both types of inhomogeneities and employ the
Q-learning algorithm developed in doubly homogeneous environments to evaluate
the long-term average reward and use the same estimator to estimate
$\eta_{i}^{\pi}$, $\eta_{t}^{\pi}$ and $\eta_{i,t}^{\pi}$. To begin with, we
introduce the following relative value function,

$\displaystyle\widetilde{Q}^{\pi}(o,a)={\mathbb{E}}^{\pi}\left[\sum_{j=t}^{\infty}(R_{j}-\eta^{\pi})|O_{t}=o,A_{t}=a\right],$

where $\eta^{\pi}$ is the average reward by letting $T\to\infty$. Notice that
$\widetilde{Q}^{\pi}$ differs from the proposed Q-function defined in Section
3 under double inhomogeneities. According to the Bellman equation,

$\displaystyle\widetilde{Q}^{\pi}(O_{t},A_{t})={\mathbb{E}}\left[R_{t}-\eta^{\pi}+\sum_{a}\pi(a|O_{t+1})\widetilde{Q}^{\pi}_{t+1}(O_{t+1},a)|O_{t},A_{t}\right].$

We also use linear sieves $\bm{\beta}_{A_{t}}^{\top}\bm{\Phi}_{L}(O_{t})$ to
approximate the Q-function. According to the Bellman equation, we estimate the
average reward $\eta^{\pi}$, as well as the parameters $\bm{\beta}_{0}$,
$\bm{\beta}_{1}$ by solving the following unbiased estimating function

$\displaystyle\begin{bmatrix}\sum_{i,t}\left[R_{i,t}-\eta^{\pi}+\sum_{a}\pi(a|O_{i,t+1})\bm{\beta}_{a}^{\top}\bm{\Phi}_{L}(O_{i,t+1})-\bm{\beta}_{A_{i,t}}^{\top}\bm{\Phi}_{L}(O_{i,t})\right]\\\
\sum_{i,t}A_{i,t}\bm{\Phi}_{L}(O_{i,t})\left[R_{i,t}-\eta^{\pi}+\sum_{a}\pi(a|O_{i,t+1})\bm{\beta}_{a}^{\top}\bm{\Phi}_{L}(O_{i,t+1})-\bm{\beta}_{1}^{\top}\bm{\Phi}_{L}(O_{i,t})\right]\\\
\sum_{i,t}(1-A_{i,t})\bm{\Phi}_{L}(O_{i,t})\left[R_{i,t}-\eta^{\pi}+\sum_{a}\pi(a|O_{i,t+1})\bm{\beta}_{a}^{\top}\bm{\Phi}_{L}(O_{i,t+1})-\bm{\beta}_{0}^{\top}\bm{\Phi}_{L}(O_{i,t})\right]\\\
\end{bmatrix}=\bm{0}.$

  2. 2.

Temporal stationary model-free OPE. The second baseline method ignores
temporal non-stationary. It applies the first baseline method to each
individual data trajectory to estimate $\eta_{i}^{\pi}$ for $i=1,2,\cdots,N$.
Denote the resulting estimator by $\widehat{\eta}_{i}^{\pi}$. Under temporal
stationarity, it sets
$\widehat{\eta}^{\pi}=\widehat{\eta}^{\pi}_{t}=\sum_{i=1}^{N}\widehat{\eta}^{\pi}_{i}/N$,
and $\widehat{\eta}^{\pi}_{i}=\widehat{\eta}^{\pi}_{i,t}$ for any $t$.

  3. 3.

Individual homogeneous model-free OPE. The third baseline method ignores
individual heterogeneity and adopts the classical backward induction algorithm
developed in the DTR literature to evaluate $\eta_{t}^{\pi}$. Denote the
resulting estimator by $\widehat{\eta}_{t}^{\pi}$. Under individual
homogeneity, it sets
$\widehat{\eta}^{\pi}=\widehat{\eta}^{\pi}_{i}=\sum_{t=1}^{T}\widehat{\eta}^{\pi}_{t}/T$,
and $\widehat{\eta}^{\pi}_{t}=\widehat{\eta}^{\pi}_{i,t}$ for any $i$.

  4. 4.

Doubly homogeneous model-based OPE. The last baseline method ignores both
types of inhomogeneities and adopts the standard model-based approach which
estimates the transition and reward functions, and constructs the resulting
value estimator based on Monte Carlo simulations as in Section 4.

In what follows, we refer to these baseline methods as B1, B2, B3, and B4. The
mean squared error (MSE) and the mean absolute error (MAE) are used as the
evaluation metrics, the true value function is calculated using 500 Monte
Carlo experiments.

Tabular setting. We first consider a tabular setting in which both the
observation and action spaces are binary, with $N=T=80$. The data generation
procedure is as follows. We set $\pi_{0}=0.6$ and $\pi_{u}=\pi_{v}=0.2$. The
transition functions are given by $p_{u_{i}}(o^{\prime}=1|u_{i})=|\sin(i)|$,
$p_{v_{t}}(o^{\prime}=1|v_{t})=|\cos(t)|$,
$p_{0}(o^{\prime}=1|o,a)=0.3\mathbb{I}(o+a\geq 1)+0.7\mathbb{I}(o+a<1)$, where
$\mathbb{I}(\cdot)$ is the indicator function. In addition, we set
$p_{u_{i}}(r|u_{i})=\phi(4\sin(i),1)$,
$p_{v_{t}}(r|v_{t})=\phi(3|\cos(2.4t)|,1)$ and $p_{0}(r|a,o)=\phi(2o+3a,1)$.
The behavior policy is $b(a=1|o)=0.3\mathbb{I}(o=0)+0.7\mathbb{I}(o=1)$; the
target policy we aim to evaluate is an observation-agnostic policy:
$\pi(a=1)=0.8$.

Figure 1: MAE and MSE of the estimated value (four targets) in the tabular
setting using our proposed two-way doubly inhomogeneous decision process model
(TWDIDP1 and TWDIDP2, model-free and model-based respectively), and four
baseline methods B1, B2, B3, and B4.

Figure 1 summarizes the MAE and the MSE of the six different estimators of
$\eta^{\pi}$, $\eta^{\pi}_{t}$, $\eta^{\pi}_{i}$, and $\eta^{\pi}_{i,t}$,
respectively. We make the following observations: (i) our proposed model-free
method outperforms other methods in all settings; (ii) our proposed model-
based method achieves a smaller MSE than other competing methods when
considering the last three estimands; (iii) the doubly homogeneous model-based
method B4 performs the worst across all scenarios.

Figure 2: MAE and MSE of the estimated value (four targets) in the continuous
observation setting using our proposed two-way doubly inhomogeneous decision
process model (TWDIDP1 and TWDIDP2, model-free and model-based respectively),
and four baseline methods B1, B2, B3, and B4.

Continuous observation space. We design a setting with a continuous
observation space. Similar to the tabular setting, we set $N=T=80$,
$\pi_{0}=0.6$ and $\pi_{u}=\pi_{v}=0.2$. The transition function consists of
$p_{u_{i}}(o^{\prime}|u_{i})=\phi(\sin(3i),1)$,
$p_{v_{t}}(o^{\prime}|v_{t})=\phi(\cos(-1.8t),1)$,
$p_{0}(o^{\prime}|o,a)=\phi(-0.25o+a,1)$. Finally, we set
$p_{u_{i}}(r|u_{i})=\phi(4\sin(i),1)$,
$p_{v_{t}}(r|v_{t})=\phi(3|\cos(2.4t)|,1)$, and $p_{0}(r|a,o)=\phi(2o+3a,1)$.
The behavior policy is $b(a=1)=0.5$; the target policy we aim to evaluate is
also a randomized policy: $\pi(a=1)=0.8$. A third order degree polynomial two-
way fixed effects model is used to estimate the Q-function.

Figure 2 presents the MAE and the MSE of the six different estimators of
$\eta^{\pi}$, $\eta^{\pi}_{t}$, $\eta^{\pi}_{i}$, and $\eta^{\pi}_{i,t}$,
respectively. Similar to the tabular setting, our proposed model-free method
outperforms other methods in all cases. In addition, the difference between
our proposed approach and the baseline methods is more obvious, compared to
the tabular setting. Our proposed model-based method also performs better than
the baseline methods in general.

##  7 Real Data Analysis

In this section, we apply our proposed method to a sepsis dataset from MIMIC-
III (Johnson et al.,, 2016), a database that contains information on critical
care patients from Beth Israel Deaconess Medical Center in Boston,
Massachusetts. As mentioned earlier, the heterogeneity in patient’s response
to treatment (Evans et al.,, 2021), along with a potential non-stationary
environments makes it difficult to consistently assess the impact of
conducting a given target policy on patient outcomes.

We focus on a subset of patients who received treatments sequentially over 20
stages. The primary outcome in this analysis is the sequential organ failure
assessment (SOFA) score (Jones et al.,, 2009), which is used to monitor the
progression of organ failure over time and measure the degree of organ
dysfunction or failure in critically ill patients. A higher SOFA score
indicates a higher risk of mortality. At any time point $t$, we consider a
binary treatment $A_{t}\in\\{0,1\\}$ where $A_{t}=1$ indicates that the
patient received an intravenous fluid intervention with a dose greater than
the median value for the group of patients being studied, and $A_{t}=0$
otherwise. In previous studies, Raghu et al., (2017) and Zhou et al., 2022b
examined joint action spaces with both vasopressors and intravenous fluid
interventions. We focus solely on intravenous fluid intervention due to the
findings of Zhou et al., 2022b , which suggest a limited impact of
vasopressors.

The following five covariates are included in the observation: gender, age,
the Elixhauser comorbidity index, weight, and the systemic inflammatory
response syndrome score.

Figure 3: The estimated value function of $\eta^{\pi}_{i}$ and
$\eta^{\pi}_{t}$ under three target policies, where $\pi_{1}$ is the always
high dose policy, $\pi_{2}$ is always low dose policy, and $\pi_{3}$ is the
tailored policy.

Three deterministic policies were evaluated using our proposed methods: (i)
always administering a high dose (${\mathbb{P}}(A=1)=1$), (ii) always
administering a low dose (${\mathbb{P}}(A=1)=0$), and (iii) administering a
low dose when the SOFA score is less than 11, and a high dose otherwise. The
third policy is tailored to the SOFA score, taking into account evidence that
a SOFA score of more than 11 is associated with a $100\%$ mortality rate
(Jones et al.,, 2009). To estimate the Q-function, we employed a third-order
degree polynomial two-way fixed effects model at each time point. The average
value estimators for the three policies are as follows: 7.26 (always high
dose), 6.85 (always low dose), and 6.51 (tailored by SOFA score). These
results indicate that the tailored policy is the most effective policy as it
yields the lowest estimated SOFA score. Figure 3 summarizes the estimated
$\eta^{\pi}_{i}$s and $\eta^{\pi}_{t}$s, clearly demonstrating that the
tailored policy outperforms the other two policies, while the always high dose
policy performs the poorest.Our conclusion is in line with these existing
results, which recommend the low dose policy over the high dose policy. It is
also consistent with physicians’ recommendations in the behavior data.

##  8 Discussion

In this study, we introduce a novel class of latent factor models for the
reward and observation transition functions. Our research also encompasses the
development of a comprehensive OPE framework, combining both model-based and
model-free approaches. To our knowledge, this is the first paper that develops
statistically sound OPE methods in offline RL with double inhomogeneities. Our
proposal relies on an additivity assumption. Without this assumption, one can
employ model-based methods to simultaneously estimate the transition function
and the random effects. Similar architectures have been considered in the
literature on factor models and multi-relational learning. We leave it for
future research.

## References

  * Angrist and Pischke, (2009) Angrist, J. D. and Pischke, J.-S. (2009).  Mostly harmless econometrics: An empiricist’s companion.  Princeton university press.
  * Arkhangelsky et al., (2021) Arkhangelsky, D., Imbens, G. W., Lei, L., and Luo, X. (2021).  Double-robust two-way-fixed-effects regression for panel data.  arXiv preprint arXiv:2107.13737.
  * Bai and Ng, (2002) Bai, J. and Ng, S. (2002).  Determining the number of factors in approximate factor models.  Econometrica, 70(1):191–221.
  * Balakrishnan et al., (2017) Balakrishnan, S., Wainwright, M. J., and Yu, B. (2017).  Statistical guarantees for the em algorithm: From population to sample-based analysis.  The Annals of Statistics, 45(1):77–120.
  * Bather, (2000) Bather, J. (2000).  Decision theory: An introduction to dynamic programming and sequential decisions.  John Wiley & Sons, Inc.
  * Bellman, (1957) Bellman, R. (1957).  A markovian decision process.  Journal of Mathematics and Mechanics, 6(5):679–684.
  * Bibaut et al., (2021) Bibaut, A., Petersen, M., Vlassis, N., Dimakopoulou, M., and van der Laan, M. (2021).  Sequential causal inference in a single world of connected units.  arXiv preprint arXiv:2101.07380.
  * Bollen, (2002) Bollen, K. A. (2002).  Latent variables in psychology and the social sciences.  Annual review of psychology, 53(1):605–634.
  * Bradley, (2005) Bradley, R. C. (2005).  Basic properties of strong mixing conditions. a survey and some open questions.  Probability surveys, 2:107–144.
  * Chakraborty and Moodie, (2013) Chakraborty, B. and Moodie, E. E. M. (2013).  Statistical methods for dynamic treatment regimes.  Springer, New York.
  * (11) Chen, E. Y., Jordan, M. I., and Li, S. (2022a).  Transferred q-learning.  arXiv preprint arXiv:2202.04709.
  * (12) Chen, E. Y., Song, R., and Jordan, M. I. (2022b).  Reinforcement learning with heterogeneous data: estimation and inference.  arXiv preprint arXiv:2202.00088.
  * Chen and Christensen, (2015) Chen, X. and Christensen, T. M. (2015).  Optimal uniform convergence rates and asymptotic normality for series estimators under weak dependence and weak conditions.  Journal of Econometrics, 188(2):447–465.
  * Chen and Qi, (2022) Chen, X. and Qi, Z. (2022).  On well-posedness and minimax optimal rates of nonparametric q-function estimation in off-policy evaluation.  In International Conference on Machine Learning, pages 3558–3582. PMLR.
  * De Chaisemartin and d’Haultfoeuille, (2020) De Chaisemartin, C. and d’Haultfoeuille, X. (2020).  Two-way fixed effects estimators with heterogeneous treatment effects.  American Economic Review, 110(9):2964–2996.
  * Dempster et al., (1977) Dempster, A. P., Laird, N. M., and Rubin, D. B. (1977).  Maximum likelihood from incomplete data via the em algorithm.  Journal of the Royal Statistical Society: Series B (Methodological), 39(1):1–22.
  * Diaconis and Freedman, (1999) Diaconis, P. and Freedman, D. (1999).  Iterated random functions.  SIAM review, 41(1):45–76.
  * Ertefaie and Strawderman, (2018) Ertefaie, A. and Strawderman, R. L. (2018).  Constructing dynamic treatment regimes over indefinite time horizons.  Biometrika, 105(4):963–977.
  * Evans et al., (2021) Evans, L., Rhodes, A., Alhazzani, W., Antonelli, M., Coopersmith, C. M., French, C., Machado, F. R., Mcintyre, L., Ostermann, M., Prescott, H. C., et al. (2021).  Surviving sepsis campaign: international guidelines for management of sepsis and septic shock 2021.  Intensive care medicine, 47(11):1181–1247.
  * Fan et al., (2020) Fan, J., Wang, Z., Xie, Y., and Yang, Z. (2020).  A theoretical analysis of deep q-learning.  In Learning for Dynamics and Control, pages 486–489. PMLR.
  * Frisch and Waugh, (1933) Frisch, R. and Waugh, F. V. (1933).  Partial time regressions as compared with individual trends.  Econometrica: Journal of the Econometric Society, pages 387–401.
  * Green et al., (2001) Green, D. P., Kim, S. Y., and Yoon, D. H. (2001).  Dirty pool.  International Organization, 55(2):441–468.
  * Grenander, (1981) Grenander, U. (1981).  Abstract inference.  Wiley Series, New York.
  * Hao et al., (2021) Hao, B., Ji, X., Duan, Y., Lu, H., Szepesvári, C., and Wang, M. (2021).  Bootstrapping fitted q-evaluation for off-policy inference.  In International Conference on Machine Learning, pages 4074–4084. PMLR.
  * Hu et al., (2022) Hu, L., Li, M., Shi, C., Wu, Z., and Fryzlewicz, P. (2022).  Doubly inhomogeneous reinforcement learning.  arXiv preprint arXiv:2211.03983.
  * Huang, (1998) Huang, J. Z. (1998).  Projection estimation in multiple regression with application to functional anova models.  The Annals of statistics, 26(1):242–272.
  * Imai and Kim, (2019) Imai, K. and Kim, I. S. (2019).  When should we use unit fixed effects regression models for causal inference with longitudinal data?  American Journal of Political Science, 63(2):467–490.
  * Imai and Kim, (2021) Imai, K. and Kim, I. S. (2021).  On the use of two-way fixed effects regression models for causal inference with panel data.  Political Analysis, 29(3):405–415.
  * Jiang and Li, (2016) Jiang, N. and Li, L. (2016).  Doubly robust off-policy value evaluation for reinforcement learning.  In International Conference on Machine Learning, pages 652–661. PMLR.
  * Johnson et al., (2016) Johnson, A. E. et al. (2016).  Mimic-iii, a freely accessible critical care database.  Scientific data, 3(1):1–9.
  * Jones et al., (2009) Jones, A. E., Trzeciak, S., and Kline, J. A. (2009).  The sequential organ failure assessment score for predicting outcome in patients with severe sepsis and evidence of hypoperfusion at the time of emergency department presentation.  Critical care medicine, 37(5):1649.
  * Kallus and Uehara, (2022) Kallus, N. and Uehara, M. (2022).  Efficiently breaking the curse of horizon in off-policy evaluation with double reinforcement learning.  Operations Research, 70(6):3282–3302.
  * Kingma and Ba, (2015) Kingma, D. and Ba, J. (2015).  Adam: A method for stochastic optimization.  In Proceedings of the 3rd International Conference for Learning Representations.
  * Kosorok and Laber, (2019) Kosorok, M. R. and Laber, E. B. (2019).  Precision medicine.  Annual review of statistics and its application, 6:263–286.
  * Kosorok and Moodie, (2015) Kosorok, M. R. and Moodie, E. E. M. (2015).  Adaptive treatment strategies in practice: planning trials and analyzing data for personalized medicine.  Society for Industrial and Applied Mathematics, Philadelphia, PA.
  * Levine et al., (2020) Levine, S., Kumar, A., Tucker, G., and Fu, J. (2020).  Offline reinforcement learning: Tutorial, review, and perspectives on open problems.  arXiv preprint arXiv:2005.01643.
  * Li et al., (2023) Li, G., Cai, C., Chen, Y., Wei, Y., and Chi, Y. (2023).  Is q-learning minimax optimal? a tight sample complexity analysis.  Operations Research.
  * Li et al., (2022) Li, M., Shi, C., Wu, Z., and Fryzlewicz, P. (2022).  Reinforcement learning in possibly nonstationary environments.  arXiv preprint arXiv:2203.01707.
  * Liao et al., (2021) Liao, P., Klasnja, P., and Murphy, S. (2021).  Off-policy estimation of long-term average outcomes with applications to mobile health.  Journal of the American Statistical Association, 116(533):382–391.
  * Liao et al., (2022) Liao, P., Qi, Z., Wan, R., Klasnja, P., and Murphy, S. A. (2022).  Batch policy learning in average reward markov decision processes.  The Annals of Statistics, 50(6):3364–3387.
  * Little and Rubin, (2019) Little, R. J. and Rubin, D. B. (2019).  Statistical analysis with missing data, volume 793.  John Wiley & Sons.
  * Liu et al., (2018) Liu, Q., Li, L., Tang, Z., and Zhou, D. (2018).  Breaking the curse of horizon: Infinite-horizon off-policy estimation.  Advances in Neural Information Processing Systems, 31.
  * Lovell, (1963) Lovell, M. C. (1963).  Seasonal adjustment of economic time series and multiple regression analysis.  Journal of the American Statistical Association, 58(304):993–1010.
  * Luckett et al., (2020) Luckett, D. J., Laber, E. B., Kahkoska, A. R., Maahs, D. M., Mayer-Davis, E., and Kosorok, M. R. (2020).  Estimating dynamic treatment regimes in mobile health using v-learning.  Journal of the American Statistical Association, 115(530):692–706.
  * McLeish, (1974) McLeish, D. L. (1974).  Dependent central limit theorems and invariance principles.  the Annals of Probability, 2(4):620–628.
  * Meinshausen and Bühlmann, (2006) Meinshausen, N. and Bühlmann, P. (2006).  High-dimensional graphs and variable selection with the lasso.  The Annals of Statistics, 34(3):1436–1462.
  * Meng and Van Dyk, (1997) Meng, X.-L. and Van Dyk, D. (1997).  The em algorithm—an old folk-song sung to a fast new tune.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 59(3):511–567.
  * Murphy, (2003) Murphy, S. A. (2003).  Optimal dynamic treatment regimes.  Journal of the Royal Statistical Society: Series B, 65(2):331–355.
  * Nickel et al., (2011) Nickel, M., Tresp, V., Kriegel, H.-P., et al. (2011).  A three-way model for collective learning on multi-relational data.  In Icml, volume 11, pages 3104482–3104584.
  * Perdomo et al., (2022) Perdomo, J. C., Krishnamurthy, A., Bartlett, P., and Kakade, S. (2022).  A sharp characterization of linear estimators for offline policy evaluation.  arXiv preprint arXiv:2203.04236.
  * Puterman, (2014) Puterman, M. L. (2014).  Markov decision processes: discrete stochastic dynamic programming.  John Wiley & Sons.
  * Qi et al., (2020) Qi, Z., Liu, D., Fu, H., and Liu, Y. (2020).  Multi-armed angle-based direct learning for estimating optimal individualized treatment rules with various outcomes.  Journal of the American Statistical Association, 115(530):678–691.
  * Qian and Murphy, (2011) Qian, M. and Murphy, S. A. (2011).  Performance guarantees for individualized treatment rules.  Annals of statistics, 39(2):1180.
  * Raghu et al., (2017) Raghu, A., Komorowski, M., Ahmed, I., Celi, L., Szolovits, P., and Ghassemi, M. (2017).  Deep reinforcement learning for sepsis treatment.  arXiv preprint arXiv:1711.09602.
  * Robins, (2004) Robins, J. M. (2004).  Optimal structural nested models for optimal sequential decisions.  In Lin, D. Y. and Heagerty, P., editors, Proceedings of the Second Seattle Symposium in Biostatistics, pages 189–326, New York. Springer.
  * Ross, (1976) Ross, S. (1976).  The arbitrage theory of capital asset pricing.  Journal of Economic Theory, 13(3):341–360.
  * Schumaker, (2007) Schumaker, L. (2007).  Spline functions: basic theory.  Cambridge university press.
  * Shi et al., (2022) Shi, C., Zhang, S., Lu, W., and Song, R. (2022).  Statistical inference of the value function for reinforcement learning in infinite-horizon settings.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 84(3):765–793.
  * Silver et al., (2016) Silver, D., Huang, A., et al. (2016).  Mastering the game of go with deep neural networks and tree search.  nature, 529(7587):484–489.
  * Singer et al., (2016) Singer, M. et al. (2016).  The third international consensus definitions for sepsis and septic shock (sepsis-3).  Jama, 315(8):801–810.
  * Sutton and Barto, (2018) Sutton, R. S. and Barto, A. G. (2018).  Reinforcement learning: An introduction.  MIT press.
  * Thomas and Brunskill, (2016) Thomas, P. and Brunskill, E. (2016).  Data-efficient off-policy policy evaluation for reinforcement learning.  In International Conference on Machine Learning, pages 2139–2148. PMLR.
  * Tropp, (2011) Tropp, J. (2011).  Freedman’s inequality for matrix martingales.  Electronic Communications in Probability.
  * Tsiatis et al., (2019) Tsiatis, A. A., Davidian, M., Holloway, S. T., and Laber, E. B. (2019).  Dynamic Treatment Regimes: Statistical Methods for Precision Medicine.  CRC press, Boca Raton, FL.
  * Uehara et al., (2022) Uehara, M., Shi, C., and Kallus, N. (2022).  A review of off-policy evaluation in reinforcement learning.  arXiv preprint arXiv:2212.06355.
  * Wang et al., (2021) Wang, J., Qi, Z., and Wong, R. K. (2021).  Projected state-action balancing weights for offline reinforcement learning.  arXiv preprint arXiv:2109.04640.
  * Wooldridge, (2010) Wooldridge, J. M. (2010).  Econometric analysis of cross section and panel data.  MIT press.
  * Wu, (1983) Wu, C. J. (1983).  On the convergence properties of the em algorithm.  The Annals of statistics, pages 95–103.
  * Xu et al., (2018) Xu, Z., Li, Z., Guan, Q., Zhang, D., Li, Q., Nan, J., Liu, C., Bian, W., and Ye, J. (2018).  Large-scale order dispatch in on-demand ride-hailing platforms: A learning and planning approach.  In Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pages 905–913.
  * Yin and Wang, (2020) Yin, M. and Wang, Y.-X. (2020).  Asymptotically efficient off-policy evaluation for tabular reinforcement learning.  In International Conference on Artificial Intelligence and Statistics, pages 3948–3958. PMLR.
  * Zhang et al., (2012) Zhang, B., Tsiatis, A. A., Laber, E. B., and Davidian, M. (2012).  A robust method for estimating optimal treatment regimes.  Biometrics, 68(4):1010–1018.
  * Zhao and Yu, (2006) Zhao, P. and Yu, B. (2006).  On model selection consistency of lasso.  The Journal of Machine Learning Research, 7:2541–2563.
  * (73) Zhou, W., Zhu, R., and Qu, A. (2022a).  Estimating optimal infinite horizon dynamic treatment regimes via pt-learning.  Journal of the American Statistical Association, accepted.
  * (74) Zhou, Y., Qi, Z., Shi, C., and Li, L. (2022b).  Optimizing pessimism in dynamic treatment regimes: A bayesian learning approach.  arXiv preprint arXiv:2210.14420.
  * Zhu et al., (2017) Zhu, R., Zhao, Y.-Q., Chen, G., Ma, S., and Zhao, H. (2017).  Greedy outcome weighted tree learning of optimal personalized treatment rules.  Biometrics, 73(2):391–400.
  * Zou, (2006) Zou, H. (2006).  The adaptive lasso and its oracle properties.  Journal of the American statistical association, 101(476):1418–1429.

Supplemental Materials for “Off-Policy Evaluation in Doubly Inhomogeneous
Environments”

Throughout this Supplemental Material, we use $c$, $c_{1}$, $c_{2}$, $C$ and
$C^{\prime}$ to denote a generic constant that can vary from line to line. In
addition, given arbitrary $\\{x_{i,t}\\}_{1\leq i\leq N,1\leq t\leq T}$, we
construct the vector $\bm{x}\in\mathbb{R}^{NT}$ such that its $((t-1)N+i)$-th
element is $x_{i,t}$. That is, $\bm{x}$ is constructed by stacking the $N$
elements of the first time point, followed by the $N$ elements of the second
time point, and so on, until the $N$ elements of the final time point. Thus,
we can express $\bm{x}$ as

$\bm{x}=(x_{1,1},x_{2,1},\dots,x_{N,1},x_{1,2},\dots,x_{N-1,T},x_{N,T})^{\top}.$

When it comes to a matrix $\bm{X}\in\mathbb{R}^{NT\times L}$, we stack it in a
similar way, i.e., the $((t-1)N+i)$th row of $\bm{X}$ is $\bm{x}_{i,t}$. In
addition, denote $n=NT$, $T_{k}=T-k+1$ and $n_{k}=NT_{k}$. Finally, we denote
$k$-th iteration quantity as $x_{k,i,t}$, and we omit $k$ if $k=1$, to
simplify the notation. To ease notation, we sometimes omit arguments for
functions of multiple arguments, e.g., $r_{1}=r_{1}(o,a)$, when there is no
confusion.

The rest of the supplementary article is organized as follows. We present
Algorithm 2 in Section A. Section B is devoted to the proof of Proposition 1.
We prove Theorem 1 in Section C. Finally, we prove Theorem 2 in Section D.

##  Appendix A Algorithm 2

Algorithm 2 Profile Least Squares Estimation Procedure.

1:function ($r_{1,i,t},A_{i,t},O_{i,t}$)

2: Set iteration counter $k\leftarrow 1$.

3: Initialize: $\bm{\theta}^{k}=\bm{\lambda}^{k}=\bm{0}$.

4: repeat

5: Update $\bm{\beta}$ and $\bm{M}$:
$\bm{\beta}^{k}=(\bm{\Phi}^{\top}\bm{S}\bm{\Phi})^{-1}\bm{\Phi}^{\top}(\bm{R}-\bm{B}\bm{\alpha}^{k}),\mbox{
and }\bm{M}^{k}=\bm{\Phi}\bm{\beta}^{k}.$

6: Update $\bm{\theta}$ and $\bm{\lambda}$: for all $i$ and $t$,
$\theta_{i}^{k}=T^{-1}\sum_{t=1}^{T}(R_{i,t}-\bm{M}^{k}_{i,t}),\mbox{ and
}\lambda_{t}^{k}=N^{-1}\sum_{i=1}^{N}(r_{1,i,t}-\bm{M}^{k}_{i,t}).$

7: $k\leftarrow k+1$

8: until convergence criterion is satisfied

##  Appendix B Proof of Proposition 1

Proof: We note that $\\{\theta_{i}\\}_{i}$ and $\\{\lambda_{t}\\}_{t}$ are not
uniquely defined. To simplify the proof, we consider one particular
parameterization, given in Equation (B.2).

We start with $k=1$. Consider the Q-function $Q^{\pi}_{i,t-1,t}$. According to
the Bellman equation, we obtain that

$\displaystyle
Q^{\pi}_{i,t-1,t}(o_{i,t-1},a_{i,t-1})=\sum_{a}{\mathbb{E}}[\pi(a|O_{i,t})Q^{\pi}_{i,t,t}(O_{i,t},a)|a_{i,t-1},o_{i,t-1}]$
$\displaystyle=$
$\displaystyle\sum_{a}{\mathbb{E}}[\pi(a|o_{i,t})(r_{1}(o_{i,t},a)+\theta_{i}+\lambda_{t})|a_{i,t-1},o_{i,t-1}]$
$\displaystyle=$
$\displaystyle\theta_{i}+\lambda_{t}+\int_{o}\sum_{a}\left[\pi(a|o)r_{1}(o,a)\left(\pi_{u}p_{u_{i}}(o)do+\pi_{v}p_{v_{t}}(o)+\pi_{0}p_{0}(o|a_{t-1},o_{t-1})\right)\right]do$
$\displaystyle=$
$\displaystyle\underbrace{\theta_{i}+\pi_{u}\int_{o}\sum_{a}\pi(a|o)r_{1}(o,a)p_{u_{i}}(o)do}_{\theta_{2,i}}+\underbrace{\lambda_{t}+\pi_{v}\int_{o}\sum_{a}\pi(a|o)r_{1}(o,a)p_{v_{t}}(o)do}_{\lambda_{2,t}}$
$\displaystyle+$
$\displaystyle\underbrace{\pi_{0}\int_{o}\sum_{a}\pi(a|o)r_{1}(o,a)p_{0}(o|a_{t-1},o_{t-1})do}_{r_{2}(o_{t-1},a_{t-1})},$

where $p_{u_{i}}(o)=\int p_{u_{i}}(o^{\prime},r|u_{i})dr$, and $p_{v_{t}}(o)$
and $p_{0}(o|\cdot,\cdot)$ is similarly defined. Now let

$\widetilde{r}_{2}(o_{t-1},a_{t-1})=\int_{o}\sum_{a}\pi(a|o)r_{1}(o,a)p_{0}(o|a_{t-1},o_{t-1})do,$

it follows that

$r_{2}(o_{t-1},a_{t-1})=\pi_{0}\widetilde{r}_{2}(o_{t-1},a_{t-1}).$

Using the same argument, we have

$\displaystyle
Q^{\pi}_{i,t-2,t}(O_{i,t-2},A_{i,t-2})=\sum_{a}{\mathbb{E}}[\pi(a|O_{i,t})Q^{\pi}_{i,t-1,t}(O_{i,t-1},a)|A_{i,t-2},O_{i,t-2}]$
$\displaystyle=$
$\displaystyle\theta_{2,i}+\pi_{u}\int_{o}\sum_{a}\pi(a|o)r_{2}(o,a)p_{u_{i}}(o)do+\lambda_{2,t}+\pi_{v}\int_{o}\sum_{a}\pi(a|o)r_{2}(o,a)p_{v_{t}}(o)do$
$\displaystyle+$
$\displaystyle\pi_{0}\int_{o}\sum_{a}\pi(a|o)r_{2}(o,a)p_{0}(o|a_{t-2},o_{t-2})do$
$\displaystyle=$
$\displaystyle\underbrace{\theta_{2,i}+\pi_{u}\pi_{0}\int_{o}\sum_{a}\pi(a|o)\widetilde{r}_{2}(o,a)p_{u_{i}}(o)do}_{\theta_{3,i}}+\underbrace{\lambda_{2,t}+\pi_{v}\pi_{0}\int_{o}\sum_{a}\pi(a|o)\widetilde{r}_{2}(o,a)p_{v_{t}}(o)do}_{\lambda_{3,t}}$
$\displaystyle+$
$\displaystyle\underbrace{\pi_{0}^{2}\int_{o}\sum_{a}\pi(a|o)\widetilde{r}_{2}(o,a)p_{0}(o|a_{t-2},o_{t-2})do}_{r_{3}(o_{t-2},a_{t-2})}.$

Similarly, define

$\displaystyle\widetilde{r}_{k}(o_{t-k+1},a_{t-k+1})=\int_{o}\sum_{a}\pi(a|o)\widetilde{r}_{k-1}(o,a)p_{0}(o|a_{t-k+1},o_{t-k+1})do,$
(B.1)

it follows that

$r_{k}(o_{t-k+1},a_{t-k+1})=\pi_{0}^{k}\widetilde{r}_{k}(o_{t-k+1},a_{t-k+1}).$

Using the same arguments, for a general $k$, we have

$\displaystyle
Q^{\pi}_{i,t-k+1,t}(O_{i,t-k+1},A_{i,t-k+1})=\sum_{a}{\mathbb{E}}[\pi(a|O_{i,t})Q^{\pi}_{i,t-k+2,t}(O_{i,t-k+2},a)|A_{i,t-k+1},O_{i,t-k+1}]$
$\displaystyle=$
$\displaystyle\theta_{k-1,i}+\pi_{u}\int_{o}\sum_{a}\pi(a|o)r_{k-1}(o,a)p_{u_{i}}(o)do+\lambda_{k-1,t}+\pi_{v}\int_{o}\sum_{a}\pi(a|o)r_{k-1}(o,a)p_{v_{t}}(o)do$
$\displaystyle+$
$\displaystyle\pi_{0}\int_{o}\sum_{a}\pi(a|o)r_{k-1}(o,a)p_{0}(o|a_{t-k+1},o_{t-k+1})do$
$\displaystyle=$
$\displaystyle\underbrace{\theta_{k-1,i}+\pi_{u}\pi_{0}^{k-1}\int_{o}\sum_{a}\pi(a|o)\widetilde{r}_{k-1}(o,a)p_{u_{i}}(o)do}_{\theta_{k,i}}+\underbrace{\lambda_{k-1,t}+\pi_{v}\pi_{0}^{k-1}\int_{o}\sum_{a}\pi(a|o)\widetilde{r}_{k-1}(o,a)p_{v_{t}}(o)do}_{\lambda_{k,t}}$
$\displaystyle+$
$\displaystyle\underbrace{\pi_{0}^{k}\int_{o}\sum_{a}\pi(a|o)\widetilde{r}_{k-1}(o,a)p_{0}(o|a_{t-k+1},o_{t-k+1})do}_{r_{k}(o_{t-k+1},a_{t-k+1})}.$
(B.2)

This completes the proof.

##  Appendix C Proof of Theorem 1

Proof: We prove Theorem 1 in this section. Recall that the proposed estimator
is Q-function based. To establish the convergence rate of the proposed value
estimator, it suffices to derive the convergence rate of the estimated
Q-function, or equivalently, the convergence rates of the estimated main
effect and two-way fixed effects.

Outline. Our proof is divided into several steps. We first show that the main
effect function is $p$-smooth under Assumption 1 in Section C.1. Together with
Assumption 2, this allows us to obtain upper bound of the model
misspecification error of the main effect function; see Equation (C.1) for
details. We next show that under the additivity assumption, the Q-function
converges exponentially fast to a constant as the number of iterations in the
backward induction algorithm (e.g., the index $k$ in Algorithm 1) approaches
infinity. In particular, we prove in Section C.2 that the main effect
converges to zero whereas the two fixed effects converge to two constants.
This ensures that the estimation error will not overly accumulate as $k$
approaches infinity. More specifically, the error bounds of the sieve
estimator and the fixed effects depend solely on $L$, $N$ and $T$, and are
independent of $k$. In Section C.3, we establish the lower bound of the
eigenvalue of the matrix $\bm{\Phi}^{T}(\bm{I}-\bm{P})\bm{\Phi}$, which allows
us to derive the error bound of the sieve estimator. Finally, we obtain the
rates of convergence of the sieve estimator and the estimated fixed effects in
Section C.4.

###  C.1 Hölder smoothness of the main effect function

We first introduce the class of $p$-smooth functions. For a given $d$-tuple
$\alpha=(\alpha_{1},\cdots,\alpha_{d})$ of non-negative integers, define
$D^{\alpha}$ as the following differential operator: for any given
$f:\mathcal{O}\to\mathbb{R}$,

$\displaystyle
D^{\alpha}f(o)=\frac{\partial^{\lVert\alpha\rVert_{1}}f(o)}{\partial
o_{1}^{\alpha_{1}}\dots\partial o_{d}^{\alpha_{d}}}.$

For a given $p>0$, let $\lfloor\cdot\rfloor$ denote the floor function such
that $\lfloor p\rfloor$ outputs the largest integer smaller than or equal to
$p$. The class of $p$-smooth functions on $\mathcal{O}$ is defined as follows:

$\displaystyle\Lambda(p,c)=\left\\{f:\sup_{\|\alpha\|_{1}\leq\lfloor
p\rfloor}\sup_{o\in\mathcal{O}}|D^{\alpha}f(o)|\leq
c,\sup_{\|\alpha\|_{1}=\lfloor
p\rfloor}\sup_{\begin{subarray}{c}o_{1},o_{2}\in\mathcal{O}\\\ o_{1}\neq
o_{2}\end{subarray}}\frac{|D^{\alpha}f(o_{1})-D^{\alpha}f(o_{2})|}{\|o_{1}-o_{2}\|_{2}^{p-\lfloor
p\rfloor}}\leq c\right\\},$

for some constant $c>0$.

We next summarize our results in Lemma 1 below.

######  Lemma 1.

Assume that there exists some constants $p$ and $C$, such that for any
$a\in\mathcal{A}$ and $o^{\prime}\in\mathcal{O}$, both $r_{1}(\cdot,a)$ and
$p_{0}(o^{\prime}|\cdot,a)$ belong to the class of $p$-smooth functions
$\Lambda(p,C)$. Then for any $t^{*}\geq 1$ and any $a$,
$r_{t^{*}}(\cdot,a)\in\Lambda(p,C^{\prime})$ for some $C^{\prime}>0$.

Proof: By the additivity assumption, for any $1\leq t\leq T$, we have
$Q^{\pi}_{i,t,t}(o,a)=\mathbb{E}(r_{1,i,t}|O_{i,t}=o,A_{i,t}=a)=r_{1}(o,a)+\theta_{i}+\lambda_{t}$.
By assumption, $r_{1}(\cdot,a)$ itself is a $p$-smooth function.

For any $t^{*}\geq 1$, notice that

$\displaystyle
Q^{\pi}_{i,t-t^{*},t}(o,a)=\sum_{a}{\mathbb{E}}[\pi(a^{\prime}|o^{\prime})Q^{\pi}_{i,t-t^{*}+1,t}(o^{\prime},a^{\prime})|a,o]$
$\displaystyle=$
$\displaystyle\int_{o^{\prime}}\sum_{a}\pi(a|o^{\prime})Q^{\pi}_{i,t-t^{*}+1,t}(o^{\prime},a^{\prime})p_{i,t}(o^{\prime}|o,a)do^{\prime}$
$\displaystyle=$
$\displaystyle\int_{o^{\prime}}\sum_{a}\pi(a|o^{\prime})Q^{\pi}_{i,t-t^{*}+1,t}(o^{\prime},a^{\prime})[\pi_{u}p_{u_{i}}(o^{\prime})+\pi_{v}p_{v_{t}}(o^{\prime})+\pi_{0}p_{0}(o^{\prime}|o,a)]do^{\prime}.$

Under the parameterization in Equation (B.2), we have

$\displaystyle
r_{t^{*}}(o,a)=\pi_{0}\int_{o^{\prime}}\sum_{a}\pi(a|o^{\prime})r_{t^{*}-1}(o^{\prime},a^{\prime})p_{0}(o^{\prime}|o,a)do^{\prime}.$

Suppose we have shown that $r_{t^{*}-1}(o,a)$ is $p$-smooth as a function of
$o$. To show $r_{t^{*}}(\cdot,a)$ is $p$-smooth, it is equivalent to show the
$p$-smoothness of
$\int_{o^{\prime}}\sum_{a}\pi(a|o^{\prime})r_{t^{*}-1}(o^{\prime},a^{\prime})p_{0}(o^{\prime}|o,a)do^{\prime}$.
This can be proven using Lemma 1 of Shi et al., (2022). Specifically, it can
be shown that, for any $d$-tuple $\alpha=(\alpha_{1},...,\alpha_{d})^{\top}$
of non-negative integers such that $\lVert\alpha\rVert_{1}\leq\lfloor
p\rfloor$,

$\left\lvert\sup_{\pi,o,a}D^{\alpha}\int_{o^{\prime}}\sum_{a}\pi(a|o^{\prime})r_{t^{*}-1}(o^{\prime},a^{\prime})p_{0}(o^{\prime}|o,a)do^{\prime}\right\rvert\leq
CR\lambda^{*}(\mathcal{O}),$

where $R=\sup_{t^{*},o,a}r_{t^{*}}(o,a)$ (which is bounded by induction) and
$\lambda^{*}(\mathcal{O})$ is the Lebesgue measure of the observation space
$\mathcal{O}$.

Moreover, it follows from the arguments in the proof of Lemma 1 of Shi et al.,
(2022) that for any $x\neq y$,

$\displaystyle\left\lvert
D^{\alpha}\int_{o^{\prime}}\sum_{a}\pi(a|o^{\prime})r_{t^{*}-1}(o^{\prime},a^{\prime})p_{0}(o^{\prime}|a,x)do^{\prime}-D^{\alpha}\int_{o^{\prime}}\sum_{a}\pi(a|o^{\prime})r_{t^{*}-1}(o^{\prime},a^{\prime})p_{0}(o^{\prime}|a,y)do^{\prime}\right\rvert$
$\displaystyle\leq CR\lambda^{*}(\mathcal{O})\lVert x-y\rVert_{2}^{p-\lfloor
p\rfloor},$

which implies that
$\int_{o^{\prime}}\sum_{a}\pi(a|o^{\prime})r_{t^{*}-1}(o^{\prime},a^{\prime})p_{0}(o^{\prime}|\cdot,a)do^{\prime}$
is in the class of functions $\Lambda(p,CR\lambda^{*}(\mathcal{O}))$. Hence,
$r_{t^{*}}(\cdot,a)$ is a $p$-smooth function, for any $1\leq t^{*}\leq t\leq
T$, $a\in\mathcal{A}$, $o\in\mathcal{O}$, and $\pi$. The proof is thus
completed. Although we focus on the parameterization in (B.2), it is easy to
show that Lemma 1 holds under other parameterizations as well, e.g., by adding
the constraints that $\sum_{i}\theta_{t^{*},i}=\sum_{t}\lambda_{t^{*},t}=0$.

Finally, combining Lemma 1 with Assumption 2(ii) yields that

$\displaystyle\max_{a,t^{*}}\inf_{\bm{\beta}\in\mathbb{R}^{L}}\sup_{o}|r_{t^{*}}(o,a)-\Phi_{L}^{\top}(o)\bm{\beta}|=O(L^{-p/d}).$
(C.1)

We use $\bm{\beta}_{t^{*},a}^{*}$ to denote the vector that the achieves the
infimum for any $a$ and $t^{*}$, and $\bm{\beta}_{t^{*}}^{*}$ to denote the
concatenation of $\\{\bm{\beta}_{t^{*},a}^{*}\\}_{a}$.

###  C.2 Q-function degeneracy

We elucidate the phenomenon of Q-function degeneracy in Lemma 2 below. As
commented in the outline, under the parameterization in (B.2), we aim to show
that the main effect function $r_{k}$ will converge exponentially fast to zero
with respect to $k$ whereas the two-way fixed effects will converge
exponentially fast to two constants. More generally, $r_{k}$ will converge
exponentially fast to a constant function under other parameterizations, e.g.,
by setting $\sum_{i}\theta_{k,i}=\sum_{t}\lambda_{k,t}=0$, based on the
observation that we can set $r_{k}\leftarrow
r_{k}+N^{-1}\sum_{i}\theta_{k,i}+T_{k}^{-1}\sum_{t}\lambda_{k,t}$,
$\theta_{k,i}\leftarrow\theta_{k,i}-N^{-1}\sum_{i}\theta_{k,i}$ and
$\lambda_{k,t}\leftarrow\lambda_{k,t}-T_{k}^{-1}\sum_{t}\lambda_{k,t}$ to
satisfy these constraints.

It implies that the Q-function will converge exponentially fast to a constant
function and that the Bellman error $\varepsilon_{k,i,t}$, for $k\geq 2$,

$\displaystyle\varepsilon_{k,i,t}\equiv\sum_{a}\pi(a|O_{i,t-k+2})Q^{\pi}_{i,t-k+2,t}(O_{i,t-k+2},a)-Q^{\pi}_{i,t-k+1,t}(O_{i,t-k+1},A_{i,t-k+1})=O(\pi_{0}^{k}),$
(C.2)

where the big-$O$ term is uniform in $i$ and $t$.

######  Lemma 2.

(i) $\sup_{o,a}|r_{k}(o,a)|=O(\pi_{0}^{k})$; (ii)
$\max_{i}|\theta_{k,i}-\theta_{k-1,i}|=O(\pi_{0}^{k})$ and
$\max_{t}|\lambda_{k,t}-\lambda_{k-1,t}|=O(\pi_{0}^{k})$; (iii) for $k\geq 2$,
the Bellman error $\varepsilon_{k,i,t}=O(\pi_{0}^{k})$.

Proof: By equation (B.2), we have for a general $k$,

$\displaystyle
Q^{\pi}_{i,t-k+1,t}(O_{i,t-k+1},A_{i,t-k+1})=\sum_{a}{\mathbb{E}}[\pi(a|O_{i,t})Q^{\pi}_{i,t-k+2,t}(O_{i,t-k+2},a)|A_{i,t-k+1},O_{i,t-k+1}]$
$\displaystyle=$
$\displaystyle\theta_{k-1,i}+\pi_{u}\int_{o}\sum_{a}\pi(a|o)r_{k-1}(o,a)p_{u_{i}}(o)do+\lambda_{k-1,t}+\pi_{v}\int_{o}\sum_{a}\pi(a|o)r_{k-1}(o,a)p_{v_{t}}(o)do$
$\displaystyle+$
$\displaystyle\pi_{0}\int_{o}\sum_{a}\pi(a|o)r_{k-1}(o,a)p_{0}(o|a_{t-k+1},o_{t-k+1})do$
$\displaystyle=$
$\displaystyle\underbrace{\theta_{k-1,i}+\pi_{u}\pi_{0}^{k-1}\int_{o}\sum_{a}\pi(a|o)\widetilde{r}_{k-1}(o,a)p_{u_{i}}(o)do}_{\theta_{k,i}}+\underbrace{\lambda_{k-1,t}+\pi_{v}\pi_{0}^{k-1}\int_{o}\sum_{a}\pi(a|o)\widetilde{r}_{k-1}(o,a)p_{v_{t}}(o)do}_{\lambda_{k,t}}$
$\displaystyle+$
$\displaystyle\underbrace{\pi_{0}^{k}\int_{o}\sum_{a}\pi(a|o)\widetilde{r}_{k-1}(o,a)p_{0}(o|a_{t-k+1},o_{t-k+1})do}_{r_{k}(o_{t-k+1},a_{t-k+1})}.$

We note that (i) and (ii) directly follow from equation (B.2) in Proposition
1. It remains to show (iii),

$\displaystyle\sum_{a}\pi(a|O_{i,t-k+2})Q^{\pi}_{i,t-k+2,t}(O_{i,t-k+2},a)=r_{k}(A_{i,t-k+1},O_{i,t-k+1})+\theta_{k,i}+\lambda_{k,t}+\varepsilon_{k,i,t},$

where

$\displaystyle\varepsilon_{k,i,t}=\sum_{a}\pi(a|O_{i,t-k+2})Q^{\pi}_{i,t-k+2,t}(O_{i,t-k+2},a)-r_{k}(A_{i,t-k+1},O_{i,t-k+1})-\theta_{k,i}-\lambda_{k,t}$
$\displaystyle=$
$\displaystyle\;\theta_{k-1,i}+\lambda_{k-1,t}+\sum_{a}\pi(a|O_{i,t-k+2})r_{k-1}(a,O_{i,t-k+2})-r_{k}(A_{i,t-k+1},O_{i,t-k+1})-\theta_{k,i}-\lambda_{k,t}$
$\displaystyle=$
$\displaystyle\;\underbrace{\theta_{k-1,i}-\theta_{k,i}}_{D_{1}}+\underbrace{\lambda_{k-1,t}-\lambda_{k,t}}_{D_{2}}+\underbrace{\sum_{a}\pi(a|O_{i,t-k+2})r_{k-1}(a,O_{i,t-k+2})-r_{k}(A_{i,t-k+1},O_{i,t-k+1})}_{D_{3}}.$

By equation (B.2), the first term is

$D_{1}=\pi_{u}\pi_{0}^{k-1}\int_{o}\sum_{a}\pi(a|o)\widetilde{r}_{k-1}(o,a)p_{u_{i}}(o)do,$

and the second term

$D_{2}=\pi_{v}\pi_{0}^{k-1}\int_{o}\sum_{a}\pi(a|o)\widetilde{r}_{k-1}(o,a)p_{v_{t}}(o)do.$

The third term $D_{3}$ can be written as

$\displaystyle\sum_{a}\pi(a|O_{i,t-k+2})r_{k-1}(a,O_{i,t-k+2})-r_{k}(A_{i,t-k+1},O_{i,t-k+1})$
$\displaystyle=\;$
$\displaystyle\sum_{a}\pi(a|O_{i,t-k+2})r_{k-1}(a,O_{i,t-k+2})-\pi_{0}^{k}\int_{o}\sum_{a}\pi(a|o)\widetilde{r}_{k-1}(o,a)p_{0}(o|A_{t-k+1},O_{t-k+1})do$
$\displaystyle=\;$
$\displaystyle\pi_{0}^{k-1}\sum_{a}\pi(a|O_{i,t-k+2})\widetilde{r}_{k-1}(a,O_{i,t-k+2})-\pi_{0}^{k}\int_{o}\sum_{a}\pi(a|o)\widetilde{r}_{k-1}(o,a)p_{0}(o|A_{t-k+1},O_{t-k+1})do.$

By Equation (B.1), it can be shown that for any $k$ and $a\in\mathcal{A}$,
$\widetilde{r}_{k}(o,a)$ is strictly bounded, i.e.,
$\sup_{o,a}|\widetilde{r}_{k}(o,a)|\leq c$. Thus we have

$|D_{1}|\leq\pi_{u}\pi_{0}^{k-1}c_{1},\mbox{
}|D_{2}|\leq\pi_{v}\pi_{0}^{k-1}c_{2},\mbox{ and }$
$\displaystyle|D_{3}|\leq\pi_{0}^{k-1}\left\lvert\sum_{a}\pi(a|O_{i,t-k+2})\widetilde{r}_{k-1}(a,O_{i,t-k+2})\right\rvert+\pi_{0}^{k}\left\lvert\int_{o}\sum_{a}\pi(a|o)\widetilde{r}_{k-1}(o,a)p_{0}(o|A_{t-k+1},O_{t-k+1})do\right\rvert$
$\displaystyle\leq$ $\displaystyle\;\pi_{0}^{k-1}(1+\pi_{0})c.$

As such, there exists a constant $C$ such that

$\displaystyle|\varepsilon_{k,i,t}|\leq|D_{1}|+|D_{2}|+|D_{3}|\leq$
$\displaystyle\;C\pi_{0}^{k},$ (C.3)

which suggests that $\varepsilon_{k,i,t}$ converges to $0$ exponentially fast
with respect to $k$. This completes the proof.

###  C.3 Eigenvalue Results

######  Lemma 3.

The minimum eigenvalue of $\bm{\Phi}^{\top}(\bm{I}-\bm{P})\bm{\Phi}/n$ is
lower bounded by $\rho_{0}/2$, with probability at least $1-n^{-c}$. Moreover,
$n^{-1}\|\bm{\Phi}\bm{\Phi}^{\top}-{\mathbb{E}}\bm{\Phi}\bm{\Phi}^{\top}\|_{2}=O_{p}(\log
n\sqrt{L/n})$.

Proof: Under Assumptions 2(iii), 3 and 4, the above results hold by directly
applying the matrix concentration inequality under nonstationary Markov
decision process developed in Lemma 2 of Li et al., (2022). We omit the
details to save the space.

###  C.4 Rates of convergence

Finally, we obtain the convergence rates of the linear sieve estimator and the
two-way fixed effects estimators in this section. We start by analyzing the
properties of $\\{\widehat{\theta}_{i}\\}_{i}$,
$\\{\widehat{\lambda}_{t}\\}_{t}$ and $\widehat{r}$ in Sections C.4.1 and
C.4.2, which we refer to as the first iteration analysis. We next proceed with
the second iteration analysis, i.e., to establish the convergence rates for
$\\{\widehat{\theta}_{2,i}\\}_{i}$, $\\{\widehat{\lambda}_{2,t}\\}_{t}$ and
$\widehat{r}_{2}$ in Sections C.4.3 and C.4.4. Finally, we extend these
results to obtain the uniform convergence rates of
$\\{\widehat{\theta}_{k,i}\\}_{k,i}$, $\\{\widehat{\lambda}_{k,t}\\}_{k,t}$
and $\\{\widehat{r}_{k}\\}_{k}$ in Sections C.4.5 and C.4.6.

####  C.4.1 First iteration analysis: linear sieve estimator

Part of this subsection’s proof utilizes the theoretical framework developed
in Shi et al., (2022) and Li et al., (2022) to establish the property of the
linear sieves estimator. Recall that
$\bm{P}=\bm{B}(\bm{B}^{\top}\bm{B})^{+}\bm{B}^{\top}$, hence

$\displaystyle\widehat{\bm{\beta}}=$
$\displaystyle(\bm{\Phi}^{\top}\underbrace{(\bm{I}-\bm{P})}_{\bm{S}}\bm{\Phi})^{-1}\bm{\Phi}^{\top}(\bm{I}-\bm{P})\bm{R}$
$\displaystyle=$
$\displaystyle(\bm{\Phi}^{\top}\bm{S}\bm{\Phi})^{-1}\bm{\Phi}^{\top}\bm{S}(\bm{B}\bm{\alpha}+\bm{M}+\bm{\varepsilon})$
$\displaystyle=$
$\displaystyle(\bm{\Phi}^{\top}\bm{S}\bm{\Phi})^{-1}\bm{\Phi}^{\top}\bm{S}(\bm{M}+\bm{\varepsilon})$
$\displaystyle=$
$\displaystyle(\bm{\Phi}^{\top}\bm{S}\bm{\Phi})^{-1}\bm{\Phi}^{\top}\bm{S}(\bm{\Phi}\bm{\beta}^{*}+\bm{M}-\bm{\Phi}\bm{\beta}^{*}+\bm{\varepsilon}),$

where the third equation using the fact that $\bm{S}\bm{B}=\bm{0}$, and the
error term

$\bm{\varepsilon}=(\varepsilon_{1,1},\varepsilon_{2,1},\dots,\varepsilon_{N,1},\varepsilon_{2,1},\dots,\varepsilon_{N-1,T},\varepsilon_{N,T})^{\top}\in\mathbb{R}^{NT}.$

Let $\\{\bm{\beta}^{*}_{a}\\}_{a}$ denote the set of vectors such that
$\sup_{o,a}|r_{1}(o,a)-\bm{\Phi}_{L}(o)^{\top}\bm{\beta}^{*}_{a}|=O(L^{-p/d})$
and $\bm{\beta}^{*}$ the concatenation of $\\{\bm{\beta}^{*}_{a}\\}_{a}$.
These vectors exist according to Equation (C.1). It follows that

$\displaystyle\widehat{\bm{\beta}}-\bm{\beta}^{*}=$
$\displaystyle(\bm{\Phi}^{\top}\bm{S}\bm{\Phi})^{-1}\bm{\Phi}^{\top}\bm{S}(\bm{M}-\bm{\Phi}\bm{\beta}^{*})+(\bm{\Phi}^{\top}\bm{S}\bm{\Phi})^{-1}\bm{\Phi}^{\top}\bm{S}\bm{\varepsilon}+\underbrace{(\bm{\Phi}^{\top}\bm{S}\bm{\Phi})^{-1}\bm{\Phi}^{\top}\bm{S}\bm{\Phi}\bm{\beta}^{*}-\bm{\beta}^{*}}_{0}$
$\displaystyle=$
$\displaystyle(\underbrace{\bm{\Phi}^{\top}\bm{S}\bm{\Phi}}_{\widehat{\bm{\Sigma}}})^{-1}\bm{\Phi}^{\top}\bm{S}(\bm{M}-\bm{\Phi}\bm{\beta}^{*})+(\bm{\Phi}^{\top}\bm{S}\bm{\Phi})^{-1}\bm{\Phi}^{\top}\bm{S}\bm{\varepsilon}$
$\displaystyle=$
$\displaystyle\underbrace{\widehat{\bm{\Sigma}}^{-1}\bm{\Phi}^{\top}\bm{S}(\bm{M}-\bm{\Phi}\bm{\beta}^{*})}_{I_{1}}+\underbrace{\widehat{\bm{\Sigma}}^{-1}\bm{\Phi}^{\top}\bm{S}\bm{\varepsilon}}_{I_{2}}.$
(C.4)

In the following, we will show that $\lVert I_{1}\rVert_{2}=O_{p}(L^{-p/d})$
and $\lVert I_{2}\rVert_{2}=O_{p}(\sqrt{L/n})$.

Error bound for $\lVert I_{1}\rVert_{2}$:

$\displaystyle\lVert\bm{\Phi}^{\top}\bm{S}(\bm{M}-\bm{\Phi}\bm{\beta}^{*})\rVert^{2}_{2}=(\bm{M}-\bm{\Phi}\bm{\beta}^{*})^{\top}\bm{S}^{\top}\bm{\Phi}\bm{\Phi}^{\top}\bm{S}(\bm{M}-\bm{\Phi}\bm{\beta}^{*})$
$\displaystyle\leq$
$\displaystyle\>\lambda_{\max}(\bm{\Phi}\bm{\Phi}^{\top})(\bm{M}-\bm{\Phi}\bm{\beta}^{*})^{\top}\bm{S}^{\top}\bm{S}(\bm{M}-\bm{\Phi}\bm{\beta}^{*})$
$\displaystyle\leq$
$\displaystyle\>\lambda_{\max}(\bm{\Phi}\bm{\Phi}^{\top})\;\lVert\bm{M}-\bm{\Phi}\bm{\beta}^{*}\rVert^{2}_{2}.$

Under Assumptions 2(iii), 3 and 4, by Lemma 3, the minimum eigenvalue of
$\widehat{\bm{\Sigma}}$ is lower bounded by $n\rho_{0}/2$, with probability at
least $1-n^{-c}$. In addition, the difference
$n^{-1}\|\bm{\Phi}\bm{\Phi}^{\top}-{\mathbb{E}}\bm{\Phi}\bm{\Phi}^{\top}\|_{2}$
decays to zero with probability at least $1-n^{-c_{1}}$. Under Assumption 1,
$p_{0}$ is bounded. This together with Assumption 3(iii) implies the
boundedness of the transition function $p$ in (2.3) see e.g., Equation (60) of
Shi et al., (2022). Together with Assumption 2(i), we obtain the boundedness
of $n^{-1}{\mathbb{E}}\|\bm{\Phi}\bm{\Phi}^{\top}\|_{2}$. As such,
$\lambda_{\max}(\bm{\Phi}\bm{\Phi}^{\top})=O(n)$ with probability at least
$1-n^{-c}$. Finally, by the definition of $\bm{\beta}^{*}$, we have
$\lVert\bm{M}-\bm{\Phi}\bm{\beta}^{*}\rVert^{2}_{2}=O(nL^{-2p/d})$. This leads
to the observation

$\displaystyle\lVert\bm{\Phi}^{\top}\bm{S}(\bm{M}-\bm{\Phi}\bm{\beta}^{*})\rVert^{2}_{2}\leq\>Cn^{2}L^{-2p/d},$

with probability at least $1-O(n^{-c})$. Thus, we obtain the following error
bound for $\lVert I_{1}\rVert_{2}$,

$\displaystyle\lVert
I_{1}\rVert_{2}\leq\lVert\widehat{\bm{\Sigma}}^{-1}\rVert_{2}\lVert\bm{\Phi}^{\top}\bm{S}(\bm{M}-\bm{\Phi}\bm{\beta}^{*})\rVert_{2}\leq\>2\rho_{0}^{-1}CL^{-p/d}.$
(C.5)

Error bound for $\lVert I_{2}\rVert_{2}$: Recall that
$\bm{\Phi}^{\top}\bm{S}\bm{\varepsilon}=\bm{\Phi}^{\top}\bm{\varepsilon}-\bm{\Phi}^{\top}\bm{P}\bm{\varepsilon}$.
We first consider $\bm{\Phi}^{\top}\bm{\varepsilon}$. Notice that the
observations across different trajectories are independent. In addition,
conditional on the two-way fixed effects, the Markov assumption is satisfied
and we have

$\displaystyle\mathbb{E}(\varepsilon_{i,t}\varepsilon_{i,j}|A_{i,t},O_{i,t})=0\mbox{
and
}\mathbb{E}(\varepsilon_{i_{1},t_{1}}\varepsilon_{i_{2},t_{2}}|A_{i_{1},t_{1}},O_{i_{1},t_{1}},A_{i_{2},t_{2}},O_{i_{2},t_{2}})=0,$
(C.6)

for any $i,t_{1},t_{2}$, $i_{1}\neq i_{2}$ and $t<j$. It follows that

$\displaystyle\mathbb{E}\lVert\bm{\Phi}^{\top}\bm{\varepsilon}\rVert^{2}_{2}=\mathbb{E}(\bm{\varepsilon}^{\top}\bm{\Phi}\bm{\Phi}^{\top}\bm{\varepsilon})=\mathbb{E}\operatorname*{trace}(\bm{\varepsilon}^{\top}\bm{\Phi}\bm{\Phi}^{\top}\bm{\varepsilon})=\mathbb{E}\operatorname*{trace}(\bm{\Phi}^{\top}\bm{\varepsilon}\bm{\varepsilon}^{\top}\bm{\Phi})$
$\displaystyle\leq\;$ $\displaystyle
C\mathbb{E}\operatorname*{trace}(\bm{\Phi}^{\top}\bm{\Phi})=O(L\lambda_{\max}({\mathbb{E}}\bm{\Phi}^{\top}\bm{\Phi}))=\;O(nL),$

where the last equality follows Assumption 2. Using the same arguments, we can
show that
${\mathbb{E}}\lVert\bm{\Phi}^{\top}\bm{S}\bm{\varepsilon}\rVert^{2}_{2}=O(nL)$.
Combined with $\lVert\widehat{\bm{\Sigma}}^{-1}\rVert_{2}=O_{p}(n^{-1})$, this
yields $\lVert I_{2}\rVert_{2}=O_{p}(\sqrt{L/n})$.

To summarize, we have shown that

$\displaystyle\lVert\widehat{\bm{\beta}}-\bm{\beta}^{*}\rVert_{2}=O_{p}(L^{-p/d})+O_{p}(\sqrt{L/n}).$
(C.7)

####  C.4.2 First iteration analysis: fixed effects estimator

Given the constraint $\sum_{i=1}^{N}\theta_{i}=\sum_{t=1}^{T}\lambda_{t}=0$,
we obtain

$\widehat{\lambda}_{t}=N^{-1}\sum_{i=1}^{N}(R_{i,t}-\widehat{r}_{1}(o_{i,t},a_{i,t})),$

and

$\widehat{\theta}_{i}=T^{-1}\sum_{t=1}^{T}(R_{i,t}-\widehat{r}_{1}(o_{i,t},a_{i,t})).$

We first focus on $\widehat{\lambda}_{t}$ and observe that

$\displaystyle\widehat{\lambda}_{t}-\lambda_{t}=N^{-1}\sum_{i=1}^{N}(R_{i,t}-\widehat{r}_{1}(o_{i,t},a_{i,t}))-\lambda_{t}$
$\displaystyle=$
$\displaystyle\;N^{-1}\sum_{i=1}^{N}(R_{i,t}+r_{1,i,t}-\widehat{r}_{1,i,t}-r_{1,i,t})-\lambda_{t}$
$\displaystyle=$
$\displaystyle\;\underbrace{N^{-1}\sum_{i=1}^{N}(R_{i,t}-r_{1,i,t})-\lambda_{t}}_{F_{1,t}}+\underbrace{N^{-1}\sum_{i=1}^{N}(r_{1,i,t}-\widehat{r}_{1,i,t})}_{F_{2,t}}.$
(C.8)

The first term $F_{1,t}$ is of the order $O_{p}(N^{-1/2})$, as it is an
average of conditionally independent mean zero random variables.

As for $F_{2,t}$, we have

$\displaystyle
F_{2,t}=\frac{1}{N}\sum_{a}\sum_{i=1}^{N}\mathbb{I}(A_{i,t}=a)\bm{\Phi}_{L}^{\top}(O_{i,t})(\widehat{\bm{\beta}}_{a}-\bm{\beta}_{a})+O_{p}(L^{-p/d}).$

A naive method is to directly apply Cauchy-Schwarz inequality to the right-
hand-side and show that it is of the order
$O(\|\widehat{\bm{\beta}}-\bm{\beta}\|_{2})$. However, this would yield a sub-
optimal bound. To motivate our analysis, we notice that $F_{2,t}$ can be upper
bounded as

$\displaystyle\frac{1}{N}\lVert\sum_{i=1}^{N}\bm{\Phi}_{L}(O_{i,t})\mathbb{I}(A_{i,t}=a)-\sum_{i=1}^{N}\mathbb{E}\bm{\Phi}_{L}(O_{i,t})\mathbb{I}(A_{i,t}=a)\rVert_{2}\lVert\widehat{\bm{\beta}}-\bm{\beta}\rVert_{2}$
$\displaystyle+\frac{1}{N}\sum_{a}|\sum_{i=1}^{N}[\mathbb{E}\mathbb{I}(A_{i,t}=a)\bm{\Phi}_{L}^{\top}(O_{i,t})](\widehat{\bm{\beta}}_{a}-\bm{\beta}_{a})|.$

Similar to Lemma 3 of Shi et al., (2022), we can show that
$N^{-1}\lVert\sum_{i=1}^{N}\bm{\Phi}_{L}(O_{i,t})-\sum_{i=1}^{N}\mathbb{E}\bm{\Phi}_{L}(O_{i,t})\rVert_{2}$
is of the order $\sqrt{N^{-1}L\log(NT)}$, with probability at least
$1-n^{-c}$. As such, under Assumption 2(iii), the first line is of the order
$\sqrt{N^{-1}L\log(NT)}(L^{-p/d}+\sqrt{L/NT})=o(L^{-p/d})+o(N^{-1/2})$. As for
the second line, notice that
$\phi_{a,t}=N^{-1}\sum_{i=1}^{N}[\mathbb{E}\mathbb{I}(A_{i,t}=a)\bm{\Phi}_{L}^{\top}(O_{i,t})]$
is a fixed vector with bounded $\ell_{2}$ norm under Assumption 2(ii).
Specifically, according to Cauchy-Schwarz inquality, we have

$\displaystyle\|\phi_{a,t}\|_{2}=\sup_{\bm{\alpha}:\|\bm{\alpha}\|_{2}=1}\left|\frac{1}{N}\sum_{i=1}^{N}\mathbb{E}\mathbb{I}(A_{i,t}=a)\bm{\Phi}_{L}^{\top}(O_{i,t})\bm{\alpha}\right|$
$\displaystyle\leq$
$\displaystyle\sup_{\bm{\alpha}:\|\bm{\alpha}\|_{2}=1}\sqrt{\frac{1}{N}\sum_{i=1}^{N}\mathbb{E}\mathbb{I}(A_{i,t}=a)[\bm{\Phi}_{L}^{\top}(O_{i,t})\bm{\alpha}]^{2}}\leq\lambda_{\max}\left[\frac{1}{N}\sum_{i=1}^{N}{\mathbb{E}}\bm{\Phi}_{L}(O_{i,t})\bm{\Phi}_{L}^{\top}(O_{i,t})\right],$

which is $O(1)$ using similar arguments in bounding $\|I_{1}\|_{2}$ in Section
C.4.1. The bias $|\phi_{a,t}^{\top}(\widehat{\bm{\beta}}_{a}-\bm{\beta}_{a})|$
is of the same order as $\|\widehat{\bm{\beta}}_{a}-\bm{\beta}_{a}\|_{2}$,
given by $O(L^{-p/d})$. The variance term, however, can be substantially
reduced from $O(n^{-1}L)$ to $O(n^{-1})$. In particular, using similar
arguments in bounding $\|I_{2}\|_{2}$ in Section C.4.1, we can show that

$\displaystyle{\mathbb{E}}|\phi_{t}^{\top}\bm{\Phi}^{\top}\bm{\varepsilon}|^{2}=O(\phi_{t}^{\top}{\mathbb{E}}\bm{\Phi}^{\top}\bm{\Phi}\phi_{t})=O(\lambda_{\max}[{\mathbb{E}}\bm{\Phi}^{\top}\bm{\Phi}])=O(n),$

where $\phi_{t}$ is obtained by concatenating the set of vectors
$\\{\phi_{a,t}\\}_{a}$. This allows us to show that the variance term is of
the order $O(n^{-1})$. As such, the second term is of the order
$O(L^{-p/d})+O_{p}(n^{-1/2})$. Hence, $F_{2,t}=O(L^{-p/d})+O_{p}(N^{-1/2})$.
Therefore,

$\displaystyle\widehat{\lambda}_{t}-\lambda_{t}=O(L^{-p/d})+O_{p}(N^{-1/2}).$
(C.9)

Notice that (C.9) only establishes the convergence rate at a single time point
$t$. Conditional on the fixed effects, the Markov assumption is satisfied and
the error term $\\{\varepsilon_{i,t}\\}_{t}$ forms a martingale difference
sequence (see e.g., Shi et al.,, 2022, Step 3 of the proof of Theorem 1).
Using the martingale concentration inequality (Tropp,, 2011), we can show with
probability at least $1-O(n^{-c})$ that
$\max_{t}|F_{1,t}|=O(N^{-1/2}\sqrt{\log(NT)})$ and
$\max_{t}|F_{2,t}|=O(L^{-p/d})+O(N^{-1/2})$. This yields the uniform
convergence rate

$\displaystyle\max_{t}|\widehat{\lambda}_{t}-\lambda_{t}|=O(L^{-p/d})+O(N^{-1/2}\sqrt{\log
n}),$ (C.10)

with probability at least $1-O(n^{-c})$.

Similarly, we can show with probability at least $1-O(n^{-c})$ that

$\displaystyle\max_{i}|\widehat{\theta}_{i}-\theta_{i}|=O(L^{-p/d})+O(T^{-1/2}\sqrt{\log
n}).$ (C.11)

####  C.4.3 Second iteration analysis: linear sieve estimator

In the second iteration (corresponding to $k=2$ in Algorithm 1), we regress
$\sum_{a}\pi(a|O_{i,t})\widehat{Q}_{i,t,t}(O_{i,t},a)$ on
$(O_{i,t-1},A_{i,t-1})$. By Proposition 1,
$Q^{\pi}_{i,t-1,t}(O_{i,t-1},A_{i,t-1})$ also has a two-way structure. It
follows from (C.2) that

$\displaystyle\sum_{a}\pi(a|O_{i,t})Q^{\pi}_{i,t,t}(O_{i,t},a)=r_{2}(O_{i,t-1},A_{i,t-1})+\theta_{2,i}+\lambda_{2,t}+\varepsilon_{2,i,t}.$
(C.12)

However, the oracle Q-function is unknown, and we only have access to its
estimator $\widehat{Q}^{\pi}_{i,t,t}$. As such, we have

$\displaystyle\sum_{a}\pi(a|O_{i,t})\widehat{Q}^{\pi}_{i,t,t}(O_{i,t},a)$
$\displaystyle=$ $\displaystyle
r_{2}(A_{i,t-1},O_{i,t-1})+\theta_{2,i}+\lambda_{2,t}+\underbrace{\varepsilon_{2,i,t}+\sum_{a}\pi(a|O_{i,t})\widehat{Q}^{\pi}_{i,t,t}(O_{i,t},a)-\sum_{a}\pi(a|O_{i,t})Q^{\pi}_{i,t,t}(O_{i,t},a)}_{\widetilde{\varepsilon}_{2,i,t},\mbox{
the adjusted error term}}$ $\displaystyle=$ $\displaystyle
r_{2}(A_{i,t-1},O_{i,t-1})+\theta_{2,i}+\lambda_{2,t}+\widetilde{\varepsilon}_{2,i,t}.$

Similarly to Equation (C.4), the linear sieve estimator
$\widehat{\bm{\beta}}_{2}$ at the second iteration satisfies

$\displaystyle\widehat{\bm{\beta}}_{2}-\bm{\beta}^{*}_{2}=\widehat{\bm{\Sigma}}_{2}^{-1}\bm{\Phi}_{2}^{\top}\bm{S}_{2}(\bm{M}_{2}-\bm{\Phi}_{2}\bm{\beta}_{2}^{*})+\widehat{\bm{\Sigma}}_{2}^{-1}\bm{\Phi}_{2}^{\top}\bm{S}_{2}\widetilde{\bm{\varepsilon}}_{2},$

where $\widehat{\bm{\Sigma}}_{2}$, $\bm{\Phi}_{2}$, $\bm{S}_{2}$ is similarly
defined as in first iteration analysis. The first term (e.g., the bias term)
is still of the order $O(L^{-p/d})$, with probability at least $1-O(n^{-c})$.
We now focus on the second term

$\displaystyle\widehat{\bm{\Sigma}}_{2}^{-1}\bm{\Phi}_{2}^{\top}\bm{S}_{2}\widetilde{\bm{\varepsilon}_{2}}=\widehat{\bm{\Sigma}}_{2}^{-1}\bm{\Phi}_{2}^{\top}\bm{S}_{2}(\bm{\varepsilon}_{2}+\bm{\varepsilon}_{1}^{\prime}+\bm{\varepsilon}_{1}^{\prime\prime}+\bm{\varepsilon}_{1}^{\prime\prime\prime})$
$\displaystyle=$
$\displaystyle\underbrace{\;\widehat{\bm{\Sigma}}_{2}^{-1}\bm{\Phi}_{2}^{\top}\bm{S}_{2}\bm{\varepsilon}_{2}}_{I_{2,1}}+\underbrace{\widehat{\bm{\Sigma}}_{2}^{-1}\bm{\Phi}_{2}^{\top}\bm{S}_{2}\bm{\varepsilon}_{1}^{\prime}}_{I_{2,2}}+\underbrace{\widehat{\bm{\Sigma}}_{2}^{-1}\bm{\Phi}_{2}^{\top}\bm{S}_{2}\bm{\varepsilon}_{1}^{\prime\prime}}_{I_{2,3}}+\underbrace{\widehat{\bm{\Sigma}}_{2}^{-1}\bm{\Phi}_{2}^{\top}\bm{S}_{2}\bm{\varepsilon}_{1}^{\prime\prime\prime}}_{I_{2,4}},$

where
$\bm{\varepsilon}_{1}^{\prime},\bm{\varepsilon}_{1}^{\prime\prime},\bm{\varepsilon}_{1}^{\prime\prime\prime}$
represent the errors resulting from the estimation of individual effects, time
effects, and the linear sieve in the first iteration respectively, i.e.,

$\displaystyle\varepsilon_{1,i,t}^{\prime}=\frac{1}{T}\sum_{t=1}^{T}\varepsilon_{1,i,t},\,\,\,\,\varepsilon_{1,i,t}^{\prime\prime}=\frac{1}{N}\sum_{i=1}^{N}\varepsilon_{1,i,t},\,\,\,\,\bm{\varepsilon}_{1}^{\prime\prime\prime}=\bm{\Phi}_{1}^{new}\widehat{\bm{\Sigma}}_{1}^{-1}\bm{\Phi}_{1}^{\top}\bm{S}_{1}(\bm{M}_{1}-\bm{\Phi}_{1}\bm{\beta}_{1}^{*}+\bm{\varepsilon}_{1}),$

where $\bm{\Phi}_{1}^{new}$ is defined in (5.1). Using similar arguments in
bounding $I_{1}$ in the first iteration analysis, we can show that $\lVert
I_{2,1}\rVert$ is of the order $O_{p}(\sqrt{L/n_{2}})$. We next analyse
$I_{2,2}=\widehat{\bm{\Sigma}}_{2}^{-1}\bm{\Phi}_{2}^{\top}\bm{S}_{2}\bm{\varepsilon}_{1}^{\prime}$.
Similarly to the first iteration analysis, we first bound
$\lVert\bm{\Phi}_{2}^{\top}\bm{\varepsilon}_{1}^{\prime}\rVert_{2}$, as
$\lVert\bm{\Phi}_{2}^{\top}\bm{P}_{2}\bm{\varepsilon}_{1}^{\prime}\rVert_{2}$
can be obtained similarly. Note that the term $\bm{\varepsilon}_{1}^{\prime}$
has a special structure, it is constructed such that the first $N$ components
of the vector are replicated $T_{2}$ times, i.e.,

$\bm{\varepsilon}_{1}^{\prime}=(\underbrace{\bm{\varepsilon}_{1N}^{{}^{\prime}\top},\bm{\varepsilon}_{1N}^{{}^{\prime}\top},\dots,\bm{\varepsilon}_{1N}^{{}^{\prime}\top}}_{T_{2}\mbox{
times}})^{\top},$

where we use $\bm{\varepsilon}_{1N}^{\prime}$ to denote the first $N$ elements
in $\bm{\varepsilon}_{2}^{\prime}$. Similarly to the first step analysis, it
suffices to bound $\bm{\Phi}_{2}^{\top}\bm{\varepsilon}_{1}^{\prime}$, which
can be written as
$\sum_{t=1}^{T_{2}}\bm{\Phi}_{2,t}^{\top}\bm{\varepsilon}_{1N}^{\prime},$
where $\bm{\Phi}_{2,t}^{\top}$ is the $t$th submatrix of
$\bm{\Phi}_{2}^{\top}$ (i.e., write $\bm{\Phi}_{2}^{\top}$ as
$[\bm{\Phi}_{2,1}^{\top},\bm{\Phi}_{2,2}^{\top}\dots,\bm{\Phi}_{2,T_{2}}^{\top}]$).
Then we have

$\displaystyle\lVert\bm{\Phi}_{2}^{\top}\bm{\varepsilon}_{2}^{\prime}\rVert_{2}=\lVert\sum_{t=1}^{T_{2}}\bm{\Phi}_{2,t}^{\top}\bm{\varepsilon}_{1N}^{\prime}\rVert_{2}$
$\displaystyle\leq$
$\displaystyle\;\sqrt{T_{2}}\sqrt{\sum_{t=1}^{T_{2}}\operatorname*{trace}(\bm{\varepsilon}_{1N}^{{}^{\prime}\top}\bm{\Phi}_{2,t}^{\top}\bm{\Phi}_{2,t}\bm{\varepsilon}_{1N}^{\prime})}$
$\displaystyle\leq$
$\displaystyle\;C\sqrt{\sum_{t=1}^{T_{2}}\operatorname*{trace}(\bm{\Phi}_{2,t}^{\top}\bm{\Phi}_{2,t})}=C\sqrt{\operatorname*{trace}(\bm{\Phi}_{2}^{\top}\bm{\Phi}_{2})}=C\sqrt{n_{2}L},$

where we apply the AM-QM inequality in the second inequality. Thus, we have
showed that $\lVert I_{2,2}\rVert=O_{p}(\sqrt{L/n_{2}})$. Applying the same
technique, $\lVert I_{2,3}\rVert$ can be also bounded by
$O_{p}(\sqrt{L/n_{2}})$. Next we provide the error bound for $\lVert
I_{2,4}\rVert$. It follows that

$\displaystyle\lVert
I_{2,4}\rVert_{2}=\lVert\widehat{\bm{\Sigma}}_{2}^{-1}\bm{\Phi}_{2}^{\top}\bm{S}_{2}\bm{\varepsilon}_{1}^{\prime\prime\prime}\rVert_{2}=\lVert\widehat{\bm{\Sigma}}_{2}^{-1}\bm{\Phi}_{2}^{\top}\bm{S}_{2}\bm{\Phi}_{1}^{new}\widehat{\bm{\Sigma}}_{1}^{-1}\bm{\Phi}_{1}^{\top}\bm{S}_{1}\bm{\varepsilon}_{1}\rVert_{2}$
$\displaystyle\leq$
$\displaystyle\;\underbrace{\lVert\widehat{\bm{\Sigma}}_{2}^{-1}\bm{\Phi}_{2}^{\top}\bm{S}_{2}\bm{\Phi}_{1}^{new}\rVert_{2}}_{O_{p}(1)\mbox{
by Assumption
4}}\lVert\widehat{\bm{\Sigma}}_{1}^{-1}\bm{\Phi}_{1}^{\top}\bm{S}_{1}\bm{\varepsilon}_{1}\rVert_{2}=O_{p}(\sqrt{L/n_{2}}),$
(C.13)

by the result in first step analysis. In summary, we have shown that the
linear sieve estimator $\widehat{\bm{\beta}}_{2}$ achieves the same the
estimation error bound, given by $O_{p}(L^{-p/d})+O_{p}(\sqrt{L/n_{2}})$.

####  C.4.4 Second iteration analysis: fixed effects estimator

The second iteration analysis for the fixed effects is very similar to the
first iteration analysis. The estimated fixed effects achieve the same uniform
error bounds as in Equations (C.10) and (C.11). We omit the details to save
space.

####  C.4.5 Kth iteration analysis: linear sieve Estimator

Notation: We use the following notation in the $k$th iteration analysis.
$\bm{S}_{k}$ is the residual maker matrix at the $k$th step, i.e.,
$\bm{S}_{k}=\bm{I}-\bm{B}_{k}(\bm{B}_{k}^{\top}\bm{B}_{k})^{+}\bm{B}_{k}^{\top}$,
$\bm{B}_{k}$ is the design matrix for the fixed effects at step $k$, i.e.,
$\bm{B}_{k}=(\bm{1}_{T-k+1}\otimes\bm{I}_{N},\bm{I}_{T-k+1}\otimes\bm{1}_{N})$,
and $\bm{\Phi}_{k}^{new}$ is the design matrix such that

$\bm{\Phi}^{new}_{k,i,t}=\frac{1}{T-k+1}\sum_{t=1}^{T-k+1}\bm{\Phi}_{k,i,t}+\frac{1}{N}\sum_{i=1}^{N}\bm{\Phi}_{k,i,t}+\sum_{a\in\mathcal{A}}\pi(a|O_{k,i,t})\bm{\Phi}_{k,i,t},$

and $\bm{\Phi}_{k}$ is defined as the matrix consisting of the first $n_{k}$
rows of matrix $\bm{\Phi}$.

We consider the scenario that the backward induction has been conducted $K-1$
times and allow $K$ to approach infinity. It follows that the estimation error
at $k$th step $E_{k,i,t}$ is

$\displaystyle
E_{k,i,t}=\widehat{\theta}_{k,i}-\theta_{k,i}+\widehat{\lambda}_{k,t}-\lambda_{k,t}+\sum_{a}\pi(a|o)(\widehat{r}_{k,i,t}-r_{k,i,t})$
$\displaystyle=$
$\displaystyle\;W_{k,i,t}+\underbrace{N^{-1}\sum_{i}^{N}(\widehat{r}_{k,i,t}-r_{k,i,t})+T_{k}^{-1}\sum_{t}^{T_{k}}(\widehat{r}_{k,i,t}-r_{k,i,t})+\sum_{a}\pi(a|o)(\widehat{r}_{k,i,t}-r_{k,i,t})}_{D_{k,i,t}},$
(C.14)

where

$\displaystyle
W_{k,i,t}=\sum_{k^{\prime}=1}^{k}\sum_{t=1}^{T_{k^{\prime}}}\varepsilon_{k^{\prime},i,t}/T_{k^{\prime}}+\sum_{k^{\prime}=1}^{k}\sum_{i=1}^{N}\varepsilon_{k^{\prime},i,t}/N,$
(C.15)

and Equation (C.14) follows from second iteration analysis, as such we have
$\bm{D}_{k}=\bm{\Phi}_{k}^{new}(\widehat{\bm{\beta}}_{k}-\bm{\beta}^{*}_{k})$.

By equation (C.4), we have at the step $K$, the error for the linear sieve
estimator $\widehat{\bm{\beta}}_{K}$ is

$\displaystyle\widehat{\bm{\beta}}_{K}-\bm{\beta}^{*}_{K}=\underbrace{\widehat{\bm{\Sigma}}_{K}^{-1}\bm{\Phi}_{K}^{\top}\bm{S}_{K}}_{\widehat{\bm{\Pi}}_{K}}\underbrace{(\bm{M}_{K}-\bm{\Phi}_{K}\bm{\beta}_{K}^{*})}_{\mbox{bias
term
}\bm{B}_{K}}+\widehat{\bm{\Sigma}}_{K}^{-1}\bm{\Phi}_{K}^{\top}\bm{S}_{K}\widetilde{\bm{\varepsilon}}_{K}$
$\displaystyle=$
$\displaystyle\;\widehat{\bm{\Pi}}_{K}\bm{B}_{K}+\widehat{\bm{\Pi}}_{K}\widetilde{\bm{\varepsilon}}_{K},$

again, the $\widetilde{\bm{\varepsilon}_{K}}$ is the adjusted error term, and
$K$ is used to address the $K$th step analysis. By Equation (C.14), the
adjusted error term can be written as

$\displaystyle\widetilde{\varepsilon}_{K,i,t}=\varepsilon_{K,i,t}+E_{K-1,i,t}=\varepsilon_{K,i,t}+W_{k-1,i,t}+D_{k-1,i,t}.$

It follows that

$\displaystyle\widehat{\bm{\beta}}_{K}-\bm{\beta}^{*}_{K}=\widehat{\bm{\Pi}}_{K}\bm{B}_{K}+\widehat{\bm{\Pi}}_{K}\widetilde{\bm{\varepsilon}}_{K}=\widehat{\bm{\Pi}}_{K}\bm{B}_{K}+\widehat{\bm{\Pi}}_{K}(\underbrace{\bm{\varepsilon}_{K}+\bm{W}_{K-1}}_{\bm{F}_{K}}+\bm{D}_{K-1})$
$\displaystyle=$
$\displaystyle\;\widehat{\bm{\Pi}}_{K}\bm{B}_{K}+\widehat{\bm{\Pi}}_{K}\bm{F}_{K}+\widehat{\bm{\Pi}}_{K}\bm{D}_{K-1}=\widehat{\bm{\Pi}}_{K}\bm{B}_{K}+\widehat{\bm{\Pi}}_{K}\bm{F}_{K}+\widehat{\bm{\Pi}}_{K}\bm{\Phi}_{K-1}^{new}(\widehat{\bm{\beta}}_{K-1}-\bm{\beta}^{*}_{K-1})$
$\displaystyle=$
$\displaystyle\;\widehat{\bm{\Pi}}_{K}\bm{B}_{K}+\widehat{\bm{\Pi}}_{K}\bm{\Phi}_{K-1}^{new}\widehat{\bm{\Pi}}_{K-1}\bm{B}_{K-1}+\widehat{\bm{\Pi}}_{K}\bm{F}_{K}+\widehat{\bm{\Pi}}_{K}\bm{\Phi}_{K-1}^{new}\widehat{\bm{\Pi}}_{K-1}\bm{F}_{K-1}+\widehat{\bm{\Pi}}_{K}\bm{\Phi}_{K-1}^{new}\widehat{\bm{\Pi}}_{K-1}\bm{D}_{K-2}$
$\displaystyle=$
$\displaystyle\;\widehat{\bm{\Pi}}_{K}\bm{B}_{K}+\widehat{\bm{\Pi}}_{K}\bm{\Phi}_{K-1}^{new}\widehat{\bm{\Pi}}_{K-1}\bm{B}_{K-1}+\widehat{\bm{\Pi}}_{K}\bm{\Phi}_{K-1}^{new}\widehat{\bm{\Pi}}_{K-1}\bm{\Phi}_{K-2}^{new}\widehat{\bm{\Pi}}_{K-2}\bm{B}_{K-2}$
$\displaystyle+\;\widehat{\bm{\Pi}}_{K}\bm{F}_{K}+\widehat{\bm{\Pi}}_{K}\bm{\Phi}_{K-1}^{new}\widehat{\bm{\Pi}}_{K-1}\bm{F}_{K-1}+\widehat{\bm{\Pi}}_{K}\bm{\Phi}_{K-1}^{new}\widehat{\bm{\Pi}}_{K-1}\bm{\Phi}_{K-2}^{new}\widehat{\bm{\Pi}}_{K-2}\bm{F}_{K-2}$
$\displaystyle+\;\widehat{\bm{\Pi}}_{K}\bm{\Phi}_{K-1}^{new}\widehat{\bm{\Pi}}_{K-1}\bm{\Phi}_{K-2}^{new}\widehat{\bm{\Pi}}_{K-2}\bm{D}_{K-3}$
$\displaystyle=$ $\displaystyle\;\dots\dots\dots\dots$ $\displaystyle=$
$\displaystyle\;\underbrace{\widehat{\bm{\Pi}}_{K}\bm{B}_{K}+\dots+\widehat{\bm{\Pi}}_{K}\bm{\Phi}_{K-1}^{new}\widehat{\bm{\Pi}}_{K-1}\dots\bm{\Phi}_{1}^{new}\widehat{\bm{\Pi}}_{1}\bm{B}_{1}}_{\mbox{systematic
bias term $I_{K,1}$}}$
$\displaystyle+\;\underbrace{\widehat{\bm{\Pi}}_{K}\bm{F}_{K}+\dots+\widehat{\bm{\Pi}}_{K}\bm{\Phi}_{K-1}^{new}\widehat{\bm{\Pi}}_{K-1}\dots\bm{\Phi}_{1}^{new}\widehat{\bm{\Pi}}_{1}\bm{F}_{1}}_{\mbox{previous
estimation error of fixed effects $I_{K,2}$}}$
$\displaystyle+\;\underbrace{\widehat{\bm{\Pi}}_{K}\bm{\Phi}_{K-1}^{new}\widehat{\bm{\Pi}}_{K-1}\dots\bm{\Phi}_{1}^{new}\widehat{\bm{\Pi}}_{1}\bm{D}_{1}}_{\mbox{previous
estimation error of Q-function $I_{K,3}$}}.$ (C.16)

By the Assumption 4 that
$\lVert[\mathbb{E}(\bm{\Phi}_{k}^{\top}\bm{S}_{k}\bm{\Phi}_{k})]^{-1}\mathbb{E}(\bm{\Phi}_{k}^{\top}\bm{S}_{k}\bm{\Phi}_{k-1}^{new})\rVert_{2}\leq\rho_{1}$
for any $k$, we also have
$\lVert\widehat{\bm{\Pi}}_{k}\bm{\Phi}_{k-1}^{new}\rVert_{2}<\rho_{1}$ with
high probability. Combined this with Lemma 2 and the first stage analysis, we
have $\lVert I_{K,1}\rVert_{2}=O_{p}(L^{-p/d})$, $\lVert
I_{K,2}\rVert_{2}=O_{p}(\sqrt{L/n_{K}})$. As for $I_{K,3}$, we have $\lVert
I_{K,3}\rVert_{2}=O_{p}\left(\rho_{1}^{K}\sqrt{L/n_{K}}\right)$ (this can be
similarly obtained by the argument in Equation (C.13)). Therefore, $\lVert
I_{K,3}\rVert_{2}=O_{p}\left(\sqrt{L/n_{K}}\right)$ if $K$ is bounded, and
$\lVert I_{K,3}\rVert_{2}=O_{p}\left(\rho_{1}^{K}\right)$ if $K$ diverges.
Thus in $K$ step, we still yield the same convergence rate for linear sieve
estimator.

####  C.4.6 Kth iteration analysis: fixed effects estimator

By Equations C.14 and C.15 and use the previous results, similarly, we have

$\displaystyle\left\lvert\widehat{\theta}_{K,i}-\theta_{K,i}\right\rvert\preceq\left\lvert\sum_{k=1}^{K}\sum_{t=1}^{T_{k}}\varepsilon_{k,i,t}/T_{k}\right\rvert.$
(C.17)

Now let $V_{k,i,t}=\sum_{t=1}^{T_{k}}\varepsilon_{k,i,t}/T_{k}$, apply Lemma
2, and Azuma Hoeffding’s inequality, and using the same argument as in first
iteration analysis, it yields the following:

$\displaystyle\mathbb{P}\left(\left\lvert\widehat{\theta}_{K,i}-\theta_{K,i}\right\rvert\geq\epsilon\right)=\mathbb{P}\left(\left\lvert\sum_{k=1}^{K}V_{k,i,t}\right\rvert\geq\epsilon\right)\leq
c/n,$ (C.18)

by setting $\epsilon=c\sqrt{T\log n}$. Combined Expression C.18 with
Bonferroni’s inequality, we immediately obtain that

$\displaystyle\max_{i}|\widehat{\theta}_{K,i}-\theta_{K,i}|=O_{p}(\sqrt{\log
n/T})+O_{p}(L^{-p/d}).$ (C.19)

One can similarly follow the way we bound
$\max_{i}|\widehat{\theta}_{K,i}-\theta_{K,i}|$ to bound the quantity
$\max_{t}|\widehat{\lambda}_{K,t}-\lambda_{K,t}|$, and it yields

$\displaystyle\max_{t}|\widehat{\lambda}_{K,t}-\lambda_{K,t}|=O_{p}(\sqrt{\log
n/N})+O_{p}(L^{-p/d}).$ (C.20)

####  C.4.7 Properties of the Final Value Estimators

Consider $\\{\widehat{\eta}_{i,t^{*}}\\}_{i,t^{*}}$ first. We aim to bound
$|\int_{o\in\mathcal{O}}[\widehat{\eta}_{i,t^{*}}^{\pi}(o)-\eta_{i,t^{*}}^{\pi}(o)]do|$,
or equivalently,

$\displaystyle\left\lvert\int_{o\in\mathcal{O}}\sum_{a}\pi(a|o)\left(\widehat{Q}_{1,1,t^{*}}(o,a)-{Q}_{i,1,t^{*}}(o,a)\right)\right\rvert.$

We first consider the scenario that $t^{*}\to\infty$, in which we can just set
$K=c\log NT$, and there is no need to do $t^{*}$ iterations to obtain the
final estimator since the approximation error between ${Q}_{i,1,t^{*}}(o,a)$
and ${Q}_{i,t^{*}-K+1,t^{*}}(O_{i,t^{*}-K+1},a)$ is order of $O(\pi_{0}^{K})$,
hence is negligible. Thus it suffices to bound the term

$\displaystyle\left\lvert\int_{o\in\mathcal{O}}\sum_{a}\pi(a|o)\left(\widehat{Q}_{1,t^{*}-K+1,t^{*}}(o,a)-{Q}_{i,t^{*}-K+1,t^{*}}(o,a)\right)\right\rvert$
$\displaystyle=$
$\displaystyle\;\left\lvert\widehat{\theta}_{K,i}-\theta_{K,i}+\widehat{\lambda}_{K,t}-\lambda_{K,t}+\bm{\varphi}^{\top}(\widehat{\bm{\beta}}_{K}-\bm{\beta}_{K}^{*})\right\rvert$

where
$\bm{\varphi}=\int_{o\in\mathcal{O}}(\pi(0|o)\bm{\Phi}_{L}^{\top}(o),\pi(1|o)\bm{\Phi}_{L}^{\top}(o))^{\top}do\in\mathbb{R}^{2L}$.

Recall that by Equations (C.19) and (C.20),

$\displaystyle\max_{i}|\widehat{\theta}_{K,i}-\theta_{K,i}|=O_{p}(\sqrt{\log
n/T})+O_{p}(L^{-p/d}),\mbox{ and }$
$\displaystyle\max_{t}|\widehat{\lambda}_{K,t}-\lambda_{K,t}|=O_{p}(\sqrt{\log
n/N})+O_{p}(L^{-p/d}).$

In addition, since $\lVert\bm{\varphi}\rVert_{2}$ is bounded by a constant
following Assumption 2, the term
$|\bm{\varphi}^{\top}(\widehat{\bm{\beta}}_{K}-\bm{\beta}_{K}^{*})|$ is
therefore order of $O_{p}(L^{-p/d})+O_{p}(\sqrt{\log n/n})$ (following the
same reasoning as we bound Equation (C.8) earlier). Thus,

$\max_{i,t}|\widehat{\eta}_{i,t}-\eta_{i,t}|=O_{p}(\sqrt{\log
n/T})+O_{p}(\sqrt{\log n/N})+O_{p}(L^{-p/d}).$

We now summarize the properties of $\widehat{\eta}_{i}^{\pi}$,
$\widehat{\eta}_{t}^{\pi}$, and $\widehat{\eta}^{\pi}$:

$\displaystyle\max_{t}|\widehat{\eta}_{t}^{\pi}-\eta_{t}^{\pi}|=O_{p}(N^{-1/2}\sqrt{\log
n})+O_{p}(L^{-p/d}),$
$\displaystyle\max_{i}|\widehat{\eta}_{i}^{\pi}-\eta_{i}^{\pi}|=O_{p}(T^{-1/2}\sqrt{\log
n})+O_{p}(L^{-p/d}),$ $\displaystyle\mbox{and
}|\widehat{\eta}^{\pi}-\eta^{\pi}|=O_{p}(\sqrt{\log n/n})+O_{p}(L^{-p/d}).$

As these three estimators are all aggregated over time and/or population, they
will converge at least as fast as $\widehat{\eta}_{i,t}^{\pi}$, while the bias
term $L^{-p/d}$ is unavoidable. In what follows, we only provide the proof for
bounding $\max_{t}|\widehat{\eta}_{t}^{\pi}-\eta_{t}^{\pi}|$ to save the
space, since the other two can be proved using similar arguments. By (C.17),
it follows that

$\displaystyle|\widehat{\eta}_{t}^{\pi}-\eta_{t}^{\pi}|=N^{-1}\left\lvert\sum_{i=1}^{N}(\widehat{\eta}_{i,t}^{\pi}-\eta_{i,t}^{\pi})\right\rvert$
$\displaystyle\preceq\;$ $\displaystyle
N^{-1}\left\lvert\sum_{i=1}^{N}(\widehat{\theta}_{K,i}-\theta_{K,i})\right\rvert+\left\lvert(\widehat{\lambda}_{K,t}-\lambda_{K,t})\right\rvert$
$\displaystyle\leq\;$ $\displaystyle
N^{-1}\left\lvert\sum_{i=1}^{N}\sum_{k=1}^{K}\sum_{t=1}^{T_{k}}\varepsilon_{k,i,t}/T_{k}\right\rvert+N^{-1}\left\lvert\sum_{i=1}^{N}\sum_{k=1}^{K}\varepsilon_{k,i,t}\right\rvert+O_{p}(L^{-p/d}).$

Thus, we have

$\displaystyle\max_{t}|\widehat{\eta}_{t}^{\pi}-\eta_{t}^{\pi}|\leq
N^{-1}\left\lvert\sum_{i=1}^{N}\sum_{k=1}^{K}\sum_{t=1}^{T_{k}}\varepsilon_{k,i,t}/T_{k}\right\rvert+\max_{t}|\widehat{\lambda}_{k,t}-\lambda_{k,t}|$
$\displaystyle=\;$ $\displaystyle O_{p}(N^{-1/2}\sqrt{\log
n})+O_{p}(L^{-p/d}).$

The scenario in which the time point $t^{*}$ is bounded is trivial, in this
case we have $K=t^{*}$, and this yields the same result. The proof is hence
completed.

##  Appendix D Proof of Theorem 2

Proof: We first establish the asymptotic property of $\widehat{\eta}^{\pi}$.
As this estimator is aggregated over both time and subjects, it suffices to
just focus on the main effects term, since we can always posit the constraint
that
$\sum_{i}(\widehat{\theta}_{k,i}-\theta_{k,i})=\sum_{t}(\widehat{\lambda}_{k,t}-\lambda_{k,t})=0$.
In what follows, we assume that $n_{K}$ is the same order of $n$ (recall that
this can be achieved by simply setting $K=c\log NT$), so that sometimes we can
use $n_{K}$ and $n$ exchangeablily.

By Equation (C.16), we have

$\displaystyle\widehat{\bm{\beta}}_{K}-\bm{\beta}^{*}_{K}=\underbrace{\widehat{\bm{\Pi}}_{K}\bm{B}_{K}+\dots+\widehat{\bm{\Pi}}_{K}\bm{\Phi}_{K-1}^{new}\widehat{\bm{\Pi}}_{K-1}\dots\bm{\Phi}_{1}^{new}\widehat{\bm{\Pi}}_{1}\bm{B}_{1}}_{\mbox{systematic
bias term $I_{K,1}$}}$
$\displaystyle+\;\underbrace{\widehat{\bm{\Pi}}_{K}\bm{F}_{K}+\dots+\widehat{\bm{\Pi}}_{K}\bm{\Phi}_{K-1}^{new}\widehat{\bm{\Pi}}_{K-1}\dots\bm{\Phi}_{1}^{new}\widehat{\bm{\Pi}}_{1}\bm{F}_{1}}_{\mbox{previous
estimation error of fixed effects $I_{K,2}$}}$
$\displaystyle+\;\underbrace{\widehat{\bm{\Pi}}_{K}\bm{\Phi}_{K-1}^{new}\widehat{\bm{\Pi}}_{K-1}\dots\bm{\Phi}_{1}^{new}\widehat{\bm{\Pi}}_{1}\bm{D}_{1}}_{\mbox{previous
estimation error of Q-function $I_{K,3}$}}.$

By the assumption that $NT\ll L^{2p/d}$, $\sqrt{NT}\lVert
I_{K,1}\rVert=o_{p}(1)$, it suffices to focus on $\sqrt{NT}(I_{K,2}+I_{K,3})$.
We first take a look at term $\widehat{\bm{\Pi}}_{K}\bm{F}_{K}$, recall that
$\widehat{\bm{\Pi}}_{K}=\widehat{\bm{\Sigma}}_{K}^{-1}\bm{\Phi}_{K}^{\top}\bm{S}_{K}$,
$\bm{F}_{K}=\bm{\varepsilon}_{K}+\bm{W}_{K-1}$, and $\bm{W}_{K}$ is
constructed by (C.15). Now let
$\widetilde{\bm{\Pi}}_{K}=\bm{\Sigma}_{K}^{-1}\bm{\Phi}_{K}^{\top}\bm{S}_{K}$,
by Lemma E.2 of Shi et al., (2022), it can be shown that
$\lVert\widehat{\bm{\Sigma}}_{K}^{-1}-\bm{\Sigma}_{K}^{-1}\rVert_{2}=O_{p}(n_{K}^{-1}\sqrt{L/n_{K}}\log
n_{K}).$ Thus we have
$\sqrt{n_{k}}\lVert(\widetilde{\bm{\Pi}}_{K}-\widehat{\bm{\Pi}}_{K})\bm{F}_{K}\rVert_{2}=o_{p}(1)$,
using the Assumption 2 (iii) that $L\ll\sqrt{N+T}/\log(N+T)$. It follows that

$\displaystyle\widetilde{\bm{\Pi}}_{K}\bm{F}_{K}=\bm{\Sigma}_{K}^{-1}\bm{\Phi}_{K}^{\top}\bm{S}_{K}\bm{F}_{K}=\bm{\Sigma}_{K}^{-1}\bm{\Phi}_{K}^{\top}\bm{F}_{K}-\bm{\Sigma}_{K}^{-1}\bm{\Phi}_{K}^{\top}\bm{B}_{K}(\bm{B}_{K}^{\top}\bm{B}_{K})^{+}\bm{B}_{K}^{\top}\bm{F}_{K},$

now let
$\widehat{\bm{\Omega}}_{K}=\bm{\Phi}_{K}^{\top}\bm{B}_{K}(\bm{B}_{K}^{\top}\bm{B}_{K})^{+}$,
and
$\bm{\Omega}_{K}=\mathbb{E}(\bm{\Phi}_{K}^{\top}\bm{B}_{K}(\bm{B}_{K}^{\top}\bm{B}_{K})^{+})$.
Similarly we have
$\lVert\bm{\Omega}_{K}-\widehat{\bm{\Omega}}_{K}\rVert_{2}=O_{p}((N+T_{K})^{-1}\sqrt{L/(N+T_{K})}\log(N+T_{K}))$,
again we can replace $\widehat{\bm{\Omega}}_{K}$ with its expectation in the
above formula, so now it yields that

$\displaystyle\sqrt{n_{k}}\widehat{\bm{\Pi}}_{K}\bm{F}_{K}=\sqrt{n_{k}}(\bm{\Sigma}_{K}^{-1}\bm{\Phi}_{K}^{\top}\bm{F}_{K}-\bm{\Sigma}_{K}^{-1}\bm{\Omega}_{K}\bm{B}_{K}^{\top}\bm{F}_{K})+o_{p}(1).$
(D.1)

We first deal with $\bm{\Phi}_{K}^{\top}\bm{F}_{K}$,

$\displaystyle\bm{\Phi}_{K}^{\top}\bm{F}_{K}=\bm{\Phi}_{K}^{\top}\bm{\varepsilon}_{K}+\bm{\Phi}_{K}^{\top}\bm{W}_{K-1}$
$\displaystyle=$
$\displaystyle\;\sum_{i^{\prime},t^{\prime}}\left[\bm{\Phi}_{K,i^{\prime},t^{\prime}}\varepsilon_{K,i^{\prime},t^{\prime}}+\left(\sum_{i}\bm{\Phi}_{K,i,t^{\prime}}/N+\sum_{t}\bm{\Phi}_{K,i^{\prime},t}/T_{K}\right)\sum_{k}^{K-1}\varepsilon_{k,i^{\prime},t^{\prime}}\right],$

again, we can replace
$\sum_{i}\bm{\Phi}_{K,i,t^{\prime}}/N+\sum_{t}\bm{\Phi}_{K,i^{\prime},t}/T_{K}$
in the above formula to its expectation
$\widetilde{\bm{\Phi}}_{K,i^{\prime},t^{\prime}}\equiv\mathbb{E}(\sum_{i}\bm{\Phi}_{K,i,t^{\prime}}/N+\sum_{t}\bm{\Phi}_{K,i^{\prime},t}/T_{K})$.
As such, Equation (D.1) now becomes