The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
210,129 characters
Post Reinforcement Learning Inference
\RUNTITLE{Post Reinforcement Learning Inference}
\TITLE{Post Reinforcement Learning Inference}
\ARTICLEAUTHORS{
\AUTHOR{Vasilis Syrgkanis\thanks{Stanford University, \EMAIL{[email removed]}. Vasilis Syrgkanis was supported by NSF Award IIS-2337916.}\quad Ruohan Zhan\thanks{University College London, \EMAIL{[email removed]}}}
\RUNAUTHOR{Syrgkanis and Zhan}
}
\ABSTRACT{
We study estimation and inference using data collected by reinforcement learning (RL) algorithms. These algorithms adaptively experiment by interacting with individual units over multiple stages, updating their strategies based on past outcomes. Our goal is to evaluate a counterfactual policy after data collection and estimate structural parameters, such as dynamic treatment effects, that support credit assignment and quantify the impact of early actions on final outcomes. These parameters can often be defined as solutions to moment equations, motivating moment-based estimation methods developed for static data. In RL settings, however, data are often collected adaptively under nonstationary behavior policies. As a result, standard estimators fail to achieve asymptotic normality due to time-varying variance. We propose a weighted
generalized method of moments (GMM)
approach that uses adaptive weights to stabilize this variance. We characterize weighting schemes that ensure consistency and asymptotic normality of the weighted
GMM estimators,
enabling valid hypothesis testing and uniform confidence region construction. Key applications include dynamic treatment effect estimation and dynamic off-policy evaluation.
}
\KEYWORDS{reinforcement learning,
GMM estimators,
adaptive weighting, asymptotic normality, strong Gaussian approximation, hypothesis testing, dynamic treatment effects, dynamic off-policy evaluation }
\maketitle
\section{Introduction}
Adaptive data collection has become a staple of the digital economy. Most major digital platforms invoke adaptive experimentation algorithms to optimize their service. Adaptive experimentation allows one to progressively update their experimentation strategy and leads to efficient sample usage \citep{chu2011contextual,agrawal2013thompson} and efficient use of the experimentation budget for increased hypothesis testing power \citep{russo2016simple}. Moreover, frequent adaptive experimentation is starting to become adopted in other domains such as personalized healthcare \citep{murphy2005experimental,offer2021adaptive}.
The popularity of adaptive experiments has increased the availability of data collected from such designs. However, adaptive data collection raises many new research challenges, especially in the case of post-collection statistical analysis. For instance, constructing confidence intervals and testing statistical significance for alternative candidate policies or for structural parameters like average treatment effects, from adaptively collected data, has been shown to be theoretically challenging.
Prior work has mostly focused on the bandit setup, where at each time the experimenter only interacts with units sampled from the environment once and observes the immediate outcome~\citep{deshpande2018accurate,zhang2021statistical,hadad2021confidence,bibaut2021post,zhan2021off}. This unfortunately cannot accommodate many applications including adaptive clinical trials and dynamic treatment regimes, where a patient often receives multiple rounds of treatments to improve the outcome~\citep{murphy2003optimal,lei2012smart}, or a digital platform interacts with their users over a sequence of multiple page visits~\citep{chen2019top}.
The goal of our work is to address this largely un-explored area of post-adaptive data collection inference, from such adaptive experiments that, within each experiment phase, involve multiple interactions with the same treated unit. In particular, we consider data that are collected from reinforcement learning~(RL) algorithms \citep{sutton1998introduction} and provide estimation and inference for structural parameters of interest under semi-parametric assumptions~\citep{neyman1979c,laan2003unified,chernozhukov2022locally}.
In RL, an adaptive experimentation algorithm (from now on ``the agent'') interacts with units that are sampled independently and identically (i.i.d.) from the environment. For each unit, the agent observes an initial state, then applies multiple treatments that cause state transitions, and finally observes an outcome at the end.\footnote{We focus on cases where only the final outcome is observed. Our framework can be generalized to settings where intermediate outcomes are also revealed.} We term the sequence of interactions with a unit as an ``episode''. For example, consider an educational platform aiming to encourage users to enroll in a course. The platform, acting as the agent, might first interact with a user on the homepage, then guide them to an enrollment page, and finally to a payment page, with each step representing a state transition. Throughout these stages, the platform can experiment with different page designs (treatments) to achieve the final goal of course enrollment.
The agent's behavior policy is typically adapted over time to incorporate learning from past realizations. In particular, we focus on a broad class of adaptive data collection mechanisms, where the agent interacts with a batch of units before updating its policy for subsequently arriving units. This setting encompasses the \emph{episodic reinforcement learning} framework \citep{neu2020unifying} and generalizes the classical bandit framework, in which the agent interacts with each unit for only a single round in an adaptive manner.
Such adaptivity progressively improves the agent's performance but results in a nonstationary behavior policy and introduces dependence between observations.
As a result, we cannot simply view the data from each episode as i.i.d.~samples and pass them to traditional estimation pipelines \citep{lewis2020double}.
Even when employing estimation techniques designed to address time-series correlation, these often require stationarity—a condition not met by most adaptive experimentation processes.
This nonstationarity causes evolving discrepancy between the behavior and target policies, an issue known as changing ``overlap'' and resulting in varied estimation variances for sequentially collected samples \citep{imbens2004nonparametric,hadad2021confidence}.
Therefore, averaging samples uniformly can be suboptimal, resulting in significant variance and a non-normal asymptotic distribution, which complicates post-experimental inference.
This problem, evident even in single-interaction bandit scenarios, has prompted recent literature to suggest re-weighting samples to stabilize time-varying variance, such that the resulting estimators are consistent and asymptotically normal~\citep{deshpande2018accurate,hadad2021confidence,zhang2021statistical,zhan2021off,bibaut2021post}.
RL settings further complicate the estimation problem in two ways.
First, exogenous random shocks during state transitions within an episode affect subsequent states and are correlated with future treatments; these shocks, being unobservable, cause unmeasured confounding~\citep{robins1986new,robins2004optimal,chakraborty2013semi}.
To address the identification issue, we follow the semi-parametric inference literature, formulating the dynamic treatment effect estimation problem as estimating the structural parameters in a structural nested mean model~\citep{robins2004optimal,lok2012impact,vansteelandt2016revisiting}.
In our model, each stage is linked to a specific structural parameter that quantifies the treatment effect at that stage. We demonstrate that these parameters are the solutions to stage-wise moment equations, derived through $g$-estimation and constructed in reverse order, from the last stage to the first. Given the nonstationarity of RL data, which often lead to time-varying estimates across units as discussed above, we apply nonuniform and adaptive weights to the unit samples to stabilize these variances.
Second, estimation with RL data often faces the ``curse of horizon'' challenge, where the overlap between behavior and target policies deteriorates across episodic stages, leading to accumulated estimation variance.
Prior work suggests heuristics like weight clipping, which---while effective in controlling variance---introduce a small bias \citep{precup2000eligibility,chen2019top}.
However, our approach, through nuanced modeling, ensures that moment equations across stages remain uncorrelated with zero covariance. This allows for the application of stage-wise adaptive weights, which stabilize the variance at each stage based on information available up to that point, effectively circumventing the problem of variance accumulation.
This work also enriches the inference literature when using adaptively collected data by providing results for GMM-estimation. Prior work mostly addresses inference in the context of M-estimation \citep{deshpande2018accurate,zhang2021statistical}, i.e. parameters that can be defined as the minimizers of a population loss function. However, structural parameters in problems defined by moment conditions, which often arise in the dynamic treatment regime settings, cannot be phrased as M-estimators (rather they can be better thought as instrumental variable problems).
\subsection{Our Contributions}
Our main contributions are outlined as follows.
First, we propose a weighted GMM-estimator using adaptive RL data to estimate structural parameters of counterfactual policy values and dynamic treatment effects.
These parameters are defined by moment conditions within a dynamic treatment regime, modeled through a semiparametric structure nested mean model for each unit's episodic data.
Specifically, for a target counterfactual policy $\pi^*$, we aim to estimate its policy value $\theta_0^*$---the expected final outcome $Y$ under policy $\pi$. We further define $\theta_j^*$ as the stage-wise dynamic treatment effect of treatment $T_j$ at stage $j$, reflecting the expected change in the final outcome $Y$ due to $T_j$, assuming future treatments follow $\pi^*$. Let $L$ be total number of stages in an episode for each unit. Our objective is to derive the structural parameter set $\theta^*=(\theta_0^*, \theta_1^*,\dots, \theta_L^*)$ from RL data, capturing both the policy value and stage-wise dynamic treatment effects.
We hereby illustrate the core concept of our approach under the episodic RL setup.
We consider $n$ sequentially collected episodes $\{Z_1, \dots, Z_n\}$. Let $\mathcal{F}_i:=\sigma(Z_{1:i-1})$ represent the $\sigma$-algebra with respect to which the random variables $Z_{1:i-1}$ are measurable.
The true parameter vector $\theta^*$ satisfies the moment condition $\mathbb{E}[m(Z_i;\theta^*) \mid \mathcal{F}_i] = 0$, where $m(\cdot)$ denotes bounded linear moment functions derived from our semiparametric model. To enhance the stability and efficiency of the estimation, we adopt the generalized method of moments (GMM) throughout this paper.\footnote{The results presented here can be directly extended to Z-estimation.} The Jacobian of the moment function, denoted by $J_i := \nabla_{\theta} m(Z_i; \theta^*)$, is nonstationary across episodes due to adaptivity. As a result, the traditional (unweighted) GMM estimator, which is based on minimizing the norm of empirical moment equations, may exhibit unstable asymptotic behavior. To address this, we introduce adaptive weights $H_i$, measurable with respect to $\mathcal{F}_i$, to stabilize the estimation.
Let $\hat{\theta}_n$ denote the weighted GMM estimator that minimizes $\left\|\frac{1}{n}\sum_{i=1}^n H_i m (Z_i;\theta)\right\|_A^2$, where $A$ is a positive-definite matrix and $\|\cdot\|_A$ is the norm induced by $A$. When $\hat{\theta}_n$ lies in the interior of $\Theta$, where $\Theta$ is the bounded parameter space of $\theta^*$, it satisfies the first-order condition: $\{\frac{1}{n}\sum_{i=1}^n H_iJ_i\}A\{\frac{1}{n}\sum_{i=1}^n H_im (Z_i;\hat\theta_n)\}=\mathbf{0}$.
By Taylor expansion, we have:
\begin{equation}
\label{eq:intro}
B_n^\top A B_n\, (\hat{\theta}_n-\theta^*)=B_n^\top A S_n, \quad\text{where}\quad B_n=-\frac{1}{n}\sum_{i=1}^nH_i\,J_i, \quad S_n =\frac{1}{n}\sum_{i=1}^n H_i\,m(Z_i;\theta^*).
\end{equation}
The property of $\hat{\theta}_n$ thus depends on the behavior of $B_n$ and $S_n$, both of which can be stabilized by the choice of weights $\{H_i\}$.
Our second contribution is to identify generic weighting schemes to achieve consistency and asymptotic normality of the weighted GMM-estimator $\hat{\theta}_n$.
The choice of weights $\{H_i\}$ should address two critical aspects: ensuring the normalizing matrix $B_n^\top A B_n$ in \eqref{eq:intro} remains well-posed to prevent the explosion of the estimated parameter $\hat{\theta}_n$, and the right-hand side $B_n^\top A S_n$ converges to a Gaussian distribution for asymptotic normality.
We show that $S_n$
is the sum of a martingale difference sequence, and we establish its convergence by leveraging recent advances in martingale limiting theorems under a crucial homoscedasticity assumption.
We note that simply providing a central limit theorem for \eqref{eq:intro} is not enough to achieve inference results on policy value $\theta_0^*$ or dynamic treatment effect $\theta_j^*$, since the normalizing matrix $B_n^\top A B_n$ is likely to not concentrate by the sheer nature of adaptivity.
Addressing this, we adapt recent work in \cite{cattaneo2022yurinskii}, which develops strong Gaussian approximation results for martingale data, to our context and show that \eqref{eq:intro} converges to a Gaussian distribution at a uniform rate. This allows us to invert the normalizing matrix $B_n^\top A B_n$ and approximate the estimation error, $\hat{\theta}_n-\theta^*$, with a Gaussian distribution, paying the way for inference results that are the primary focus of this paper.
Our third contribution is to substantiate the generic weighting schemes and offer specific weighting choices to achieve consistency and asymptotic normality for a wide range of RL algorithms with polynomially decaying exploration rates.
We focus on bilinear parametrization, which includes common setups such as categorical treatment and polynomial scalar treatment. We show that the adaptive weights for consistency can be directly derived from the collected data.
Further, we identify the ``oracle weights'' to achieve strong Gaussian approximation. While these weights
may require additional estimations,
the necessary quantities to be estimated only depend on the state transition dynamics, which are independent of the behavior policy; this allows for their consistent estimation through online regression techniques, which we term ``feasible weights''.
We prove that the GMM-estimator, when applied with these feasible weights,
preserves strong Gaussian approximation, thereby facilitating post-RL inference.
Finally, we apply our estimation and inference framework to high-dimensional Markovian models.
We offer estimation guarantees for feasible weights and show that the adaptively weighted GMM-estimator, when using these feasible weights, achieves strong Gaussian approximation. This approximation enjoys a uniform convergence rate of $O(n^{-(1-L\alpha)/10})$, where $L$ represents the episode length, and $\alpha$ indicates the exploration decaying rate of the RL agent.
We further substantiate our theoretical guarantees with numerical evidence.
Our findings show that the standard GMM-estimator, when unweighted, deviates from asymptotic normality, resulting in either insufficient coverage or overly conservative inferences for structural parameters. Conversely, our weighted GMM-estimator, applied with feasible weights, consistently achieves near-perfect coverage, tighter confidence intervals, and more precise estimates for evaluation policy values. Its performance closely matches that of estimations under oracle weights, which, though ideal, require inaccessible true structural parameters and are generally impractical.
Furthermore, we demonstrate our method's robustness under model misspecification, employing a polynomial model to approximate the true exponential model within our semiparametric framework. Even under misspecification, our weighted GMM-estimator with feasible weights outperforms the standard approach, offering more accurate policy value estimates and maintaining near-nominal coverage as approximation complexity increases. This robustness highlights our method's suitability for real-world scenarios, where model misspecification is often inevitable.
\section{Setup}
This section establishes our problem framework by detailing the data generating process. We start by describing the stochastic control process, which involves multiple stages of interaction with individual units. Subsequently, we define the structural parameter estimation problem within this framework.
\subsection{Episodic Potential Outcome}
We start with describing the stochastic control process to roll out an episode for a unit and the corresponding potential outcome.
Consider each unit starts from an initial state~$S_1$ i.i.d.~sampled from a fixed distribution $P_S$.
An RL agent then runs an episode of $L$ stages for each unit.
At each stage $j \in [1:L]$, the agent observes the unit’s current state $S_j \in \mathcal{S} \subset \mathbb{R}^{d_s}$ and then assigns a treatment $T_j \in \mathcal{T} \subset \mathbb{R}^{d_\tau}$ (also known as an action or intervention), where $d_s$ and $d_\tau$ denote the dimensions of the state and treatment, respectively. The treatment transitions the unit to a new state $S_{j+1}$ for the next stage. The long-term outcome of interest, denoted by $\ensuremath{Y} \in \mathcal{Y} \subset \mathbb{R}$, is observed only at the end of the episode $(S_1, T_1, S_2, \dots, S_L, T_L)$.
Following the potential outcome framework, we denote the potential outcome under a sequence of treatment assignments $\{\tau_1,\dots, \tau_L: \tau_j\in\mathcal{T}\}$ as $\ensuremath{Y}(\tau_{1:L})$.\footnote{We use $A_{n_1:n_2}$ as a short hand for the set $\{A_{n_1},\dots, A_{n_2}\}$ for any indexed variable $A_i$ and any two integers $n_1\leq n_2$.}
A dynamic treatment assignment policy $\pi = \{\pi_1, \dots, \pi_L\}$ specifies a sequence of treatment rules applied across stages of an episode. Specifically, each $\pi_j$ determines the treatment $T_j$ at stage $j$ based on the history $(S_{1:j}, T_{1:j-1})$, with $\pi_j : \mathcal{S}^j \times \mathcal{T}^{j-1} \rightarrow \mathcal{T}$. We overload the notation and use $\ensuremath{Y}(\pi_{1:L})$ to denote the potential outcome under policy $\pi$.
We introduce an assumption critical to our analysis: all confounders influencing both the treatment assignment and the final outcome are observed by the states.
\begin{assumption}[Sequential Conditional Exogeneity]
\label{assump:exogeneity}
The data generating process satisfies
\begin{equation*}
\ensuremath{Y}(T_{1:j-1}, \pi_{j:L})\, {\perp\!\!\!\perp}\, T_j \mid T_{1:j-1}, S_{1:j}, \quad \mbox{for any}~ j=1,\dots,L.
\end{equation*}
\end{assumption}
This assumption is an extension of the unconfoundedness assumption in causal inference literature to settings where a unit undergoes multiple stages of treatment \citep{rosenbaum1983central,robins2004optimal}.
\begin{remark}
A Markov decision process (as depicted in Figure~\ref{fig:cg}) satisfies Assumption~\ref{assump:exogeneity}. In contrast, a partially observable Markov decision process (POMDP) violates this assumption. By POMDP, we refer to settings where the actual treatment received depends on some latent state; this could occur if the behavior policy relies on unobserved information beyond $S_t$, or if there is non-compliance. However, such scenarios do not arise in our analysis framework, as we assume full knowledge of the behavior policy.
\end{remark}
\subsection{Inference Goal: Value of Dynamic Treatment Assignment Policy}
Given a policy $\pi^*$ to be evaluated, the goal in this paper is to estimate its policy value, which is defined as the mean outcome under this policy, outlined below:
\begin{equation*}
\theta^*_0 := \mathbb{E}\left[Y(\pi^*_{1:L})\right].
\end{equation*}
Inference on $\theta^*_0$ allows one to argue whether the observed outcome achieved under the behavior policy is statistically larger than the outcome under some simple baseline policies, and therefore whether the data collection agent (for example, from some RL algorithm) leads to any statistically significant benefits as compared to these baseline policies.
\begin{figure}
\centering
\begin{subfigure}[t]{0.45\textwidth}
\hspace*{-4.5cm}
\centering
\begin{tikzpicture}[scale=1.7]
\node (S) at (0,0) {$S_1$};
\node (T1) at (1,1) {$T_1$};
\node (X) at (2,0) {$S_{2}$};
\node (T2) at (3,1) {$T_{2}$};
\node (Y) at (4,0) {$Y$};
\draw[thick, ->] (S) -- (T1);
\draw[thick, ->] (S) -- (X);
\draw[thick, ->] (T1) -- (X);
\draw[thick, ->] (X) -- (Y);
\draw[thick, ->] (X) -- (T2);
\draw[thick, ->] (T2) -- (Y);
\hspace{10mm}
\end{tikzpicture}
\hspace{-3cm}
\caption{\footnotesize{Causal graph for observed data}}
\label{fig:cg}
\end{subfigure}
\begin{subfigure}[t]{0.45\textwidth}
\hspace*{-.7cm}
\centering
\begin{tikzpicture}[scale=1.7]
\node (S) at (0,0) {$S_1$};
\node (T1) at (0.67, 1) {$T_1$};
\node (pi1) at (1.67, 1) {$\pi_1(S_1)$};
\node (X) at (2.33,0) {$S_2(\pi_1)$};
\node (T2) at (2.8,1) {$T_2$};
\node (pi2) at (4,1) {$\pi_2(S_2(\pi_1))$};
\node (Y) at (4.63,0) {$Y(\pi_{1:2})$};
\draw[thick, ->] (S) -- (T1);
\draw[thick, ->] (S) -- (X);
\draw[thick, ->, dashed] (S) to[bend right] (pi1);
\draw[thick, ->, dashed] (pi1) -- (X);
\draw[thick, ->] (X) -- (Y);
\draw[thick, ->] (X) -- (T2);
\draw[thick, ->, dashed] (X) to[bend right] (pi2);
\draw[thick, ->, dashed] (pi2) -- (Y);
\end{tikzpicture}
\caption{\footnotesize{Intervention graph for counterfactual data}}
\label{fig:swig}
\end{subfigure}
\vspace{1em}
\caption{
Causal graphs illustrating the sequential conditional exogeneity assumption for a two-period setting, showing both observed data and counterfactual data for each unit under an alternative policy $\pi$.
Solid arrows represent the observed treatment assignments and state transitions in the collected data,
while dotted arrows indicate counterfactual treatment assignments and state transitions under policy $\pi$.
}
\end{figure}
We introduce the following definition to attribute the long-term outcome to intermediate treatments at each stage.
\begin{definition}[Blip function, \cite{robins2004optimal}]
Given an evaluation policy $\pi^*$, let the ``blip function'' $\gamma_j^{(\pi^*)}$ at stage $j$ be defined as:
\begin{align*}
\gamma_j^{(\pi^*)}(S_{1:j}, T_{1:j}) := \mathbb{E}\left[\ensuremath{Y}\left(T_{1:j}, \pi^*_{j+1:L}\right) - \ensuremath{Y}\left(T_{1:j-1}, 0,\pi^*_{j+1:L}\right) \mid S_{1:j}, T_{1:j}\right].
\end{align*}
\end{definition}
In other words, the blip effect $\gamma_j^{(\pi^*)}(\cdot)$, given the realized history $(S_{1:j}, T_{1:j-1})$ and conditional on assigned treatment $T_{i,j}$, captures the effect of assigning treatment $T_{i,j}$ in the current period relative to a baseline control policy $0 \in \mathcal{T}$ that applies the status quo treatment, assuming future treatments follow policy $\pi^*$. From now on,
we fix the evaluation policy $\pi^*$ and omit the superscript in the blip function, i.e. we use $\gamma_j(\cdot)$ to denote $\gamma_j^{(\pi^*)}(\cdot)$ for notation convenience.
\begin{remark}
This blip function $\gamma_j(\cdot)$ quantifies the treatment effect of $T_j$ at stage $j$ on the long-term outcome $Y$. This concept is analogous to the ``advantage'' function used in the RL literature, which captures the discrepancy in Q-value functions from taking alternative actions \citep{baird1995residual,shi2022statistically}.
\end{remark}
To identify the policy value of $\pi^*$ from the observed data, we introduce a semi-parametric structural assumption on the blip functions,
following \cite{lewis2020double,robins2004optimal}.
\begin{assumption}[Linear blip assumption]
\label{assump:linear_blip_function}
The blip function has the linear form:
$\gamma_j(S_{1:j}, T_{1:j}) = \phi_j( S_{1:j}, T_{1:j})^\top \theta_j^* $
for some unknown structural parameter $\theta_j^*\in\mathbb{R}^d$ and a known feature mapping function $\phi_j:\mathcal{S}^j\times \mathcal{T}^j\rightarrow\mathbb{R}^d$ that satisfies $\phi_j(s_{1:j}, (\tau_{1:j-1},0))=\mathbf{0}$ for any $s_{1:j}\in\mathcal{S}^j$ and $\tau_{1:j-1}\in\mathcal{T}^{j-1}$.
As shorthand notation, let $\Phi_j:=\phi_j(S_{1:j}, T_{1:j}) - \phi_j(S_{1:j}, T_{1:j-1}, \pi^*(S_{1:j}, T_{1:j-1}))$.
\end{assumption}
Note that even when the blip functions are assumed to be parametrically linear, we do not make any assumption on the conditional mean of the evaluation policy outcome $\mathbb{E}\left[\ensuremath{Y}(\pi^*_{1:L})\mid S_1\right]$, which can be potentially nonlinear in the initial state $S_1$. In later sections we
shall use moment equations to identify the unknown structural parameter $\theta_j$, and this assumption implies that the moment conditions are linear in target parameters.
This assumption allows for a wide range of model classes including high-dimensional Markovian models substantiated in Section \ref{sec:partial_linear_model}.
To this end, we can explicitly characterize the policy value $\theta^*_0$ for evaluation policy $\pi^*$, following results in \cite{robins2004optimal} -- which we adapt to our setup and provide its proof (along with all other proofs) in the appendix for completeness.\footnote{This result can be viewed as an analogue of the \emph{performance difference lemma} that is frequently used in the reinforcement learning literature \cite{kakade2002}.}
\begin{lemma}[Expected Outcome via Blip Functions]
\label{lemma:identification_policy_value}
Under Assumptions~\ref{assump:exogeneity}~\&~\ref{assump:linear_blip_function}, it holds that
\[
\mathbb{E}\left[\ensuremath{Y}(T_{1:j-1}, \pi^*_{j:L})\mid S_{1:j}, T_{1:j}\right] = \mathbb{E}\left[Y -\sum_{j'=j}^L\Phi_{j'}^\top \theta_{j'}^* \mid S_{1:j}, T_{1:j}\right].
\]
Hence, the evaluation policy value is $\theta^*_0= \mathbb{E}\left[\ensuremath{Y}-\sum_{j=1}^L\Phi_j^\top \theta_j^*\right]$.
\end{lemma}
\subsection{Research Question}
We conclude this section by formally describing our problem within the context formalized above.
\paragraph{Adaptive Episodic RL Data.}
We consider the setting where an RL agent sequentially rolls out episodes for $n$ units, each over $L$ stages. Before episode~$i$ begins, the agent updates the treatment assignment policy $\pi_i := \{\pi_{i,1}, \dots, \pi_{i,L}\}$ based on data from previous episodes. Let $Z_i := \{S_{i,1}, T_{i,1}, \dots, S_{i,L}, T_{i,L}, Y_i\}$ denote the data collected during episode~$i$. The distribution of $Z_i$ generally depends on the past realizations $Z_{1:i-1}$. The full dataset consists of observations $\{Z_i\}_{i=1}^n$. Our results directly extend to batched episodic RL settings, where the RL agent may interleave interactions with units within a batch and simultaneously progress these units through their respective episodes. This extension is possible because our technical analysis is based on filtrations defined with respect to unit-level episodes. As long as the behavior policy remains fixed during the rollout of each individual unit, similar filtrations can be constructed based on episode completion, and our results continue to hold.
We assume that the behavior policy $\pi_i$ used to generate data for episode~$i$ is known and recorded by the RL agent. Without this information, recovering the behavior policy from the data is generally difficult or even statistically impossible, since each $\pi_i$ may correspond to only one observed episode. This assumption is also common in prior work on inference with adaptively collected bandit data \citep{hadad2021confidence,zhang2021statistical}.
\paragraph{Goal.} Given the data, we aim to estimate and perform inference on the structural parameters in the parameter space $\Theta\subset \mathbb{R}^{1+dL}$:
\[\theta^* = (\theta_0^*, \theta_1^* ,\dots,\theta_L^*)\in\Theta,\]
where $\theta_0^*$ corresponds to the value of the evaluation policy $\pi^*$, and $\{\theta_j^* \in \mathbb{R}^d:j\in[1:L]\}$ denote the dynamic treatment effect parameters that appear in the linearization of the blip functions.
In particular, we assume that $\theta^*$ is an interior point of $\Theta$.
\paragraph{Notation.}
We use $i$ to index episode and $j$ to index stage within an episode. Define the $\sigma$-fields $\mathcal{F}_{i} := \sigma(\{Z_{1:i-1}\})$ and $\mathcal{F}_{i,j} := \sigma(\{Z_{1:i-1}, S_{i, 1:j-1}, T_{i,1:j-1}\})$ with the convention that $\mathcal{F}_{i,1}\equiv \mathcal{F}_{i}$ and $\mathcal{F}_{i,0}\equiv \mathcal{F}_{i}$.
We also introduce an augmented $\sigma$-field: $\mathcal{F}_{i,j}^+$, which extends $\mathcal{F}_{i,j}$ by including the $\sigma$-algebra generated by the next state $S_{i,j}$; this additional information about state $S_{i,j}$ is available prior to the treatment assignment by the agent for unit $i$ at stage $j$.
We use $\mathbb{E}_{i}[\cdot]$ as a short hand for the conditional expectation $\mathbb{E}[\cdot\mid \mathcal{F}_{i}]$ and use $\mathbb{E}_{i,j}[\cdot]$ for $\mathbb{E}[\cdot \mid \mathcal{F}_{i,j}]$, $\mathbb{E}_{i,j}^+[\cdot]$ for $\mathbb{E}[\cdot \mid \mathcal{F}_{i,j}^+]$.
Similarly, we define $\mathrm{Var}_{i}(\cdot), \mathrm{Var}_{i,j}(\cdot)$ and $\mathrm{Var}_{i,j}^+(\cdot)$ as shorthand for the conditional variances $\mathrm{Var}(\cdot\mid \mathcal{F}_{i}), \mathrm{Var}(\cdot\mid \mathcal{F}_{i,j})$ and $\mathrm{Var}(\cdot\mid \mathcal{F}_{i,j}^+)$, respectively. For a random vector $X$, we define $\mathrm{Var}_i(X)=\mathbb{E}_i[XX^\top] - \mathbb{E}_i[X]\mathbb{E}_i[X]^\top$, $\mathrm{Var}_{i,j}(X)=\mathbb{E}_{i,j}[XX^\top] - \mathbb{E}_{i,j}[X]\mathbb{E}_{i,j}[X]^\top$ and $\mathrm{Var}_{i,j}^+(X) = \mathbb{E}_{i,j}^+[XX^\top]-\mathbb{E}_{i,j}^+[X] \mathbb{E}_{i,j}^+[X]^\top$.
\section{Identification and Estimation}
In this section, we provide conditional moment equations for identifying the target parameters. We then propose our adaptively weighted method of moments (AW-GMM) for parameter estimation using the adaptively collected RL data.
\subsection{Identification via Moment Equations}
We construct moment equations to identify $\theta^*$, as formalized by the lemma below. This approach is known as the $g$-estimation framework, typically used in the context of structural nested mean models~\citep{robins2004optimal}.
However, most prior work focuses on data collected by a fixed behavior policy \citep{robins2004optimal,chakraborty2013semi,lewis2020double}, while here we apply the $g$-estimation approach to adaptive RL data and show that the true parameter vector $\theta^*$ also satisfies the moment restrictions proposed in \cite{robins2004optimal}.
\begin{lemma}[Identification of Blip Functions]
\label{lemma:identification_parameter}
For unit $i$ and stage $j$, define
\begin{equation}
\Phi_{i,j}:=\phi_j(S_{i, 1:j}, T_{i,1:j}) - \phi_j(S_{i, 1:j}, T_{i,1:j-1}, \pi^*(S_{i, 1:j}, T_{i,1:j-1})),
\end{equation}
and
\begin{equation}
\label{eq:residual}
R_{i,j}:= \ensuremath{Y}_i - \sum_{j'=j}^L \Phi_{i,j'}^\top \theta_{j'}^*.
\end{equation}
Here, $ R_{i,j}$ represents the \textbf{residual} by subtracting the blip effects of the future treatments $T_{i,j:L}$ from the final outcome~$\ensuremath{Y}_i$
and adding the blip effects of the treatments assigned by the evaluation policy $\pi^*$.
Given known functions $\psi_j:\mathcal{S}^j\times \mathcal{T}^j\rightarrow\mathbb{R}^{d'}$ with $d'\geq d$, write
the instruments $\Psi_{i,j}:=\psi_j(S_{i,1:j}, T_{i,1:j})$, and define $\bar{\Psi}_{i,j}:=\mathbb{E}^+_{i,j}[\Psi_{i,j}]$.
Under Assumptions~\ref{assump:exogeneity}~\&~\ref{assump:linear_blip_function}, the true parameter vector $\theta^*$ satisfies:
\begin{equation}
\label{eq:conditional_moment_identification}
\mathbb{E}_{i,j}^+\left[
R_{i,j}
\left( \Psi_{i,j} - \bar{\Psi}_{i,j} \right)
\right]= 0, \quad \forall i\in\{1:n\}, j\in\{0:L\},
\end{equation}
with the convention that $\Psi_{i,0}=1$, $\bar{\Psi}_{i,0}=0$, and $\Phi_{i,0}=1$.
\end{lemma}
\begin{remark}
One approach to constructing instruments is to set $\Psi_{i,j} = \Phi_{i,j}$ by defining
\[\psi_j(s_{1:j}, \tau_{1:j})=\phi_j(s_{1:j}, \tau_{1:j}) - \phi_j(s_{1:j}, \tau_{1:j-1},\pi^*(s_{1:j}, \tau_{1:j-1}))\]
for any $(s_{1:j},\tau_{1:j})\in\mathcal{S}^j\times \mathcal{T}^j$. When additional information is available, it is possible to include more instruments (such that $\Psi_{i,j}$ has higher dimension than $\Phi_{i,j}$) to improve estimation efficiency.
\end{remark}
The moment conditions in the above lemma can be expressed more compactly as a single vector of moment constraints. Given a unit $i$ with episodic data $Z_i$,
define $\beta_i:=(\beta_{i,0},\beta_{i,1},\dots, \beta_{i,L})$, where each component is given by $\beta_{i,j}=Y_i(\Psi_{i,j} -\bar \Psi_{i,j})$; define
\begin{align*}
J_i =- \begin{bmatrix}
(\Psi_{i,0}-\bar\Psi_{i,0})\Phi_{i,0}^\top
&\quad(\Psi_{i,0}-\bar\Psi_{i,0})\Phi_{i,1}^\top &\quad\dots &\quad(\Psi_{i,0}-\bar\Psi_{i,0})\Phi_{i,L}^\top \\
\mathbf{0} &\quad(\Psi_{i,1}-\bar\Psi_{i,1}) \Phi^\top_{i,1} &\quad \dots &\quad(\Psi_{i,0}-\bar\Psi_{i,0})\Phi^\top_{i,L}\\
\vdots &\quad \ddots &\quad\ddots &\quad\vdots\\
\mathbf{0} &\quad \dots &\quad \mathbf{0} &\quad (\Psi_{i,L}-\bar\Psi_{i,L})\Phi^\top_{i,L}
\end{bmatrix}.
\end{align*}
Then
for any parameter vector $\theta$, the moment function $m(Z_i;\theta)$ can be written as:
\begin{equation}
m(Z_i;\theta) =~ \beta_i + J_i \theta.
\label{eq:moment}
\end{equation}
In particular, the moment condition in \eqref{eq:moment} can be constructed directly from observed data. The key quantity, $\bar\Psi_{i,j} = \mathbb{E}_{i,j}^+[\Psi_{i,j}]$, is the conditional expectation over the treatment assignment $T_{i,j}$ under the behavior policy. Since the behavior policy is assumed to be known, this expectation can be computed exactly. As a result, the moment function $m(Z_i; \theta)$ is also computable from the observed data.
Finally, Lemma \ref{lemma:identification_parameter} implies that the true model parameter vector $\theta^*$ solves the expected moment equation:
\begin{align}
\label{eq:moment_condition}
\bar\beta_i+\bar J_i \theta^* =\mathbf{0}, \quad \forall i = 1,\dots, n,
\end{align}
where $\bar\beta_i:=(\bar\beta_{i,0},\bar\beta_{i,1},\dots, \bar\beta_{i,L})$ for $\bar\beta_{i,j}=\mathbb{E}_{i,j}^+\left[Y_i(\Psi_{i,j} -\bar \Psi_{i,j})\right]$, and
\begin{align*}
\bar{J}_i :=- \begin{bmatrix}
\mathbb{E}_{i,0}^+[(\Psi_{i,0}-\bar{\Psi}_{i,0})\Phi_{i,0}^\top]
&\quad
\mathbb{E}_{i,0}^+[(\Psi_{i,0}-\bar{\Psi}_{i,0})\Phi_{i,1}^\top]
&\quad
\dots
&\quad
\mathbb{E}_{i,0}^+[(\Psi_{i,0}-\bar{\Psi}_{i,0})\Phi_{i,L}^\top] \\
\mathbf{0} &\quad
\mathbb{E}_{i,1}^+[(\Psi_{i,1}-\bar{\Psi}_{i,1})\Phi_{i,1}^\top] &\quad
\dots &\quad
\mathbb{E}_{i,1}^+[(\Psi_{i,1}-\bar{\Psi}_{i,1})\Phi_{i,L}^\top]\\
\vdots &\quad
\ddots &\quad \ddots &\quad \vdots\\
\mathbf{0} &\quad \dots &\quad \mathbf{0} &\quad \mathbb{E}_{i,L}^+[(\Psi_{i,L}-\bar{\Psi}_{i,L})\Phi_{i,L}^\top]
\end{bmatrix}.
\end{align*}
This observation will become handy in our subsequent analysis.
\subsection{Adaptively Weighted Generalized Method of Moments (AW-GMM)}
Given the moment conditions in \eqref{eq:moment_condition}, one can estimate the target parameter $\theta^*$ using moment-based methods such as Z-estimation or the generalized method of moments (GMM). To improve stability and efficiency (particularly when $\Psi$ has higher dimension than $\Phi$), we focus on GMM throughout this paper. The results also apply directly to Z-estimation. Specifically, the standard GMM estimator $\tilde{\theta}_n$ minimizes the norm of the empirical moment equations:
\begin{equation*}
\tilde{\theta}_n \in \argmin_{\theta\in\Theta}~\left\|\frac{1}{n}\sum_{i=1}^n m (Z_i;\theta)\right\|_A^2 =
\left\{\frac{1}{n}\sum_{i=1}^n m (Z_i;\theta)\right\}^\top A \left\{\frac{1}{n}\sum_{i=1}^n m (Z_i;\theta)\right\},
\end{equation*}
where $A$ is a positive-definite matrix and $\|\cdot\|_A$ denotes the norm induced by $A$. When $\tilde{\theta}_n$ lies in the interior of $\Theta$, it satisfies the first-order condition: $\{\frac{1}{n}\sum_{i=1}^n J_i\}A\{\frac{1}{n}\sum_{i=1}^n m (Z_i;\tilde\theta_n)\}=\mathbf{0}$.
To build more intuition, consider the case where the sample average Jacobian $\frac{1}{n}\sum_{i=1}^n J_i$ is invertible. Then the first-order condition implies that the GMM estimator solves the empirical moment equations such that: $\frac{1}{n}\sum_{i=1}^n m (Z_i;\tilde\theta_n)=\mathbf{0}$. The asymptotic behavior of $\tilde{\theta}_n$ is typically analyzed through a Taylor expansion around the true parameter $\theta^*$, which links the normalized estimation error $J_i(\tilde{\theta}_n-\theta^*)$ with an empirical influence function $m(Z_i; \theta^*)$:
\begin{align}
\frac{1}{n} \sum_{i=1}^n J_i (\tilde{\theta}_n - \theta^*) = \frac{1}{n} \sum_{i=1}^n m(Z_i; \tilde\theta_n) -\frac{1}{n} \sum_{i=1}^n m(Z_i; \theta^*) =-\frac{1}{n} \sum_{i=1}^n m(Z_i; \theta^*)
\end{align}
However, since the behavior policy is adaptively evolving over time, the Jacobian $J_i$ is likely to diverge,
and the variance of the empirical influence functions on the right-hand-side might fluctuate and never converge.
This risk of instability motivates us to modify the GMM estimator by applying \emph{non-uniform} and \emph{time-varying} adaptive weights, which are carefully chosen to stabilize the variance of the empirical influence functions $m(Z_i;\theta^*)$ over time,
thereby preserving the desirable asymptotic characteristics of the GMM-estimator.
We hereby propose our \emph{adaptively weighted} generalized method of moments (AW-GMM) estimator $\hat{\theta}_n$, which minimizes the norm of a \emph{non-uniform} average of empirical moment equations:
\begin{align}
\label{eq:gmm_estimator_linear}
\hat{\theta}_n \in \argmin_{\theta\in\Theta}~\mathcal{L}_n(\theta):=\left\|\frac{1}{n}\sum_{i=1}^n H_i\,m (Z_i;\theta)\right\|_A^2,
\end{align}
where $H_{i}$ denotes the time-varying weighting matrix designed to counterbalance the potential divergence in the empirical influence functions. Specifically, we consider a block-diagonal structure for $H_{i}$: $H_i:= \mbox{diag}\left\{H_{i,0}, H_{i,1}, \dots, H_{i,L}\right\}$, where $H_{i,0}\in\mathbb{R}$ and $\{H_{i,j}\in \mathbb{R}^{d'\times d'}$ for each $j\in [1:L]\}$; each diagonal block $H_{i,j}$ is adapted to filtration $\mathcal{F}_{i,j}^+$ and
stabilizes the
empirical influence functions
of data $Z_i$ involved in estimating the structural parameter $\theta_j^*$.
Similar weighting schemes have been introduced in the bandit setups \citep{deshpande2018accurate,hadad2021confidence,zhang2021statistical,bibaut2021post,zhan2021off}, though none applies to RL data with more than one stages.
The technical challenge in identifying proper weights for RL data is two-fold.
First, it is crucial that these weights remain stable and do not become degenerate or explode, otherwise the resulting estimator would suffer from diverging variance. This is known as the ``curse of horizon'' challenge, which arises from the deteriorating overlap between the behavior and evaluation policies over lengthy episodes. Previous RL literature on weighting-based offline evaluation has often turned to heuristic solutions, such as weight clipping, to manage variance, albeit at the expense of introducing bias \citep{precup2000eligibility,chen2019top}. We show that by carefully designing our weighting matrix $H_i$, our estimator not only controls estimation variance but also achieves asymptotic unbiasedness.
Second, we want to construct these weights using information only up to current observations, following the practice of constructing weights for adaptive bandit data \citep{hadad2021confidence}. By doing so, we can employ martingale limiting theorems to show
that the weighted GMM-estimator admits asymptotic normality and achieves post-RL inference.
\section{Consistency of AW-GMM Estimation}
\label{sec:consistency}
We begin by introducing a general weighting scheme that ensures consistent estimation of target parameters from adaptively collected RL data. We then apply this framework to specific RL algorithms later in the section. Recall that the weighted GMM estimator $\hat{\theta}_n$ minimizes the empirical loss given by:
\begin{equation}
\label{eq:gmm_loss}
\mathcal{L}_n(\theta):=\left\|\frac{1}{n}\sum_{i=1}^n H_im (Z_i;\theta)\right\|_A^2 =\left\|\frac{1}{n}\sum_{i=1}^nH_i(\beta_i+J_i\theta)\right\|_A^2.
\end{equation}
By the conditional moment equation \eqref{eq:moment_condition} and weight $H_{i,j}$ adapting to $\mathcal{F}_{i,j}^+$, we have the true $\theta^*$ satisfies $H_i(\bar\beta_i+\bar J_i\theta^*)=\mathbf{0}$ and thus is a minimizer of the norm of averaged conditional moments:
\begin{equation}
\label{eq:conditional_moment}
\mathcal{L}(\theta):=\left\|\frac{1}{n}\sum_{i=1}^nH_i(\bar\beta_i+\bar J_i\theta)\right\|_A^2,
\end{equation}
with $\mathcal{L}(\theta^*)=0$.
To establish that $\hat{\theta}_n$ is consistent for $\theta^*$, it suffices to show two conditions: (i) for every $\theta\in\Theta$, the GMM loss $ \mathcal{L}_n(\theta)$ converges to the expected loss $ \mathcal{L}(\theta)$; and (ii) the loss $\mathcal{L}(\theta)$ has a unique minimizer. The following weighting property ensures both conditions are satisfied.
\begin{property}[$(\alpha_1, \alpha_2)$-Regularizing Weights]
\label{property:weight_regularizing}
We set $H_{i,0}=1$ for $i\in[1:n]$. Given a matrix $A$, we use $\lambda_{\min}(A)$ to denote its smallest eigenvalue and $\mathrm{Tr}(A)$ to denote its trace (sum of diagonal entries).
For each $j\in[1:L]$, the weights $\{H_{i,j}\}_{i=1}^n$ are adapted to filtration $\{\mathcal{F}_{i,j}^+\}_{i=1}^n$ and satisfy that:
\begin{enumerate}[(a)]
\item $\lambda_{\min}\left(\left\{\frac{1}{n}\sum_{i=1}^n \mathbb{E}_{i,j}\left[H_{i,j} \mathrm{Cov}^+_{i,j}(\Psi_{i,j},\Phi_{i,j})\right]\right\}^\top \left\{\frac{1}{n}\sum_{i=1}^n \mathbb{E}_{i,j}\left[H_{i,j} \mathrm{Cov}^+_{i,j}(\Psi_{i,j},\Phi_{i,j})\right]\right\}
\right) \geq ~ c_1^2\cdot n^{-\alpha_1}$;
\item $ \frac{1}{n}\sum_{i=1}^n \mathrm{Tr}\left(H_{i,j}H_{i,j}^\top\right)\leq c_1^{-1}\cdot n^{\alpha_1}$;
\item $ \frac{1}{n} \sum_{i=1}^n \mathrm{Tr}\left(\mathbb{E}[H_{i,j} \mathrm{Var}_{i,j}^+(\Psi_{i,j}) H_{i,j}^\top]\right) \leq c_1^{-1}n^{\alpha_2}$.
\end{enumerate}
for universal constants $c_1>0$, $\alpha_1\in[0,\frac{1}{L+1})$, and $\alpha_2 \in[0, \frac{1}{L-1}]$.\footnote{Boundedness of $\Psi$ yields $\mathrm{Tr}(H_{i,j}\mathrm{Var}_{i,j}^+(\Psi_{i,j})H_{i,j}) \lesssim \mathrm{Tr}(H_{i,j}H_{i,j})$. As a result, we have $ \frac{1}{n} \sum_{i=1}^n \mathrm{Tr}\left(\mathbb{E}[H_{i,j} \mathrm{Var}_{i,j}^+(\Psi_{i,j}) H_{i,j}^\top]\right)\lesssim\frac{1}{n}\sum_{i=1}^n \mathrm{Tr}\left(H_{i,j}H_{i,j}^\top\right)$. To make Property~\ref{property:weight_regularizing}(c) nontrivial at the presence of Property~\ref{property:weight_regularizing}(b), we consider $\alpha_2\leq \alpha_1$.}
\end{property}
In particular, Properties \ref{property:weight_regularizing}(a) \& \ref{property:weight_regularizing}(b) ensure that $ \mathcal{L}(\theta)$ has a unique minimizer with high probability, and Property \ref{property:weight_regularizing}(c) ensures that $\mathcal{L}_n(\theta)$ uniformly converges to $ \mathcal{L}(\theta)$ for any bounded $\theta\in\Theta$.
\begin{theorem}
\label{thm:consistency}
Suppose Assumptions~\ref{assump:exogeneity}~\&~\ref{assump:linear_blip_function}
hold and the weights $\{H_i\}$ satisfy the $(\alpha_1, \alpha_2)$-regularizing Property \ref{property:weight_regularizing}. Also assume that $\mathcal{S},\mathcal{T},\mathcal{Y}$ are bounded. The AM-GMM estimator $\hat{\theta}_n$ specified in \eqref{eq:gmm_estimator_linear} achieves the following consistency result for $j\in[0:L]$:
\begin{equation}
\mathbb{E}\left[\left\|\hat{\theta}_{n,j}-\theta_{j}^*\right\|_2\right]=O\left(n^{\frac{(L-j+1)(\alpha_1+\alpha_2)-1}{2}}\right),\quad\mbox{and}\quad \mathbb{E}\left[\left\|\hat{\theta}_{n,j}-\theta_{j}^*\right\|^2_2\right] =O\left(n^{(2L-2j+1)\alpha_1+\alpha_2-1}\right).
\end{equation}
\end{theorem}
Consider that the data is collected by an RL agent with a fixed amount of exploration, such as the $\epsilon$-greedy RL algorithms with a constant $\epsilon$. Then the uniform weighting with $H_i$ being the identity matrix satisfies Property \ref{property:weight_regularizing} with $\alpha_1=\alpha_2=0$, and Theorem \ref{thm:consistency} recovers the $n^{-1/2}$ convergence rate, consistent with results in \citep{lewis2020double} for i.i.d.~data. However, our result is stronger in the sense that we allow for adaptivity in the experiments, as opposed to the i.i.d.~settings in \citep{lewis2020double}.
\subsection{Consistency Weights for Decaying Exploration under Bilinear Features}
\label{sec:consistency_application}
We now provide explicit weighting schemes for common RL algorithms with polynomially decaying exploration rates.
We particularly focus on scenarios where the instrument $\psi$ adopts the feature map $\phi$ such that $\Psi_{i,j}=\Phi_{i,j}$, and the feature map $\phi$ can be expressed as a bilinear form of $\phi_j(S_{i,1:j}, T_{i,1:j}) = \chi_{j}(S_{i,1:j}, T_{i,1:j-1}) \otimes \mu(T_{i,j})$,
with $\otimes$ being the Kronecker product.\footnote{We remind that for two vectors $v\in R^p$ and $u\in R^d$, the Kronecker product is the vector whose entries contain the product of all pairs of entries of the two vectors, i.e., $v\otimes u = \text{vec}( uv^\top) = (v_1 u_1, \ldots, v_1 u_d, \ldots, v_p u_1, \ldots, v_p u_d)$.} Such feature maps accommodate a wide family of categorical treatments and continuous treatments (we defer the examples to the end of this section).
Without loss of generality, we rename $\mu(T_{i,j})$, the transformation of treatment, as $T_{i,j}$ and focus on the below form throughout:
\begin{assumption}[Bilinear Feature Map]\label{ass:bilinear} The feature map $\phi_j$ of the blip function takes the form:
\begin{equation}
\label{eq:feature_mapping_bilinear}
\phi_j(S_{i,1:j}, T_{i,1:j}) =(T_{i,j} - \pi^*(S_{i,1:j}, T_{i,1:j-1}))\otimes \chi_{j}(S_{i,1:j}, T_{i,1:j-1}).
\end{equation}
\end{assumption}
Let $d_\tau$ be the dimension of treatment $T_{i,j}$, and we use $X_{i,j}$ as a shorthand for $\chi_j(S_{i,1:j}, T_{i,1:j-1})$. Thus we can also write that:\footnote{Here we leverage the property that for two vectors $u\in \mathbb{R}^p$ and $v\in \mathbb{R}^d$, $v\otimes u=(I_d \otimes u) \cdot v$, where $I_d \otimes u = \text{diagonal}(u, \ldots, u)$. Hence, our feature map can equivalently be written as $\Phi_{i,j} = (I_{d_{\tau}} \otimes \Psi_{i,j}) \cdot T_{i,j}$.}
\begin{align*}
\Phi_{i,j} =
(T_{i,j} - \pi^*(S_{i,1:j}, T_{i,1:j-1}))\otimes X_{i,j} = (I_{d_{\tau}} \otimes X_{i,j}) \cdot (T_{i,j} - \pi^*(S_{i,1:j}, T_{i,1:j-1}))
\end{align*}
To identify the structure parameter $\theta^*$, one should expect sufficient overlap condition in the behavior policy and the state transition dynamics, due to co-linearities in the linear system that identifies the structural parameters. The following assumption formalizes the overlap condition.
\begin{assumption}[Overlap]
\label{assump:overlap}
Let $c>0$ be some universal constant. For each unit $i$ at stage $j$,
\begin{enumerate}[(a)]
\item The behavior policy satisfies that,
\begin{equation}
\label{eq:require_observational_policy}
\mathrm{Var}^+_{i,j}(T_{i,j})\succeq c_j\cdot i^{-\alpha}\cdot I_{d_\tau},
\end{equation}
for constants $c_j>0$ and $\alpha\in[0,1)$. We refer to $\alpha$ as behavior exploration rate.
\item The state transition satisfies that $\|X_{i,j}\|^2_2\in[c,c^{-1}]$ and $\mathbb{E}_{i,j}[X_{i,j}X_{i,j}^\top]\succeq c\cdot I$ almost surely.
\end{enumerate}
\end{assumption}
Assumption \ref{assump:overlap}(a) imposes an explicit rate on the decay of the amount of randomization (exploration), which is employed by the behavior policy of the RL agent as a function of the number of observations it has collected so far. Similar assumption is also required in ex post inference when using bandit data~\citep{hadad2021confidence,zhan2021off}. In particular, when $\alpha=0$, this becomes a \emph{relaxed} version (by allowing for the dependence among observations) of the commonly made ``overlap'' condition in the causal inference literature on i.i.d.~samples \citep{imbens2004nonparametric}.
Assumption \ref{assump:overlap}(b) says that the state transition dynamic, which is independent of the behavior policy, is non-degenerate. Specifically, conditional on the past states and actions $(S_{i,1:j-1}, T_{i,1:j-1})$, the next state feature $X_{i,j}$ is not collinear; it can be achieved if the exogenous noise injected at each stage has a full-rank covariance matrix.
With these, we instantiate the weighting choices that satisfy Property \ref{property:weight_regularizing}, yielding the corollary of Theorem \ref{thm:consistency}.
\begin{corollary}
\label{cor:consistency}
Suppose Assumptions~\ref{assump:exogeneity},~\ref{assump:linear_blip_function},~\ref{ass:bilinear},~\&~\ref{assump:overlap} hold.
\begin{enumerate}
\item If $H_{i,j}$ is the identity matrix (i.e. no weighting), then Property \ref{property:weight_regularizing} is satisfied with $\alpha_1=2\alpha$, $\alpha_2=0$, and we have:
\begin{equation*}
\mathbb{E}\left[\left\|\hat{\theta}_{n}-\theta^*\right\|_2\right] =O\left(n^{\frac{2(L-j+1)\alpha-1}{2}}
\right) .
\end{equation*}
\item The weighting scheme $H_{i,j} = (I_{d_\tau}\otimes X_{i,j})\mathrm{Var}^+_{i,j}(T_{i,j})^{-1/2}
(I_{d_\tau}\otimes X_{i,j})^\top $ satisfies Property \ref{property:weight_regularizing} with $\alpha_1=\alpha$ and $\alpha_2=0$, and
we have:
\begin{equation*}
\mathbb{E}\left[\left\|\hat{\theta}_{n}-\theta^*\right\|_2\right] =O\left(n^{\frac{(L-j+1)\alpha-1}{2}}
\right).
\end{equation*}
\end{enumerate}
\end{corollary}
\begin{remark}
Corollary~\ref{cor:consistency} highlights the potential improvements enabled by a weighted GMM estimator. The weighting scheme in Corollary~\ref{cor:consistency}(2) depends on $\mathrm{Var}^+_{i,j}(T_{i,j})$, which can be computed exactly given the known behavior policy. Incorporating these weights ensures consistency even under faster rates of exploration decay. For example, in bandit settings, uniform weighting restricts the exploration rate $\alpha$ (as defined in the lower bound on treatment variance in \eqref{eq:require_observational_policy}) to be at most $\frac{1}{2}$, whereas adaptive weighting allows $\alpha$ to reach $1$.
\end{remark}
We conclude this section by instantiating common feature maps and RL algorithms that satisfy \eqref{eq:feature_mapping_bilinear}.
\begin{example}[Categorical Treatment] Consider the action space of the RL agent contains one control action and $K$ treatment action. Omitting the episode index $i$ for simplicity, let $T_{j}=\mathbf{0}$ denote the control action being taken, and $T_{j}=e_k$ denote the $k$-th treatment being taken, at the $j$-th period of unit $i$, where $e_k$ is the one-hot vector for category $k$. Without loss of generality, the blip function can be represented as
\begin{align}
\label{eq:categorical_treatment}
(\theta_j^*)^\top\phi_j(S_{1:j}, T_{1:j}) = \sum_{k=1}^K (\theta_{j,k}^*)^\top \chi_{j}(S_{1:j}, T_{1:j-1}) \mathbf{1}[T_{j}=e_k],
\end{align}
where $(\theta_{j,k}^*)^\top \chi_{j}(S_{1:j}, T_{1:j-1})$ characterizes the heterogeneous treatment effect of the $k$-th treatment at stage~$j$. Note that \eqref{eq:categorical_treatment} is equivalent to \eqref{eq:feature_mapping_bilinear} by expanding the Kronecker product in \eqref{eq:feature_mapping_bilinear}, i.e.
\begin{align*}
(\theta_j^*)^\top\phi_j(S_{1:j}, T_{1:j}) = (\theta_{j}^*)^\top (
T_{j}\otimes \chi_{j}(S_{1:j}, T_{1:j-1}) ).
\end{align*}
To meet the overlap condition on behavior policy in Assumption \ref{assump:overlap}, it suffices to run an $\epsilon$-greedy RL algorithm, where $\epsilon$ decays over time. At each decision time $t$, the agent chooses the treatment that was best performing on historical data for that stage with probability $1-\epsilon_t$ with $ \epsilon_t= t^{-\alpha}$; alternatively, the agent opts for a random treatment.
\end{example}
\begin{example}[Continuous or Binary Treatment Vector] Consider the case when the treatment is a continuous (or binary) vector in $\mathbb{R}^{d_\tau}$, with the blip function represented as
\begin{align}
\label{eq:continuous_treatment}
(\theta_j^*)^\top \phi_j(S_{1:j}, T_{1:j}) = \sum_{k=1}^{d_\tau} (\theta_{j,k}^*)^\top \chi_{j}(S_{1:j}, T_{1:j-1}) T_{j,k}.
\end{align}
Each coordinate can be viewed as a separate continuous treatment applied to the unit, where different treatments can be applied simultaneously to each unit.
Here, $(\theta_{j,k}^*)^\top \chi_{j}(S_{1:j}, T_{1:j-1})$ characterizes the heterogeneous marginal treatment effect from the $k$-th coordinate of the treatment vector $T_{j}$, equivalently, the $k$-th treatment applied to the unit. Similarly, Equation~\eqref{eq:continuous_treatment} falls in the functional form prescribed in Equation~\eqref{eq:feature_mapping_bilinear}.
To satisfy the overlap of behavior policy required in Assumption \ref{assump:overlap}, it suffices to add an exogenous exploration noise $\varepsilon_t\in \mathbb{R}^{d_\tau}$ to treatment assignment at time $t$, where $\varepsilon_t$ has covariance matrix $t^{-\alpha}I_{d_\tau}$. If each of the treatments are binary, then it suffices to independently randomize the assignment of each of the simultaneous binary treatments.
\end{example}
\begin{example}[Polynomial Scalar Treatment] Consider the case when a single scalar treatment $\tau_{i,j}$ can be applied at stage $j$ for unit $i$. We can express non-linear effects of the scalar treatment by considering a fixed expansion to a set of engineered treatment features $\mu(\tau_{i,j})$, as follows:
\begin{align}
\label{eq:polynomial_treatment}
(\theta_j^*)^\top \phi_j(S_{1:j}, T_{1:j}) = \sum_{k=1}^{d_\tau} (\theta_{j,k}^*)^\top \chi_{j}(S_{1:j}, T_{1:j-1}) \mu_k(\tau_{i.j})
\end{align}
which is equivalent to \eqref{eq:feature_mapping_bilinear} by setting $T_{i,j}$ in \eqref{eq:feature_mapping_bilinear} to $\mu(\tau_{i,j})$. For instance, $\mu(\tau_{i,j})$ could be chosen to be a high-degree polynomial, i.e. $\mu(\tau_{i,j}) = (\tau_{i,j}, \tau_{i,j}^2, \tau_{i,j}^3, \ldots, \tau_{i,j}^{d_\tau})$. In the context of a pricing application, one can view $\tau_{i,j}$ as the price offering for some product to a buyer $i$, and the outcome as revenue (i.e. whether there was a purchase times the purchase price). Aggregate revenue would typically be some bell-shaped curve, which can be well approximated by a third or fourth degree polynomial of price (a typical choice in empirical work). In such a revenue model, the quantities $(\theta_{j,k}^*)^\top \chi(S_{1:j}, T_{1:j-1})$ correspond to the heterogeneous coefficients in this parametric revenue model, reflecting heterogeneous price elasticities based on the current state of the buyer.
\end{example}
\section{Asymptotic Normality of AW-GMM Estimation}
\label{sec:normality}
In this section, we identify a family of weighting schemes, under which the AW-GMM estimator is asymptotically normal. These weights perfectly stabilize the variance of the empirical influence function in the asymptotic regime. Building on this, we prove strong Gaussian approximation results and characterize the uniform convergence rate.
These results are practically appealing, as they enable the construction of uniformly valid confidence regions for parameters of interest over a large class of RL algorithms.
Recall that the AW-GMM estimator $\hat{\theta}_n$ minimizes the empirical loss $\mathcal{L}_n(\theta)$ while the true parameter $\theta^*$ minimizes the expected loss $\mathcal{L}(\theta)$:
\begin{equation*}
\hat\theta_n\in\argmin_{\theta\in\Theta} \mathcal{L}_n(\theta) = \left\|\frac{1}{n}\sum_{i=1}^n H_i\, \left(\beta_i+J_i\theta\right)\right\|_A^2\quad \mbox{and}\quad \theta^*\in\argmin_{\theta\in\Theta} \mathcal{L}(\theta)=\left\|\frac{1}{n}\sum_{i=1}^n H_i\left( \bar{\beta}_i+\bar{J}_i\theta\right)\right\|_A^2.
\end{equation*}
The key idea is to relate the estimation error $(\hat{\theta}_n - \theta^*)$ to a martingale difference sequence (MDS), which allows us to apply martingale theory to characterize the asymptotic distribution of $\hat{\theta}_n$.
As established in Section~\ref{sec:consistency}, the estimator $\hat{\theta}_n$ is consistent for $\theta^*$. Since $\theta^*$ lies in the interior of $\Theta$, $\hat{\theta}_n$ will also lie in the interior with high probability. In this case, $\hat{\theta}_n$ satisfies the first-order condition of the empirical loss $\mathcal{L}_n(\theta)$:
\begin{equation}
\nabla_\theta \mathcal{L}_n(\hat\theta_n) =\left(\frac{1}{n}\sum_{i=1}^n H_i J_i\right)^\top A\left(\frac{1}{n}\sum_{i=1}^n H_i\, (\beta_i+J_i\hat\theta_n)\right)=\mathbf{0}.
\end{equation}
Define $\xi_i := \beta_i+J_i\theta^* $. Then after some algebra, we have:
\begin{align}
\label{eq:link_error_mds}
-\left(\frac{1}{n}\sum_{i=1}^n H_i J_i\right)^\top A\left(\frac{1}{n}\sum_{i=1}^n H_i \xi_i\right)
=& \left(\frac{1}{n}\sum_{i=1}^n H_i J_i\right)^\top A\left(\frac{1}{n}\sum_{i=1}^n H_iJ_i(\hat\theta_n-\theta^*)\right).
\end{align}
Since $\theta^*$ satisfies the conditional moment equations $\bar{\beta}_i+\bar J_i\theta^*=\mathbf{0}$, the sequence $\{H_i \xi_i\}$ forms an MDS under suitable conditions on the weights $H_i$, as formally established in Section~\ref{sec:mds}. By linking $\{H_i \xi_i\}$ to the estimation error $(\hat{\theta}_n - \theta^*)$, we derive a strong Gaussian approximation for $\hat{\theta}_n$ and develop corresponding inference results in Section~\ref{sec:strong_gaussian}. We then apply these results to common RL algorithms in Section~\ref{sec:normality_application}.
\subsection{Martingale from Adaptive RL Data}
\label{sec:mds}
We now formalize the martingale difference sequence structure, which plays a central role in our analysis.
The true parameter $\theta^*$ satisfies the conditional moment equations, so that $ \bar \beta_i + \bar J_i\theta^* =\mathbf{0}$. Define
\begin{equation}
\xi_i := \beta_i +J_i\theta^* = (\xi_{i,0}, \xi_{i,1}, \dots, \xi_{i,L}),\quad \mbox{with}\quad \xi_{i,j}=R_{i,j}(\Psi_{i,j} - \bar{\Psi}_{i,j}), \forall j\in[0:L],
\end{equation}
where $R_{i,j}$, defined in Lemma~\ref{lemma:identification_parameter},
captures the residual component of $Y_i$ after removing the treatment effects from $T_{i,j:L}$ and adding the effects under the evaluation policy $\pi^*_{j:L}$. The sequence $\{\xi_i\}$ forms a martingale difference sequence, as established in the lemma below.
\begin{lemma}
\label{lemma:mds-1}
Suppose Assumptions~\ref{assump:exogeneity} and~\ref{assump:linear_blip_function} hold.
Let $H_i$ be a block-diagonal adaptive weighting matrix, where the $j$-th diagonal block $H_{i,j}$ is measurable with respect to $\mathcal{F}^+_{i,j}$.
Then the sequence $\{H_i\xi_i\}_{i=1}^n$ forms a martingale difference sequence adapted to the filtration $\{\mathcal{F}_i\}_{i=1}^n$.
\end{lemma}
Equation~\eqref{eq:link_error_mds} shows that the estimation error $(\hat{\theta}_n-\theta^*)$ is closely linked to the sum $\sum_{i=1}^nH_i\xi_{i}$.
To analyze the asymptotic properties of $\hat{\theta}_n$, it thus suffices to understand the MDS, $\{H_i\xi_i\}$. The following result highlights the key components that characterize its asymptotic behavior.
\begin{proposition}[Martingale CLT, \cite{hall2014martingale}]
\label{prop:martingale_clt_2}
Let $\{\phi_{t}, \mathcal{F}_{t}\}_{t=1}^\top$, with $\phi_t\in \mathbb{R}$, be a square-integrable scalar martingale difference sequence. Suppose that the two conditions below are satisfied,
\begin{itemize}
\item[(a)] conditional variance convergence: $\sum_{t=1}^\top \mathbb{E}[\phi_{t}^2|\mathcal{F}_{t-1}]\xrightarrow{p} \eta^2$ for some a.s. finite r.v. $\eta^2$;
\item[(b)] conditional higher-moment decay: $\sum_{t=1}^\top \mathbb{E}[\phi_{t}^4|\mathcal{F}_{t-1}]\xrightarrow{p}0$.
\end{itemize}
Then, $\sum_{t=1}^\top \phi_{t}\Rightarrow Z$, where the random variable Z has characteristic function $\mathbb{E}[\exp(-\frac{1}{2}\eta^2t^2)]$.
\end{proposition}
Proposition~\ref{prop:martingale_clt_2} illuminates how weights $H_i$ can be designed to improve the asymptotic properties of the MDS sum $\sum_{i=1}^n H_i\xi_i$, and thus our AW-GMM estimation. A key quantity is the sum of conditional variances, $\sum_{i=1}^n \mathrm{Var}_{i}(H_i\xi_{i})$: when it converges to the identity matrix, this MDS sum is asymptotically normal.
\subsubsection{Variance and Covariance under Homoscedasticity.}
We now take a closer look at the conditional variance $\mathrm{Var}_{i}(H_i\xi_{i})$, which forms the basis for our subsequent weighting strategy.
In full generality, one could design the weights $H_i$ to approximate $\mathrm{Var}_i(\xi_i)^{-1/2}$, so that $\mathrm{Var}_i(H_i \xi_i)$ is close to the identity matrix. However, this approach would require estimating $\mathrm{Var}_i(\xi_i)$, which in turn depends on the state transition dynamics and is thus an impractical task, particularly in high-dimensional continuous state spaces.
To address this challenge, we introduce a homoscedasticity assumption on the outcome process. Under this assumption, the conditional variance matrix $\mathrm{Var}_i(\xi_i)$ becomes block-diagonal: the covariances across different stages (corresponding to the off-diagonal blocks) are zero. As a result, it suffices to design weights conditioning on the realized state $S_{i,j}$ to stabilize the stage-specific conditional variances of $\xi_{i,j}$, thereby eliminating the need to estimate the state transition dynamics.
\begin{assumption}[Homoscedasticity of residuals with respect to current treatment]
\label{assump:homoscedasticity}
The conditional variance of the residuals $\mathrm{Var}(R_{i,j}\mid T_{i,j}, \mathcal{F}_{i,j}^+)$ is independent of $T_{i,j}$, i.e. $\mathrm{Var}(R_{i,j}\mid T_{i,j}, \mathcal{F}_{i,j}^+)=\mathrm{Var}_{i,j}^+(R_j)$.
Moreover, for the residual $R_{i,j}$ defined in Lemma \ref{lemma:identification_parameter}, there exists universal constants $\sigma^2, M^2$ such that $R_{i,j}^2\leq M^2$ and $\mathrm{Var}_{i,j}^+(R_{i,j}) \geq \sigma^2$ almost surely.
\end{assumption}
As shown in Appendix~\ref{appendix:proof_mds}, Assumption~\ref{assump:homoscedasticity} is equivalent to $\mathbb{E}_{i,j}^+[R_{i,j}^2 (\Psi_{i,j} - \bar{\Psi}_{i,j})]=0$, while
Lemma~\ref{lemma:identification_parameter} states that $\mathbb{E}_{i,j}^+[R_{i,j} (\Psi_{i,j}-\bar{\Psi}_{i,j})]=0$, reflecting that residuals are uncorrelated with the current treatment conditional on past information.
Assumption~\ref{assump:homoscedasticity} strengthens this by requiring the same uncorrelated property for the residual variance. It is satisfied under the additive, rank-preserving version of the structural nested mean model \cite[Chapter 14.5]{miguel2023causal}, but is more general.
It implies that, conditional on the past, variance of $R_{i,j}$ is independent of future treatment assignments and is fully explained by exogenous factors.
Assumption~\ref{assump:homoscedasticity} holds, for instance, if $R_{i,j}\, {\perp\!\!\!\perp}\, T_{i,j}\mid \mathcal{F}^+_{i,j}$, which is natural if $R_{i,j}$ equals to the potential outcome $\ensuremath{Y}(T_{i,1:j-1}, \pi^*_{j:L})$ in distribution given the past.\footnote{This equality is stronger than Lemma~\ref{lemma:identification_policy_value}, which requires only equality in expectation.}
Given Assumption~\ref{assump:homoscedasticity}, we layout the block diagonal structure of the conditional variance $\mathrm{Var}_i(H_i\xi_i)$.
\begin{lemma}
\label{lemma:mds}
Suppose Assumptions~\ref{assump:exogeneity}, \ref{assump:linear_blip_function}~\&~\ref{assump:homoscedasticity} hold.
Let $H_i$ be a block-diagonal adaptive weighting matrix, where the $j$-th diagonal block $H_{i,j}$ is measurable with respect to $\mathcal{F}^+_{i,j}$.
Then the conditional variance matrix $\mathrm{Var}_i(H_i\xi_i)=\mathbb{E}_i[H_i\xi_i\xi_i^\top H_i^\top]$ is block-diagonal, such that:
\begin{enumerate}
\item[(a)] its $(j_1, j_2)$ off-diagonal block
$\mathrm{Cov}_i(H_{i,j_1}\xi_{i,j_1}, \xi_{i,j_2}H_{i,j_2}^\top)=0$, for any $j_1\neq j_2\in[0:L]$;
\item[(b)] its $j$-th diagonal block $ \mathrm{Var}_{i}(H_i\xi_{i,j}) =\mathbb{E}_i\left[\mathbb{E}^+_{i,j}\left[ R_{i,j}^2 \right] \cdot H_{i,j}\,\mathrm{Var}^+_{i,j}(\Psi_{i,j})\,
H_{i,j}^\top\right]j}^\top}$, for any $j\in[0:L]$.
\end{enumerate}
\end{lemma}
This lemma inspires us to decompose the weight $H_{i,j}$ into two components: one to stabilize $\mathrm{Var}_{i,j}^+(\Psi_{i,j})$ that is known exactly, as it depends on the known behavior policy, and another to stabilize $\mathbb{E}^+_{i,j}\left[ R_{i,j}^2 \right]$ that can be estimated from data.
\subsection{Strong Gaussian Approximation}
\label{sec:strong_gaussian}
We now describe the construction of weights that ensure asymptotic normality of the estimator. This builds on the connection between the MDS $\{\xi_i\}$ and the estimation error $(\hat{\theta}_n - \theta^*)$, as established in Equation~\eqref{eq:link_error_mds}:
\begin{align}
\label{eq:link_mds_error_2}
B_n^\top A\left(\frac{1}{n}\sum_{i=1}^n H_i\xi_i\right)
=& B_n^\top A B_n(\hat\theta_n-\theta^*),\quad \mbox{where}\quad B_n:=-\frac{1}{n}\sum_{i=1}^n H_i J_i.
\end{align}
To establish the asymptotic normality of $\hat{\theta}_n$, it suffices to construct weights such that $\sum_{i=1}^n H_i \xi_i$ is asymptotically normal. Proposition~\ref{prop:martingale_clt_2} indicates that this requires stabilizing the conditional variance $\mathrm{Var}_i(H_i \xi_i)$. Lemma~\ref{lemma:mds} further shows that $\mathrm{Var}_i(H_i \xi_i)$ is block-diagonal, so it is enough to stabilize each stage-specific variance $\mathrm{Var}_i(H_{i,j} \xi_{i,j})$ such that $\sum_{i=1}^n \mathrm{Var}_i(H_{i,j} \xi_{i,j}) \to I$ for each stage $j$.
Lemma \ref{lemma:mds} shows that:
\begin{align}
\label{eq:residual_definition}
\mathrm{Var}_{i}(H_{i,j}\xi_{i,j})=\mathbb{E}_{i} \left[H_{i,j}\,(F_{i,j} \cdot\mathrm{Var}^+_{i,j}(\Psi_{i,j}))\,
H_{i,j}^\top\right],\quad\mbox{where}\quad F_{i,j}:=\mathbb{E}^+_{i,j}[R_{i,j}^2].
\end{align}
We thus can split $H_{i,j}$ into two parts:
\[
H_{i,j} = \hat{F}_{i,j}^{-\frac{1}{2}}\cdot W_{i,j},
\]
where the first part $\hat{F}_{i,j}^{-\frac{1}{2}}$ standardizes $F_{i,j}$ and is estimated using data available up to $\mathcal{F}_{i,j}^+$. Section \ref{sec:estimate_f} provides further discussion on estimating $\hat{F}_{i,j}$; the second part $W_{i,j}$ is designed to satisfy a stabilizing property below that controls the time-varying variance of $\mathrm{Var}^+_{i,j}(\Psi_{i,j})$.
\begin{property}[$(\alpha_1, \alpha_3)$-Stabilizing Weights]
\label{property:weight_stabilizing}
Set $W_{i,0}=1$ for all $i\in[1:n]$.
Given any $j\in[1:L]$, the weights $\{W_{i,j}\}_{i=1}^n$ are adapted to the filtration $\{\mathcal{F}_{i,j}^+\}_{i=1}^n$ and satisfy the following for a universal constant $c_2>0$ independent of the behavior policy:
\begin{enumerate}[(a)]
\item $\lambda_{\min}\left(\left\{\frac{1}{n}\sum_{i=1}^n \mathbb{E}_{i,j}\left[W_{i,j} \mathrm{Cov}^+_{i,j}(\Psi_{i,j},\Phi_{i,j})\right]\right\}^\top \left\{\frac{1}{n}\sum_{i=1}^n \mathbb{E}_{i,j}\left[W_{i,j} \mathrm{Cov}^+_{i,j}(\Psi_{i,j},\Phi_{i,j})\right]\right\}
\right) \geq ~ c_2^2\cdot n^{-\alpha_1}$;
\item $\frac{1}{n}\sum_{i=1}^n\mathrm{Tr}\left(\mathbb{E}\left[ W_{i,j}^\top W_{i,j}\right]\right) \leq c_2^{-1} \cdot n^{\alpha_1}$;
\item
$
\mathrm{Tr}\left( W_{i,j}\mathrm{Var}^+_{i,j}(\Psi_{i,j})W_{i,j}^\top\right)\leq c_2^{-1},\\
\frac{1}{n}
\sum_{i=1}^n \mathbb{E}
\left\|\mathbb{E}_{i,j}[W_{i,j}\mathrm{Var}^+_{i,j}(\Psi_{i,j})W_{i,j}^\top] - \mathbb{E}[W_{i,j}\mathrm{Var}_{i,j}^+(\Psi_{i,j})W_{i,j}^\top]
\right\|_{2} \leq c_2^{-1} \cdot n^{-\alpha_3},\\
\lambda_{\min}\left(\frac{1}{n} \mathbb{E}\left[\sum_{i=1}^n W_{i,j}\mathrm{Var}^+_{i,j}(\Psi_{i,j})W_{i,j}^\top\right]\right) \geq c_2.
$
\end{enumerate}
for a universal positive constant $c_2$, with $\alpha_1\in[0,\frac{1}{L+1}), \alpha_3>0$.
\end{property}
Property \ref{property:weight_stabilizing} looks
similar to Property \ref{property:weight_regularizing} but is stronger. In particular, to achieve asymptotic normality, Property \ref{property:weight_stabilizing}(c) further requires that after weight application, the resulted MDS variance is asymptotically standardized and stabilized around its expectation.
\begin{theorem}[Strong Gaussian Approximation]
\label{thm:be_feasible_2}
Suppose Assumptions~\ref{assump:exogeneity}, \ref{assump:linear_blip_function}, and \ref{assump:homoscedasticity} hold.
Suppose the sets $\mathcal{S},\mathcal{T},\mathcal{Y}$ are bounded.
Suppose the weights $H_{i,j}$ decompose as $H_{i,j} = \hat{F}_{i,j}^{-1/2} W_{i,j}$, where $W_{i,j}$ satisfies Property~\ref{property:weight_stabilizing}. The term $\hat{F}_{i,j}$ is an estimator of $F_{i,j}$, which is defined in Equation~\eqref{eq:residual_definition} and is bounded in $[\sigma^2, M^2]$ under Assumption~\ref{assump:homoscedasticity}. $\hat{F}_{i,j}$ is constructed using information available up to $\mathcal{F}^+_{i,j}$, and we assume that it satisfies the following condition:
\begin{equation}
\label{eq:f_convergence}
\frac{1}{n} \sum_{i=1}^n \mathbb{E}\left[ \left| \hat{F}_{i,j} - F_{i,j} \right| \right] = O(n^{-\gamma_F}) \mbox{ and }\hat{F}_{i,j} \in [\sigma^2, M^2].
\end{equation}
Let $\hat{\theta}_n$ be the AW-GMM estimator defined in \eqref{eq:gmm_estimator_linear}. Define $\widehat{\Xi}_n$ as a block-diagonal matrix with the $j$-th block equal to $n^{-1}\sum_{i=1}^n W_{i,j}\mathrm{Var}^+_{i,j}(\Psi_{i,j})W_{i,j}^\top$. Then there exists a Gaussian random variable $Z_\xi \sim \mathcal{N}(0, \widehat{\Xi}_n)$ such that for any convex set $C \in \mathcal{C}$ (the set of convex subsets of $\Theta$),
\begin{align}
\left| \mathbb{P}\left( \sqrt{n} (\hat{\theta}_n -\theta^*)\in C\right)-\mathbb{P}\left( (B_n^\top A B_n)^{\dagger}B_n^\top A Z_\xi\in C\right)\right|
= O \Big( L^2\cdot n^{-\min\big(\frac{1-\alpha_1}{12}, \frac{1-(L+1)\alpha_1}{2}, \frac{\alpha_3}{5}, \frac{\gamma_F}{5}\big)}\Big),
\end{align}
where for any matrix $M$, we use $M^\dagger$ to denote its Moore–Penrose pseudo-inverse.
\end{theorem}
\begin{remark}
Theorem~\ref{thm:be_feasible_2} not only establishes the asymptotic normality of $\hat{\theta}_n$, but also provides an upper bound on the uniform convergence rate of its strong Gaussian approximation. This result is stronger than the central limit theorem for $\hat{\theta}_n$, which can be achieved under a weaker condition: when $\hat{F}_{i,j}$ estimates $F_{i,j}$ consistently (the $L_1$ convergence is not required as in \eqref{eq:f_convergence}), one can invoke Proposition~\ref{prop:martingale_clt_2} to show that,\footnote{We omit such a proof for conciseness and since it is a weaker result than our uniform convergence.}
\begin{equation}
\label{eq:clt_feasible}
\sqrt{n} B_n^\top A B_n(\hat{\theta}_n-\theta^*) \Rightarrow \mathcal{N}(0,B_n^\top A \widehat{\Xi}_nA B_n ).
\end{equation}
However, the normalizing matrix $B_n^\top A B_n$ may not necessarily converge, particularly in the adaptive settings where the behavior policy is evolving over time; we thus cannot rely on \eqref{eq:clt_feasible} to construct confidence intervals for any arbitrary,
{data-independent,} projection of $\theta^*$ (for example, the treatment effect $\theta^*_j$ at stage $j$), unless we use the uniform convergence results provided in Theorem \ref{thm:be_feasible_2}.
\end{remark}
\subsubsection{Post-RL Inference.}
Note that both $B_n$ and $\widehat{\Xi}_n$ are measurable from the data. Theorem~\ref{thm:be_feasible_2} implies the following Gaussian approximation,
\begin{align}
\label{eq:strong_gaussian_approx}
\sqrt{n} (\hat{\theta}_n-\theta^*) \approx \mathcal{N}\left( 0,M_n
\right),\quad \mbox{where}\quad M_n:=(B_n^\top A B_n)^{\dagger}B_n^\top A\widehat{\Xi}_n
A B_n(B_n^\top A B_n)^{\dagger}.
\end{align}
This approximation enables the construction of confidence regions, formalized by the two corollaries below.
\begin{corollary}[Uniform Confidence Intervals]
\label{cor:uniform_ci}
Suppose the conditions in Theorem \ref{thm:be_feasible_2} hold.
For a confidence level $a\in (0, 1)$, for any projection $\ell\in\Theta$, define the confidence interval as
\[
CI_a=\bigg[\ell^\top \hat{\theta}_n\pm
n^{-\frac{1}{2}} q_{\ell,\frac{a}{2}}
\bigg],
\]
where $q_{\ell,\frac{a}{2}}$ is the $\frac{a}{2}$-quantile of $ \mathcal{N}\left(0,
\ell^\top M_n \ell
\right)$. It holds that:
\[
\sup_{\ell\in\mathbb{R}^{1+dL}}\left|\mathbb{P}(\ell^\top \theta\in CI_a) - (1-a)\right|=O\Big( L^2\cdot n^{-\min\left( \frac{1-\alpha_1}{12}, \frac{1-(L+1)\alpha_1}{2},\frac{\alpha_3}{5},\frac{\gamma_F}{5} \right)}\Big).
\]
\end{corollary}
\begin{corollary}[Simultaneous Confidence Band]
\label{cor:uniform_cb}
Suppose the conditions in Theorem \ref{thm:be_feasible_2} hold.
For a confidence level $a\in (0, 1)$, define the confidence band as
\[
CB_a = \left[\hat{\theta}_n \pm n^{-\frac{1}{2}} q_{\infty,1-a}\cdot \hat{d}_n^{-\frac{1}{2}}\right],
\]
where $\hat{d}_n$ is the matrix vector that concatenates the diagonal blocks of $M_n$, and $q_{\infty,1-a}$ is the $(1-a)$-quantile of $\|\mathcal{N}(0,(\hat{D}_n^{\dagger})^{\frac{1}{2}}
M_n
(\hat{D}_n^{\dagger})^{\frac{1}{2}})\|_\infty$ for $\hat{D}_n:=\mbox{diag}\{\hat{d}_n\}$.\footnote{
{By $\|N(0,A)\|_{\infty}$, we denote the random variable that corresponds to the maximum absolute value of any entry in a vector that is distributed according to $N(0,A)$.}} It holds that,
\[
\left|\mathbb{P}(\theta\in CB_a)-( 1-a)\right|= O\Big( L^2\cdot n^{-\min\left( \frac{1-\alpha_1}{12}, \frac{1-(L+1)\alpha_1}{2}, \frac{\alpha_3}{5}, \frac{\gamma_F}{5} \right)}\Big).
\]
\end{corollary}
\subsubsection{Estimating $F_{i,j}$.}
\label{sec:estimate_f}
We now provide a generic framework for estimating $F_{i,j}$ in a way that satisfies the condition in Equation~\eqref{eq:f_convergence}. Recall that $F_{i,j} := \mathbb{E}^+_{i,j}[R_{i,j}^2]$ is, by definition, a function of $(S_{i,1:j}, T_{i,1:j-1})$. If we further assume that this function is invariant to the behavior policy, then it can be consistently estimated using the realized data.
This condition holds, for instance, when the residual $R_{i,j}$ corresponds to the counterfactual outcome $Y(T_{i,1:j-1}, \pi^*_{j:L-1})$, which arises in our application to the Markovian model discussed in Section~\ref{sec:partial_linear_model}. In such cases, we denote this function by $f_j$:
\begin{assumption}\label{assump:invariant_fj}
The conditional expectation function $\mathbb{E}[R_{i,j}^2 \mid S_{i,1:j} = s_{1:j},\ T_{i,1:j-1} = \tau_{1:j-1}]$
is invariant to the episode index $i$ and is thus denoted by $f_j$ as a function of $(s_{1:j}, \tau_{1:j-1})$:
\begin{equation}
\label{eq:define_f}
f_j(s_{1:j}, \tau_{1:j-1}):=\mathbb{E}[R^2_{i,j}\mid S_{i, 1:j}=s_{1:j}, T_{i,1:j-1}=\tau_{1:j-1}].
\end{equation}
\end{assumption}
Let $\hat{f}_{i,j}$ be a predictor of $f_j$, fitted on the previous $i-1$ episodes and satisfying the estimation rate $i^{-\gamma_F}$ for $\gamma_F\in(0,1)$:
\begin{equation}
\label{eq:estimate_f}
\|\hat{f}_{i,j}- f_j\|_{1,\infty} = O(i^{-\gamma_F}),
\end{equation}
where we introduce the norm $\|\cdot\|_{1,\infty}$ for any $f$ defined on domain $\mathcal{S}^{j_1}\times\mathcal{T}^{j_2}\rightarrow \mathbb{R}$:
$\|f\|_{1,\infty} := \mathbb{E}_{s_1\sim P_S}\big[\sup_{(s_{2:j_1},\tau_{1:j_2})\in\mathcal{S}^{j_1-1}\times \mathcal{T}^{j_2}} |f(s_{1:j_1}, \tau_{1:j_2})|\big].$
Then
we have Condition \eqref{eq:f_convergence} satisfied:
\begin{align}
& \frac{1}{n}\sum_{i=1}^n\mathbb{E}[|\hat F_{i,j}-F_{i,j}|]=\frac{1}{n}\sum_{i=1}^n \mathbb{E}[|\hat{f}_{i,j}(S_{i,1:j}, T_{1,j-1})- f_j(S_{i,1:j}, T_{1,j-1})|] \leq \frac{1}{n}\sum_{i=1}^n \|\hat{f}_{i,j}- f_j\|_{1,\infty}
\nonumber\\
&= O\left( \frac{1}{n}\sum_{i=1}^n i^{-\gamma_F}\right)
\leq O\left( \Big(\frac{1}{n}\sum_{i=1}^{n} i^{-1} \Big)^{\gamma_F}\right) =
O((\log n/n)^{\gamma_F}),\label{eq:f_cumulative_convergence}
\end{align}
where the last inequality is due to that $x^{-\gamma_F}$ is a concave function with $\gamma_F\in(0,1)$.
Estimating $f_j$ to achieve the rate in \eqref{eq:estimate_f} can be addressed by online learning algorithms \citep{rakhlin2014online,daskalakis2022fast}, and in particular \citeauthor{daskalakis2022fast} propose a learning algorithm to achieve fast rates of convergence in nonparametric online regression. Later in Section~\ref{sec:partial_linear_model} we show how to estimate $f_{j}$ for high-dimensional Markovian models via Lasso (simpler than proposed algorithm in \cite{daskalakis2022fast}) and establish the guaranteed convergence rate in \eqref{eq:estimate_f}.
\subsection{Normality Weights for Decaying Exploration under Bilinear Features}
\label{sec:normality_application}
We now instantiate weight construction on asymptotic normality for common RL algorithms.
Similar to Section \ref{sec:consistency_application}, we adopt features as instruments and consider the bilinear feature map in Assumption~\ref{ass:bilinear}:
\begin{equation*}
\Phi_{i,j} = (T_{i,j} - \pi^*(S_{i,1:j}, T_{i,1:j-1})) \otimes X_{i,j}
\end{equation*}
where recall that $X_{i,j}$ serves as a shorthand for $\chi_j(S_{i,1:j}, T_{i,1:j-1})$.
Under the overlap condition specified in Assumption \ref{assump:overlap}, we provide weighting choices in the following corollary to achieve strong Gaussian approximation in Theorem \ref{thm:be_feasible_2}.
\begin{corollary}
\label{cor:normality}
Suppose Assumptions~\ref{assump:exogeneity},~\ref{assump:linear_blip_function},~\ref{ass:bilinear},~\ref{assump:overlap},~\&~\ref{assump:homoscedasticity} hold.
Define the random matrix:
\begin{equation}
V_{i,j}:=\mathbb{E}_{i,j}[X_{i,j}X_{i,j}^\top].
\end{equation}
Let weights $H_{i,j} := \hat{F}_{i,j}^{-1/2}W_{i,j}$, where $\hat{F}_{i,j}$ satisfies Condition \eqref{eq:f_convergence} and
\begin{equation}
\label{eq:w_normality}
W_{i,j}:=(I_{d_\tau}\otimes\hat{V}_{i,j}^{-1/2}) (I_{d_\tau}\otimes X_{i,j})\mathrm{Var}^+_{i,j}(T_{i,j})^{-1/2}
(I_{d_\tau}\otimes X_{i,j})^\top\cdot\|X_{i,j}\|_2^{-2},
\end{equation}
where $\|X_{i,j}\|_2$ denotes the $\ell_2$ norm of the vector $X_{i,j}$;
$ \hat{V}_{i,j}$ is positive definite adapted to $\mathcal{F}_{i,j}$, satisfies
$c^{-1}\cdot I \succeq\hat{V}_{i,j}\succeq c\cdot I$ and approximates $V_{i,j}$ (where $c$ and $V_{i,j}$ are introduced in Assumption \ref{assump:overlap}) with
\begin{equation}
\frac{1}{n}
\sum_{i=1}^n \mathbb{E}\left[\left\|\hat{V}_{i,j}- V_{i,j}
\right\|_{2}\right] = O(n^{-\gamma_V}).
\label{eq:v_convergence}
\end{equation}
Then, this $W_{i,j}$ satisfies Property \ref{property:weight_stabilizing} with $\alpha_1=\alpha$ and $\alpha_3=\gamma_V$, and thus we have the uniform Gaussian approximation rate in Theorem \ref{thm:be_feasible_2} to be $O(L^{2}n^{-\min(\frac{1-\alpha}{12},\frac{1-(L+1)\alpha}{2},\frac{\gamma_F}{5},\frac{\gamma_V}{5})})$.
\end{corollary}
The choice of $W_{i,j}$ in \eqref{eq:w_normality} requires estimating $V_{i,j}=\mathbb{E}_{i,j}[X_{i,j}X_{i,j}^\top]$.
When $j=1$, we have $V_{i,1}=\mathbb{E}\left[X_{i,1}X_{i,1}^\top\right]$, with $\{X_{i,1}\}$ being i.i.d.; thus we can estimate $V_{i,1}$ consistently via the sample-mean estimator and achieve $O(n^{-1/2})$ estimation rate. We hence focus on estimating $V_{i,j}$ for stage $j\geq 2$, where $V_{i,j}$ involves the expectation over $S_{i,j}$ conditional on $(S_{i,1:j-1}, T_{i,1:j-1})$. The randomness in $S_{i,j}$ arises from the state transition process, which is independent of the behavior policy. As a result, we can express $V_{i,j}$ as the output of a function $\nu_j(\cdot)$, with input $(S_{i,1:j-1}, T_{i,1:j-1})$. This function $\nu_j$ depends only on the stage index $j$ but not on the episode index $i$, i.e.,
\begin{equation}
\label{eq:define_v}
\nu_j(s_{1:j-1}, \tau_{1:j-1}):=\mathbb{E}[X_{i,j}X_{i,j}^\top\mid S_{i, 1:j-1}=s_{1:j-1}, T_{i,1:j-1}=\tau_{1:j-1}].
\end{equation}
Estimating $\nu_j$ can be similarly addressed by the online learning literature as we have discussed for estimating $f_j$ in Section \ref{sec:estimate_f}, which shall be instantiated with Lasso estimator for the high-dimensional Markovian models in the next section. Below we provide guarantees for generic estimators.
Let $\hat{\nu}_{i,j}$ be a generic predictor of $\nu_j$, fitted on the previous $i-1$ episodes and satisfying the estimation rate $i^{-\gamma_V}$, with $\gamma_V\in(0,1)$:
\begin{equation}
\label{eq:estimate_v}
\|\hat{\nu}_{i,j}- \nu_j\|_{1,\infty} = O(i^{-\gamma_V}).
\end{equation}
Above, we reload the norm notation $\|\cdot\|_{1,\infty}$ for any $v$ defined on domain $\mathcal{S}^{j_1}\times\mathcal{T}^{j_2}\rightarrow \mathbb{R}^{d_1\times d_2}$:
$\|v\|_{1,\infty} := \mathbb{E}_{s_1\sim P_S}\big[\sup_{(s_{2:j_1},\tau_{1:j_2})\in\mathcal{S}^{j_1-1}\times \mathcal{T}^{j_2}} \|v(s_{1:j_1}, \tau_{1:j_2})\|_{2}\big]$, where $\|\cdot\|_{2}$ denotes the matrix operator norm induced by vector $2$-norm. Then we have Condition \eqref{eq:v_convergence} satisfied:
\begin{align}
\label{eq:nu_cumulative_convergence}
&\frac{1}{n}
\sum_{i=1}^n \mathbb{E}\left[\left\|\hat{V}_{i,j}- V_{i,j}
\right\|_{2}\right]=\frac{1}{n}\sum_{i=1}^n \mathbb{E}[\|\hat{\nu}_{i,j}(S_{i,1:j-1}, T_{1,j-1})- \nu_j(S_{i,1:j-1}, T_{1,j-1})\|_2] \nonumber\\
& \leq \frac{1}{n}\sum_{i=1}^n \|\hat{\nu}_{i,j}- \nu_j\|_{1,\infty}=
O((\log n/n)^{\gamma_V}).
\end{align}
To this end, we provide a full theory with the generic nuisance estimators.
\begin{corollary}
\label{cor:full_normality}
Suppose Assumptions~\ref{assump:exogeneity},~\ref{assump:linear_blip_function},~\ref{ass:bilinear},~\ref{assump:overlap},~\ref{assump:homoscedasticity},~\&~\ref{assump:invariant_fj} hold.
Define functions $f_j$ and $\nu_j$ as in Equations~\eqref{eq:define_f} and \eqref{eq:define_v} respectively.
Suppose $\hat{f}_{i,j}$ and $\hat{\nu}_{i,j}$ are fitted using the first $i-1$ episodes and satisfy
$ \|\hat{f}_{i,j}- f_j\|_{1,\infty} = O(i^{-\gamma_F})$ and $\|\hat{\nu}_{i,j}- \nu_j\|_{1,\infty} = O(i^{-\gamma_V})$ respectively.
Set the weights $H_{i,j}$ as in Corollary \ref{cor:normality}, where we use $\hat{F}_{i,j}:=\hat{f}_{i,j}(S_{i,1:j}, T_{i,1:j-1})$ and $\hat{V}_{i,j}:=\hat{\nu}_{i,j}(S_{i,1:j-1}, T_{i,1:j})$. Then the uniform Gaussian approximation rate in Theorem \ref{thm:be_feasible_2} holds with rate $\widetilde O(L^{2}n^{-\min(\frac{1-\alpha}{12},\frac{1-(L+1)\alpha}{2},\frac{\gamma_F}{5},\frac{\gamma_V}{5})})$, where we use $\widetilde O(\cdot)$ to omit logarithm terms.
\end{corollary}
This result is a direct application of Corollary \ref{cor:normality}, Equations \eqref{eq:f_cumulative_convergence} and \eqref{eq:nu_cumulative_convergence}; we omit its proof for brevity.
\section{Application to High-dimensional Markovian Models}
\label{sec:partial_linear_model}
In this section, we instantiate the generic inference framework we present in Section \ref{sec:normality} to high-dimensional Markovian models. In particular, we consider the baseline policy of ``no treatment'' as the evaluation policy, i.e.~$0_{1:L}$, where the treatment $0_j$ represents a pre-existing status quo treatment that would have been applied at stage $j$ in the absence of the RL experimentation process.
Estimating the value of a no-treatment policy is highly valuable in practical RL scenarios. It enables verification, with statistical confidence, that the deployed RL policy resulted in a statistically significant higher reward compared to consistently applying the status quo treatment at each stage.
This allows to test an alternative policy, which introduces a new innovation, product feature, or treatment.
We consider high-dimensional Markovian models that satisfy the high-level assumptions on the data-generating process (DGP) introduced in previous sections. Our focus is on a specific instantiation of the estimation framework in which the features $\Phi_{i,j}$ are used as instruments $\Psi_{i,j}$, and the norm matrix $A$ in the GMM estimation is set to the identity. While our analysis centers on this setting, extensions to more general cases are straightforward. We introduce concrete statistical estimation algorithms for weight construction and provide the corresponding strong Gaussian approximation rates. Finally, we conduct simulations to empirically validate the inferential ability and robustness of our method in this setting.
\subsection{Data Generating Process}
We start by formalizing the DGP for high-dimensional Markovian models.
Let $\mathcal{S}, \mathcal{T}$ be the bounded state and treatment spaces.
The RL agent sequentially rolls out one episode per unit. Each episode $i$ consists of a length-$L$ trajectory of high-dimensional states, where only a sparse subset of state components influences the final outcome. Specifically, the blip functions depend only on low-dimensional sub-vectors of the state. To facilitate the weight construction discussed in Sections~\ref{sec:consistency_application} and~\ref{sec:normality_application}, we assume the blip function at each stage $j$ takes a bilinear form and is written as $T_{i,j}\otimes\chi(S_{i,j,\Omega})$ for Lipschitz and bounded $\chi(\cdot)$, where $\Omega$ denotes the set of relevant state coordinates. We use $d_\Omega$ and $d_s$ to denote the dimensions of $S_{i,j,\Omega}$ and $S_{i,j}$, respectively.
For each episode $i$, the behavior policy $\pi^{\text{obs}}_i$ is decided based on the previous $i-1$ episodes. This policy $\pi^{\text{obs}}_i$ is assumed to be known and assigns the treatment via $T_{i,j} := \pi^{\text{obs}}_i(S_{i,j}, \zeta_{i,j})$, where $\zeta_{i,j}$ is an i.i.d.~bounded noise term.
We consider the behavior policy satisfies Assumption \ref{assump:overlap}(a) such that:
\begin{equation}
\mathrm{Var}^+_{i,j}(T_{i,j}) \succsim i^{-\alpha}\cdot I_{d_\tau},
\label{eq:plmm_behavior_policy_rate}
\end{equation}
for an exploration rate $\alpha$.
State transitions are governed by unknown matrices $(A_j, B_j, M_j)$, which capture the effects of the current treatment, the current state, and the initial state, respectively, on the next state. Let $\eta_{i,j}$ denote the i.i.d.~mean-zero bounded noise associated with the state transitions.
The final outcome $Y_i$ is modeled as a sparse linear function of the last treatment, the last state, and the initial state, perturbed by an i.i.d.~bounded noise term $\epsilon_i$.
Model \ref{algo:plmdgp} summarizes the DGP.
\begin{model}
\DontPrintSemicolon
\KwIn{Time horizon $L$; behavior policy $\pi^{\text{obs}}_i$}
\KwOut{Episodic data $(S_{i,1}, T_{i,1}, \ldots, S_{i,L}, T_{i,L}, Y_i)$}
Observe $S_{i,1} \stackrel{\text{i.i.d.}}{\sim} P_S$
\For{$j=\{1,\dots, L-1\}$}{
Assign $T_{i,j} \gets \pi^{\text{obs}}_i(S_{i,j}, \zeta_{i,j})$, for $\zeta_{i,j} \stackrel{\text{i.i.d.}}{\sim} P_{\zeta,j}$
Observe $S_{i,j+1} \gets A_{j+1}^\top( T_j\otimes\chi(S_{i,j, \Omega})) + B_{j+1}^\top S_{i,j} +M_{j+1}^\top S_{i,1} + \eta_{i,j+1}$, for $\eta_{i,j+1} \stackrel{\text{i.i.d.}}{\sim} P_{\eta,j+1}$}
Assign $T_{i,L} \gets \pi^{\text{obs}}_i({S}_{i,L}, \zeta_{i,L})$, for $\zeta_{i,L} \stackrel{\text{i.i.d.}}{\sim} P_{\zeta, L}$
Observe $\ensuremath{Y}_i \gets \alpha^\top (T_{i,L}\otimes\chi(S_{i,L,\Omega}))) + \beta^\top S_{i,L} + \kappa_L^\top S_{i,1} + \epsilon_i$, for $\epsilon_i \stackrel{\text{i.i.d.}}{\sim} P_{\epsilon}$
\caption{{\sc PLM DGP}: high-dimensional Markovian Data Generating Process}
\label{algo:plmdgp}
\end{model}
\begin{lemma}
\label{lemma:partial_linear_model}
Consider the no-treatment $0_{1:L}$ policy as the evaluation policy. The high-dimensional Markovian model specified in Model \ref{algo:plmdgp} satisfies Assumptions \ref{assump:exogeneity}, \ref{assump:linear_blip_function}, \ref{ass:bilinear}, \& \ref{assump:homoscedasticity} .
\end{lemma}
Given the DGP outlined in Model \ref{algo:plmdgp} , we can expand the outcome $\ensuremath{Y}$ recursively:
\begin{align}
\label{eq:final_unroll_plm}
&\quad \quad \ensuremath{Y}_i = \sum_{j'=j}^L\theta_{j'}^\top ( T_{i,j'}\otimes\chi(S_{i,j',\Omega})) + \beta_j^\top S_{i,j} + \kappa_j^\top S_{i,1} +
\epsilon_{i,j}\\
\mbox{where}\quad & \theta_L^*:=\alpha,\quad \beta_L:=\beta; \quad \mbox{For} ~j=L,\cdots,2: \nonumber\\
& \quad \beta_{j-1} := B_{j}\beta_j,\quad \theta_{j-1}^*:= A_{j}\beta_j, \quad\kappa_{j-1} := \kappa_j + M_{j}\beta_j, \quad \epsilon_{i,j}= \sum_{j'=j}^{L-1}\beta_{j'+1}^\top \eta_{i,j'} + \epsilon_i.\nonumber
\end{align}
Our goal is to estimate the structural parameters $\{\theta_0^*, \theta_1^*,\dots, \theta_L^*\}$, where
$\theta_0^*:=\mathbb{E}[Y(0_{1:L})]=\beta_1^\top\mathbb{E}[S_{i,1}] + \kappa_1^\top\mathbb{E}[ S_{i,1}]$ represents the expected outcome under the baseline policy,
and $\theta_{1:L}^*$ reflect the stage-wise dynamic treatment effects in Model \ref{algo:plmdgp} outlined in \eqref{eq:final_unroll_plm}.
\subsection{Constructing the Weights}
\label{sec:estimation_plm}
We now provide details on the weight construction for Model \ref{algo:plmdgp}.
With the bilinear blip functions,
Corollary~\ref{cor:consistency} provides analytical weight construction to achieve consistency.
To achieve the strong Gaussian approximation, Corollary~\ref{cor:full_normality} shows that we need to estimate functions $f_j$ and
$\nu_j$
to satisfy Conditions \eqref{eq:estimate_f} \& \eqref{eq:estimate_v} respectively. Recall that in general settings, $f_{j}$ and $\nu_{j}$ can be estimated by online learning algorithms, while here under the high-dimensional Markovian models, this estimation procedure can be greatly simplified.
\subsubsection{Estimating $f_{j}$.}
We first show that Model \ref{algo:plmdgp} satisfies the below homoscedasticity property, which satisfies Assumptions \ref{assump:homoscedasticity} \& \ref{assump:invariant_fj} and makes estimating $f_{j}$ a regular estimation problem.
\begin{lemma}[Homoscedasticity]
\label{lemma:homoscedasticity_plm}
Consider the no-treatment $0_{1:L}$ policy as the evaluation policy. For Model~\ref{algo:plmdgp},
the residual $R_{i,j}$, when conditioning on the $(S_{i,1:j}, T_{1:j-1})$, is independent from the behavior policy and regardless of the realization of $(S_{i,1:j}, T_{1:j-1})$, has the same variance, denoted as $\sigma_j^2$.
\end{lemma}
With this property, we can decompose $F_{i,j}$ into two parts:
\begin{equation}
\label{eq:decomposition_of_f}
F_{i,j}:=\mathbb{E}^+_{i,j}\left[R_{i,j}^2\right]= G_{i,j}^2 + \sigma_j^2, \quad\mbox{with}\quad G_{i,j}:=\mathbb{E}^+_{i,j}\left[R_{i,j}\right]\quad \mbox{and}\quad \sigma_j^2=\sigma_{i,j}^2:=\mathrm{Var}^+_{i,j}(R_{i,j}).
\end{equation}
Note that $G_{i,j}$, when conditioning on $(S_{i, 1:j}, T_{i,1:j-1})$, does not depend on the behavior policy. Thus we can view $G_{i,j}$ as a function $g_j$ of $(S_{i,1:j}, T_{i,1:j-1})$, which only depends on the stage index $j$ but not the specific episode index $i$, that is,
\begin{equation}
\label{eq:g_hdmm}
g_j(S_{i,1:j}, T_{i,1:j-1}):=G_{i,j} = \mathbb{E}[R_{i,j}\mid S_{i, 1:j}, T_{i,1:j-1}]= S_{i,j}^\top\beta_j + S_{i,1}^\top \kappa_j,
\end{equation}
where the last equation is due to that $R_{i,j}=\beta_j^\top S_{i,j} + \kappa_j^\top S_{i,1} + \epsilon_{i,j}$ by Equation~\eqref{eq:final_unroll_plm}.
We can estimate $g_j$ by regressing $\{R_{i,j}\}$ on $\{(S_{i,1}, S_{i,j})\}$ to get $\hat{g}_{i,j}$; we then estimate
$\sigma_{j}^2$ by its sample variance $\hat\sigma_{i,j}^2$ up till episode~$i$;
finally we set $\hat{f}_{i,j}(\cdot):=\mathrm{Clip}_{[\sigma^2, M^2]}\left(\hat{g}_{i,j}(\cdot)^2 + \hat\sigma_{i,j}^2\right)$, where $(\sigma^2,M^2)$ denotes the specified range for $F_{i,j}$.
Here, for any real numbers $a < b$ and any $x \in \mathbb{R}$, the clipping operator is defined as:
\begin{equation}
\label{eq:clip_scalar}
\text{Clip}_{[a,b]}\left(x\right):=\max(\min(x,b), a).
\end{equation}
Algorithm \ref{algo:f} summarizes the estimation process. Note that
we cannot observe the true residual $R_{i,j} =Y_{i}-\sum_{j'>j}\Phi_{i,j'}^\top\theta^*_j $, which requires the knowledge of the true structure parameter $\theta^*$. However, we can estimate it using an approximated $\hat{\theta}$ (for example, the consistent estimate from Section \ref{sec:consistency}). With that, we obtain the estimation rate of $\hat{f}_{i,j}$ as $\|\hat f_{i,j}-f_j\|_{1,\infty}=O(i^{-\frac{1-L\alpha}{2}})$; see details in Appendix~\ref{appendix:proof_hdmm_gaussian}.
\setcounter{algocf}{0}
\begin{algorithm}
\KwIn{The first $i-1$ episodes $\{(S_{i',1}, T_{i',1}, \ldots, S_{i',L}, T_{i',L}, Y_{i'})\}_{i'=1}^{i-1}$; stage index $j$; estimated $\hat{\theta}_{1:L}$ ; specified range $[\sigma^2, M^2]$ for $F_{i,j}$.}
\KwOut{Estimated $\hat{f}_{i,j}$.}
Set approximated residual as $\hat{R}_{i',j}=Y_{i'} - \sum_{j'=j}^L(T_{i',j'}\otimes \chi(S_{i',j',\Omega}))^\top\hat{\theta}_{j'}$ for $i'\in [1:i-1]$.
Estimate $g_{j}$, defined in \eqref{eq:g_hdmm}, by regressing $\{\hat{R}_{i',j}\}_{i'=1}^{i-1}$ on $\{(S_{i',1}, S_{i',j})\}_{i'=1}^{i-1}$ using Lasso regression:
\begin{equation}
(\hat{\beta}_{i,j}, \hat{\kappa}_{i,j})= \argmin_{\tilde{\beta}, \tilde{\kappa}} \frac{1}{i-1}\sum_{i'=1}^{i-1} \left(\hat{R}_{i',j}-\tilde{\beta}^\top S_{i',j} - \tilde{\kappa}^\top S_{i',1}\right)^2 + \lambda_g (\|\tilde{\beta}\|_1+ \|\tilde{\kappa}\|_1)\label{eq:lasso_g}.
\end{equation}
Set the estimate as $\hat{g}_{i,j}(s_1, s_j)=s_j^\top \hat{\beta}_{i,j} + s_1^\top \hat{\kappa}_{i,j}$.
Calculate $
\hat{\sigma}_{i,j}^2 := \frac{1}{i-1}\sum_{i'=1}^{i-1}\left( \hat{R}_{i',j}- \hat{g}_{i,j}(S_{i',1}, S_{i',j} )\right)^2.
$
Set $\hat{f}_{i,j}(\cdot)=\text{Clip}_{[\sigma^2, M^2]}\left(
\hat{g}_{i,j}(\cdot)^2+\hat{\sigma}_{i,j}^2\right)$ with clipping operator defined in \eqref{eq:clip_scalar}.
\caption{Estimating $f_j$ under Model \ref{algo:plmdgp}}
\label{algo:f}
\end{algorithm}
\subsubsection{Estimating $\nu_{j}$.} We now estimate
\[\nu_j(s_{1:j-1},\tau_{1:j-1}):=\mathbb{E}\left[\chi(S_{i,j,\Omega})\chi(S_{i,j,\Omega})^\top\mid S_{i,1:j-1}=s_{1:j-1}, T_{i,1:j-1}=\tau_{1:j-1}\right].\]
Note that when $j=1$, we have $\nu_{1}\equiv\mathbb{E}\left[\chi(S_{i,1,\Omega})\chi(S_{i,1,\Omega})^\top\right]$, with $\{S_{i,1}\}$ being i.i.d.; thus we can estimate $\nu_{1}$ consistently via the sample-mean estimator and achieve $O(n^{-1/2})$ estimation rate. We hence focus on estimation of $\nu_{j}$ for stage $j\geq 2$.
Let $h_{j}(S_{i,1:j-1}, T_{i,1:j-1})$ denote the conditional expectation of $S_{i,j}$ when restricted to coordinates $\Omega$~(since only those coordinates enter into $\psi_j$ and contribute to $V_j$):
\begin{align}
& h_{j}(S_{i,1:j-1}, T_{i,1:j-1}):=\mathbb{E}_{i,j}[S_{i,j,\Omega}]=A_{j,\Omega}^\top\left(T_{i,j-1}\otimes\chi(S_{i,j-1, \Omega})\right) +B_{j-1,\Omega}^\top S_{i,j-1} + M_{j,\Omega}^\top S_{i,1},\label{eq:h_hdmm}
\end{align}
Then $\nu_j$ can be induced by $h_j$:
$
\nu_j(\cdot)
=~ \mathbb{E}_{i,j}\left[\chi\left(h_{j}(\cdot) + \eta_{i,j,\Omega}\right)\chi\left(h_{j}(\cdot) + \eta_{i,j,\Omega}\right)^\top\right],
$
for i.i.d.~exogenous noise $\eta_{i,j,\Omega}$ during the state transition.
We therefore first regress $\{S_{i,j,\Omega}\}$ on $\{S_{i,1}, S_{i,j-1},T_{i,j-1}\}$ to get estimate $\hat{h}_{i,j}$, based on which we obtain estimate $\hat \nu_{i,j}$, as summarized in Algorithm \ref{algo:v}.
Note that in Algorithm~\ref{algo:v}, we also apply the clipping operator when obtaining $\hat{\nu}_{i,j}$ to ensure its eigenvalues lie within the specified range for the true matrix $V_{i,j}$.
Here, For any real numbers $a < b$ and a symmetric matrix $X = U \Lambda U^\top \in \mathbb{R}^{d \times d}$, where $U$ is an orthogonal matrix and $\Lambda = \text{diag}(\lambda_1, \dots, \lambda_d)$ is the diagonal matrix of eigenvalues, the clipping operator is defined as:
\begin{equation}
\label{eq:clip_matrix}
\text{Clip}_{[a,b]}(X):= U\mbox{diag}\{
\text{Clip}_{[a,b]}(\lambda_1),\dots,
\text{Clip}_{[a,b]}(\lambda_d)U^\top
\},
\end{equation}
where recall $\text{Clip}_{[a,b]}(\lambda) := \max(\min(\lambda, b), a)$ for any scalar $\lambda \in \mathbb{R}$ as defined \eqref{eq:clip_scalar}.
The estimation rate of $\hat{\nu}_{i,j}$ satisfies $\|\hat\nu_{i,j}-\nu_j\|_{1,\infty}=O(i^{\frac{-\gamma(1-2\alpha)}{2}})$ (proof deferred to Appendix \ref{appendix:proof_hdmm_gaussian}).
\begin{algorithm}
\KwIn{Data $\{S_{i',1}, T_{i',1}, \ldots, S_{i',L}, T_{i',L}, Y_{i'}\}_{i'=1}^i$; stage index $j$; specified range $[c,c^{-1}]$ for eigenvalues of $V_{i,j}$.}
\KwOut{Estimated $\hat{\nu}_{i,j}$.}
Estimate $h_{j}$, defined in \eqref{eq:h_hdmm}, by regressing $\{S_{i',j, \Omega}\}_{i'=1}^i$ on $\{(S_{i',1}, S_{i',j-1}, T_{i',j-1})\}_{i'=1}^{i-1}$ using Lasso regression:
\begin{align}
(\hat{A}_{i,j,k}, \hat{B}_{i,j,k}, \hat{M}_{i,j,k}) &= \argmin_{(\tilde{A}_k, \tilde{B}_k, \tilde{M}_k)} \frac{1}{i-1}\sum_{i'=1}^{i-1} \left(S_{i',j,k}- \tilde{A}_k^\top \Phi_{i',j-1}- \tilde{B}_k^\top S_{i',j-1} - \tilde{M}_k^\top S_{i',1})\right)^2\nonumber \\
&\quad\quad\quad\quad\quad\quad\quad\quad + \lambda_k ( \|\tilde{A}_k\|_1 +\|\tilde{B}_k\|_1 +\|\tilde{M}_k\|_1 ), \quad\quad \mbox{for}\quad k\in\Omega.\label{eq:lasso_h}
\end{align}
Set the estimate $\hat{h}_{i,j}(s_1, s_{j-1},\tau_{j-1})=\hat{A}_{i,j,\Omega}^\top\left(\tau_{j-1}\otimes\chi(s_{j-1, \Omega})\right) +\hat B_{i,j,\Omega}^\top s_{j-1} +\hat{M}_{i,j,\Omega}^\top s_1$.
Define $\hat\nu_{i,j}(\cdot)=\mathrm{Clip}_{[c,c^{-1}]}\big(\frac{1}{i-1}\sum_{i'=1}^{i-1} \chi\big(\hat{h}_{i,j}(\cdot) + S_{i',j,\Omega}- \hat{h}_{i,j}(S_{i',1:j-1}, T_{i,1:j-1})\big)\chi\big(\hat{h}_{i,j}(\cdot) + S_{i',j,\Omega}- \hat{h}_{i,j}(S_{i',1:j-1}, T_{i,1:j-1})\big)^\top\big)$, with clipping operator defined in \eqref{eq:clip_matrix}.
\caption{Estimating $\nu_{j}$ under Model \ref{algo:plmdgp}}
\label{algo:v}
\end{algorithm}
\subsubsection{Putting everything together.}
\label{sec:combine_all_plmm}
With the estimated $\hat{f}_{i,j}$ and $\hat{\nu}_{i,j}$ in the previous sections, we are ready to construct weights to achieve strong Gaussian approximation in Corollary \ref{cor:full_normality}. We summarize the full steps in
Algorithm \ref{algo:plmdgp_weights} and provide the strong Gaussian approximation rate in Corollary \ref{cor:plmm_normality_rate}.
\begin{algorithm}
\DontPrintSemicolon
\KwIn{Data $\mathcal{D}=\{S_{i,1}, T_{i,1}, \ldots, S_{i,L}, T_{i,L}, Y_i\}_{i=1}^n$; true $\theta^*$ range $[-U,U]$; true $F_{i,j}$ range $[\sigma^2, M^2]$.}
\KwOut{Consistent Estimation $\hat{\theta}^{(C)}_{n}$ and Asymptotically Normal Estimation $\hat{\theta}_n^{(N)}$}
\textsc{//Part I: Construct consistent estimation}
Set $H_{i,0}^{(C)}:=1,\forall i=1,\dots, n$.
\For{$i=\{1,\dots, n\}, j=\{1,\dots, L\}$}{
Set $H_{i,j}^{(C)}:= (I_{d_\tau}\otimes \chi(S_{i,j,\Omega}))\mathrm{Var}^+_{i,j}(T_{i,j})^{-1/2}
(I_{d_\tau}\otimes \chi(S_{i,j,\Omega}))^\top$.
}
\For{$i=\{1,\dots, n\}, j=\{1,\dots, L\}$}{
Let $\tilde{\theta}^{(C)}_{i,j}$ solve the AW-GMM estimator \eqref{eq:gmm_estimator_linear} using
data $\{Z_{i'}\}_{i'=1}^{i-1}$ and weights $\{H_{i',j}^{(C)}\}_{i'=1}^{i-1}$.
}
Set the consistent estimation $\hat{\theta}^{(C)}_{n}$ as $\tilde{\theta}^{(C)}_{n,L}$.
\hrulefill
\textsc{//Part II: Construct asymptotically normal estimation}
\For{$i=\{1,\dots, n\}$}{
Set $\hat{\sigma}_{i,0}^2 := \frac{1}{n_{i,j}^F}\sum_{i'\in\mathcal{I}_{i,j}^F}\left(Y_{i'}-\sum_{j=1}^L
(T_{i',j}\otimes\chi(S_{i',j,\Omega}))^\top
\tilde{\theta}_{i,j}^{(C)}\right)^2$.
Set $H_{i,0}^{(N)}=\hat{\sigma}_{i,0}^{-1}$.
}
\For{$i=\{1,\dots, n\}$, $j=\{1,\dots, L\}$}{
Obtain $\hat{f}_{i,j}$ via Algorithm \ref{algo:f}$(\mathcal{D}, i,j,\tilde{\theta}_{i,j}^{(C)}, [\sigma^2, M^2])$.
Obtain $\hat{\nu}_{i,j}$ via Algorithm \ref{algo:v}$(\mathcal{D}, i,j)$.
Set $H_{i,j}^{(N)}=\hat{f}_{i,j}(S_{i,1:j}, T_{i,1:j-1})^{-1/2}(I_{d_\tau}\otimes\hat{\nu}_{i,j}(S_{i,1:j-1}, T_{i,1:j-1})^{-1/2}) (I_{d_\tau}\otimes \chi(S_{i,j,\Omega}))\mathrm{Var}^+_{i,j}(T_{i,j})^{-1/2}
(I_{d_\tau}\otimes \chi(S_{i,j,\Omega}))^\top\cdot\|\chi(S_{i,j,\Omega})\|_2^{-2}$.
}
Set the asymptotically normal estimation $\hat{\theta}^{(N)}_{n}$ as the solution to the AW-GMM estimator \eqref{eq:gmm_estimator_linear} using the full data $\mathcal{D}$ and weights $\{H_{i,j}^{(N)}\}_{i=1}^n$.
\caption{Post RL Estimation and Inference under Model \ref{algo:plmdgp}}
\label{algo:plmdgp_weights}
\end{algorithm}
\begin{corollary}
\label{cor:plmm_normality_rate}
Consider Model \ref{algo:plmdgp} and the no-treatment $0_{1:L}$ policy as the evaluation policy. Suppose that the nuisance components $g_j$ and $h_{j}$ are estimated via Lasso regression as in Lemma \ref{lemma:estimation_rate_plmm}. Under Assumption \ref{assump:overlap} with behavior exploration rate $\alpha\in[0,\frac{1}{L+1})$,
the AW-GMM estimation $\hat{\theta}_n^{(N)}$ given by Algorithm \ref{algo:plmdgp_weights} is asymptotically normal with strong Gaussian approximation rate of $\widetilde O\big(L^2 n^{-\min\big(\frac{1-\alpha}{12}, \frac{1-L\alpha}{10}, \frac{1-(L+1)\alpha}{2}\big)}\big)$.
\end{corollary}
\begin{figure}[t]
\centering
\includegraphics[width=\textwidth]{plots/gmm/cov_gmm.pdf}
\includegraphics[width=\textwidth]{plots/gmm/ci_gmm.pdf}
\caption{
Inference results of AW-GMM Estimations with different weights across varying sample size. Error bars are $95\%$ confidence intervals derived from $10^3$ simulations.
Results under \textsc{Oracle}~weights are shown in dashed line to indicate that oracle weighting requires knowledge of ground truth structure parameters and thus cannot be applied in practice.
AW-GMM Estimations with \textsc{Oracle} and \textsc{Feasible} weights meet nominal coverage, while the \textsc{Naive}~and \textsc{Consistent}~are either under- or over-coverage.
}
\label{fig:cov}
\end{figure}
\begin{figure}[t]
\centering
\includegraphics[width=\textwidth]{plots/gmm/tstats_hist_gmm.pdf}
\caption{Histogram of studentized statistics from Gaussian approximation \eqref{eq:strong_gaussian_approx} at sample szie $n=5\times 10^3$. Numbers are aggregated from $10^3$ simulations. AW-GMM Estimations with \textsc{Oracle} and \textsc{Feasible} weights are asymptotically normal.}
\label{fig:tstat}
\end{figure}
\begin{figure}[t]
\centering
\includegraphics[width=\textwidth]{plots/gmm/mse_gmm.pdf}
\includegraphics[width=\textwidth]{plots/gmm/bias_gmm.pdf}
\caption{Estimation results of AW-GMM Estimations with different weights across varying sample size. Error bars are $95\%$ confidence intervals derived from $10^3$ simulations. Results under \textsc{Oracle}~weights are shown in dashed line to indicate that oracle weighting requires knowledge of ground truth structure parameters and thus cannot be applied in practice. AW-GMM Estimations with \textsc{Oracle}~and \textsc{Feasible}~weights provide more accurate estimations for policy value $\theta_0^*$.}
\label{fig:estimation}
\end{figure}
\begin{figure}[t]
\centering
\includegraphics[width=\textwidth]{plots/gmm/misspec_last_all_gmm.pdf}
\caption{Estimation and inference results of AW-GMM Estimations with different weights under mis-specification at sample size $n=5\times 10^3$. We use polynomial approximations from degrees 1 to 5 for exponential feature mappings.
Error bars are $95\%$ confidence intervals derived from $10^3$ simulations.
AW-GMM Estimations under \textsc{Feasible} weights achieve nominal coverage and tight confidence intervals for degrees above one, consistently offering more accurate estimations with lower MSE and bias across all approximation degrees.}
\label{fig:misspecifcation}
\end{figure}
\subsection{Numerical Experiments}
We finally present empirical evidence supporting our method's effectiveness in scenarios with both correct and mis-specified feature mapping function $\phi_j$. Our findings illustrate that applying weights enhances outcomes in both cases: under correct specification, it empirically confirms our method's consistency and asymptotic normality across all structural parameters. In cases of mis-specification, it notably improves the accuracy of estimating the evaluation policy value.
We study a two-stage high-dimensional Markovian model (with $L=2$), as outlined in Model \ref{algo:plmdgp}, focusing on binary treatment scenarios. Data collection is performed by an $\epsilon$-greedy episodic RL agent, which sequentially rolls out an episode for each unit. The agent's behavior policy undergoes batch updates, with each batch including $100$ units; the exploration amount $\epsilon$ decays at a polynomial rate such that for any given batch index $b$, the exploration $\epsilon$ is set to be $b^{-0.5}/2$, satisfying Assumption \ref{assump:overlap}. Throughout the experiment, a total of $5000$ units are collected.
We consider high dimensional states and low dimensional features, with sparse linear models for both the state transition and the final outcome.
Specifically, the state space is in $\mathbb{R}^{10}$ with $S_{i,j}= (S_{i,j,1},\dots, S_{i,j,10})$; only the first coordinate is informative, with others being noises.
We focus on GMM-estimations using four different weighting schemes:
\begin{itemize}
\item \textsc{Naive}: Standard GMM estimation with no weights applied.
\item \textsc{Consistent}: AW-GMM estimation $\hat{\theta}^{(C)}_{n}$ under consistency weights as in Algorithm \ref{algo:plmdgp_weights}.
\item \textsc{Oracle}: AW-GMM estimation with oracle weights, using ground truth $F_{i,j}$ and $V_{i,j}$ in Corollary \ref{cor:normality}.
\item \textsc{Feasible}: AW-GMM estimation $\hat{\theta}_n^{(N)}$ under asymptotic normality weights as in Algorithm~\ref{algo:plmdgp_weights}.
\end{itemize}
\medskip
Note that the \textsc{Oracle}~weights are infeasible, as they require knowledge of the true data-generating process. Conversely, the \textsc{Naive}, \textsc{Consistent}, and \textsc{Feasible}~weights can either be directly computed or estimated from the data. Below we show that \textsc{Feasible}~weighting scheme performs comparably to the \textsc{Oracle}~and significantly outperforms the other two in achieving asymptotic normality for post-RL inference; the \textsc{Feasible}~also demonstrates robustness under misspecification.
\subsubsection{Estimation and Inference Validity.}
We first consider cases with correctly specified feature mapping.
Define the true feature mapping as $\phi_j(S_{i,1:j}, T_{i,1:j})=T_{i,j}\cdot (S_{i,j,1}, S_{i,j,1}^2)\in\mathbb{R}^2$.
Our estimand, $\theta^*=(\theta_0^*, \theta_{1,1}^*, \theta_{1,2}^*, \theta_{2,1}^*, \theta_{2,2}^*)$, includes the evaluation policy value $\theta_0^*$ and structure parameters $(\theta_{1,1}^*, \theta_{1,2}^*)$ and $(\theta_{2,1}^*, \theta_{2,2}^*)$, indicating the effect of the treatment at each stage.
Figure \ref{fig:cov} shows that
estimations with \textsc{Oracle}~and \textsc{Feasible}~weights achieve nominal coverage, contrasting with the \textsc{Naive}~no-weighting or \textsc{Consistent}~weights, which either have low coverage or are overly conservative. In particular, the \textsc{Naive}~estimator's confidence intervals for the policy value $\theta_0^*$ and the first-stage treatment effect $(\theta_{1,1}^*,\theta_{1,2}^*)$ even widen with increased sample size.
Figure \ref{fig:tstat} further shows that the studentized statistics of the AW-GMM estimators, derived from Eq.\eqref{eq:strong_gaussian_approx}, under \textsc{Oracle}~and \textsc{Feasible}~weights conform to asymptotic normality, unlike those from other weighting schemes.
Furthermore, Figure \ref{fig:estimation} shows that while all estimators yield similar results for the stage-wise treatment effect estimation, those with \textsc{Oracle}~and \textsc{Feasible}~weights achieve higher accuracy in estimating the policy value $\theta_0^*$ with smaller MSE and bias.
\subsubsection{Robustness under Misspecification.} Transitioning to cases of mis-specification, we consider a true feature mapping defined as $\phi_j(S_{i,1:j}, T_{i,1:j})=\exp(T_{i,j}S_{i,j,1}/2) -1$, unknown to the AW-GMM estimators.
These estimators then use polynomial approximations with degree $d\in[1:5]$ instead of the true $\phi_j$, i.e., we use $\hat{\phi}^{(d)}_j(S_{i,1:j}, T_{i,1:j})=T_{i,j}(S_{i,j,1}, \dots, S_{i,j,1}^d)\in\mathbb{R}^d$ as the approximated feature mapping for each degree $d$.
Under this scenario, the stage-wise structural parameters lose their causal interpretation, yet estimating the evaluation policy value $\theta_0^*$ remains relevant. Figure \ref{fig:misspecifcation} shows that AW-GMM estimations under \textsc{Feasible}~weights outperform those with \textsc{Naive}~no-weighting or \textsc{Consistent}~weights, demonstrating better coverage, narrower confidence intervals, and reduced MSE and bias.
Interestingly, with \textsc{Feasible}~weights, improvements in inference validity (coverage) and efficiency (confidence interval length) plateau for approximation degrees $d \geq 3$, while estimation quality (MSE and bias) slightly deteriorates at $d=5$.
The decrease in estimation accuracy at higher approximation degrees is due to the requirement to estimate more structural parameters from the same sample size, which instead complicates the estimation process.
This observation highlights a balance between inference and estimation: without knowing the precise feature mapping, selecting an approximation with appropriate complexity is crucial in optimizing both estimation accuracy and inferential robustness.
\bibliographystyle{informs2014}
\bibliography{reference}
\newpage
\begin{APPENDICES}
\section{Proofs of Main Lemmas}
\subsection{Proof of Lemma \ref{lemma:identification_policy_value}}
\label{appendix:proof_identification_policy_value}
We follow the proof pattern for Lemma 6 in \cite{lewis2020double}. Note that for the observed $\ensuremath{Y}$, we always have $\ensuremath{Y}\equiv \ensuremath{Y}(T_{1:L})$. For any stage $j$, we have
\begin{align*}
\theta_j^\top \phi_j(S_{1:j}, T_{1:j})\stackrel{(i)}{=} \gamma(S_{1:j}, T_{1:j}) \stackrel{(ii)}{=}\mathbb{E}\left[ \ensuremath{Y}(T_{1:j}, \pi^*_{j+1:L}) - \ensuremath{Y}(T_{1:j-1}, 0, \pi^*_{j+1:L}) \mid S_{1:j}, T_{1:j}\right] ,
\end{align*}
where (i) is by Assumption \ref{assump:linear_blip_function}, (ii) is by the definition of blip function. Moreover, note that:
\begin{align*}
\theta_j^\top \phi_j(S_{1:j}, T_{1:j-1}, \pi^*(S_{1:j}, T_{1:j-1})) =~& \gamma(S_{1:j}, T_{1:j-1}, \pi^*(S_{1:j}, T_{1:j-1}))\\
=~& \mathbb{E}\left[ \ensuremath{Y}(T_{1:j}, \pi^*_{j+1:L}) - \ensuremath{Y}(T_{1:j-1}, 0, \pi^*_{j+1:L}) \mid S_{1:j}, T_{1:j-1}, T_j=\pi^*(S_{1:j}, T_{1:j-1})\right]\\
=~& \mathbb{E}\left[ \ensuremath{Y}(T_{1:j-1}, \pi^*_{j:L}) - \ensuremath{Y}(T_{1:j-1}, 0, \pi^*_{j+1:L}) \mid S_{1:j}, T_{1:j-1}, T_j=\pi^*(S_{1:j}, T_{1:j-1})\right]\\
\stackrel{(iii)}{=}~& \mathbb{E}\left[ \ensuremath{Y}(T_{1:j-1}, \pi^*_{j:L}) - \ensuremath{Y}(T_{1:j-1}, 0, \pi^*_{j+1:L}) \mid S_{1:j}, T_{1:j-1}, T_j\right],
\end{align*}
where (iii) follows since the counterfactual outcomes are independent of the value of $T_j$ conditional on $S_{1:j}, T_{1:j-1}$ by Sequential Conditional Exogeneity Assumption~\ref{assump:exogeneity}. Subtracting the two equalities, and by the definition of $\Phi_j$, we derive that:
\begin{align}
\mathbb{E}\left[ \ensuremath{Y}(T_{1:j}, \pi^*_{j+1:L}) - \ensuremath{Y}(T_{1:j-1}, \pi^*_{j:L}) \mid S_{1:j}, T_{1:j}\right] = \theta_j^\top \Phi_j
\end{align}
Then we have
\begin{align*}
\mathbb{E}\left[ \ensuremath{Y}-\ensuremath{Y}(T_{1:j-1}, \pi^*_{j:L})|S_{1:j}, T_{1:j} \right] =~& \mathbb{E}\left[ \ensuremath{Y}(T_{1:L})-\ensuremath{Y}(T_{1:j-1}, \pi^*_{j:L}) \mid S_{1:j}, T_{1:j} \right]\\
\stackrel{(i)}{=}~& \sum_{j'=j}^L\mathbb{E}\left[ \ensuremath{Y}(T_{1:j'}, \pi_{l'+1:L})-\ensuremath{Y}(T_{1:j'-1}, \pi_{l':L}) \mid S_{1:j}, T_{1:j} \right]\\
=~& \sum_{j'=j}^L\mathbb{E}\left[ \mathbb{E}\left[\ensuremath{Y}(T_{1:j'}, \pi_{l'+1:L})-\ensuremath{Y}(T_{1:j'-1}, \pi_{l':L}) \mid X_{1:j'}, T_{1:j'} \right] \mid S_{1:j}, T_{1:j} \right]_{1:j} }\\
=~&\sum_{j'=j}^L\mathbb{E}\left[ \theta_{j'}^\top \Phi_{j'} \mid S_{1:j}, T_{1:j} \right],
\end{align*}
where (i) uses a telescoping sum.
Rearranging the above, we have
\begin{align*}
\mathbb{E}\left[ \ensuremath{Y}(T_{1:j-1}, \pi^*_{j:L}) \mid S_{1:j}, T_{1:j} \right] = \mathbb{E}\left[ \ensuremath{Y} -\sum_{j'=j}^L\theta_{j'}^\top \Phi_{j'} \mid S_{1:j}, T_{1:j} \right].
\end{align*}
\subsection{Proof of Lemma \ref{lemma:identification_parameter}}
\label{appendix:proof_identification_parameter}
We adapt the proof pattern of Lemma 7 in \cite{lewis2020double} to the RL data.
\begin{align*}
\mathbb{E}_{i,j}^+\left[R_{i,j}\, \left(\Psi_{i,j} - \bar{\Psi}_{i,j}\right)\right] =~& \mathbb{E}_{i,j}^+\left[ \mathbb{E}\left[R_{i,j} \mid \mathcal{F}_{i,j}^+, T_{i,j}\right]\, \left(\Psi_{i,j} - \bar{\Psi}_{i,j}\right)\right]_{i,j}}}
\end{align*}
Moreover, note that by Lemma~\ref{lemma:identification_policy_value} and the definition of $R_{i,j}$, we have:
\begin{align*}
\mathbb{E}\left[R_{i,j} \mid \mathcal{F}_{i,j}^+, T_{i,j}\right] = \mathbb{E}\left[\ensuremath{Y}(T_{1:j-1}, \pi^*_{j:L}) \mid \mathcal{F}_{i,j}^+, T_{i,j}\right]
\end{align*}
Moreover, by Assumption~\ref{assump:exogeneity}, we have:
\begin{align*}
\mathbb{E}\left[\ensuremath{Y}(T_{1:j-1}, \pi^*_{j:L}) \mid \mathcal{F}_{i,j}^+, T_{i,j}\right] = \mathbb{E}\left[\ensuremath{Y}(T_{1:j-1}, \pi^*_{j:L}) \mid \mathcal{F}_{i,j}^+\right] = \mathbb{E}_{i,j}^+\left[\ensuremath{Y}(T_{1:j-1}, \pi^*_{j:L})\right]
\end{align*}
Combining the last three equations, we conclude that:
\begin{align*}
\mathbb{E}_{i,j}^+\left[R_{i,j}\, \left(\Psi_{i,j} - \bar{\Psi}_{i,j}\right)\right]
=~& \mathbb{E}_{i,j}^+\left[ \mathbb{E}_{i,j}^+\left[\ensuremath{Y}(T_{1:j-1}, \pi^*_{j:L}) \right]\, \left(\Psi_{i,j} - \bar{\Psi}_{i,j}\right)\right]_{i,j}}} \\
=~& \mathbb{E}_{i,j}^+\left[\ensuremath{Y}(T_{1:j-1}, \pi^*_{j:L}) \right]\, \mathbb{E}_{i,j}^+\left[\Psi_{i,j} - \bar{\Psi}_{i,j}\right] = 0
\end{align*}
where the last equality holds, since by definition $\bar{\Psi}_{i,j}:=\mathbb{E}_{i,j}^+[\Psi_{i,j}]$.
\subsection{Proof of Lemma \ref{lemma:mds-1}}
\label{appendix:proof_mds-1}
We show that $\sum_{i=1}^n H_i\xi_{i}$ is a sum of martingale difference sequence adapted to filtration $\{\mathcal{F}_i\}$.
By Lemma~\ref{lemma:identification_parameter} and the definition of $\xi_{i,j}$, we have for any $j\in \{0,\dots,L\}$
\begin{align*}
\mathbb{E}_{i,j}^+\left[H_{i,j} \xi_{i,j} \right] = \mathbb{E}_{i,j}^+\left[ H_{i,j}\, R_{i,j}\, \left(\Psi_{i,j} - \bar{\Psi}_{i,j}\right)\right] \stackrel{(i)}{=} H_{i,j}\, \mathbb{E}_{i,j}^+\left[ R_{i,j}\, \left(\Psi_{i,j} - \bar{\Psi}_{i,j}\right)\right] \stackrel{(ii)}{=} 0
\end{align*}
where $(i)$ follows by construction that $H_{i,j}$ is measurable with respect to $\mathcal{F}^+_{i,j}$ and $(ii)$ follows by Lemma~\ref{lemma:identification_parameter}. Therefore we have:
\begin{align*}
\mathbb{E}_i\left[H_{i,j}\xi_{i,j} \right] = \mathbb{E}_i \left[\mathbb{E}_{i,j}^+\left[H_{i,j}\xi_{i,j} \right]\right]\right]}=0,\quad \forall j\in[0:L],
\end{align*}
implying $ \mathbb{E}_i\left[H_{i}\xi_{i} \right]=0$.
\subsection{Proof of Lemma \ref{lemma:mds}}
\label{appendix:proof_mds}
We first prove that $\mathbb{E}[R_{i,j}^2\mid T_{i,j}, \mathcal{F}_{i,j}^+]$ does not depend on $T_{i,j}$ such that
$\mathbb{E}[R_{i,j}^2\mid T_{i,j}, \mathcal{F}_{i,j}^+]=\mathbb{E}[R_{i,j}^2\mid \mathcal{F}_{i,j}^+]=\mathbb{E}_{i,j}^+[R_{i,j}^2]$.
Define $\Bar{R}_{i,j}:=\mathbb{E}[R_{i,j}\mid T_{i,j}, \mathcal{F}^+_{i,j}]$. We have:
\begin{align*}
\mathbb{E}[R_{i,j}^2\mid T_{i,j}, \mathcal{F}_{i,j}^+] = \mathbb{E}[(R_{i,j}-\Bar{R}_{i,j})^2\mid T_{i,j}, \mathcal{F}_{i,j}^+] + \Bar{R}_{i,j}^2 = \mathrm{Var}(R_{i,j}\mid T_{i,j}, \mathcal{F}_{i,j}^+) + \Bar{R}_{i,j}^2.
\end{align*}
By Assumption \ref{assump:homoscedasticity}, $\mathrm{Var}(R_{i,j}\mid T_{i,j}, \mathcal{F}_{i,j}^+) \equiv\mathrm{Var}(R_{i,j}\mid \mathcal{F}_{i,j}^+)$ does not depend on $T_{i,j}$. Moreover, Appendix \ref{appendix:proof_identification_parameter} proves that $\Bar{R}_{i,j}$ is independent of $T_{i,j}$ and thus $\Bar{R}_{i,j}=\mathbb{E}_{i,j}^+[R_{i,j}]$. We therefore have
\begin{equation}
\mathbb{E}[R_{i,j}^2\mid T_{i,j}, \mathcal{F}_{i,j}^+]\equiv \mathbb{E}_{i,j}^+[R_{i,j}^2], \label{eq:residual_independence}
\end{equation}
which does not depend on $T_{i,j}$.
\paragraph{Part (a).} Consider the conditional covariance. Without loss of generality, assume $j_1<j_2$. Since $H_{i,j}$ are measurable with respect to $\mathcal{F}^+_{i,j}$, we can write:
\begin{align*}
\mathbb{E}_{i,j_2}^+\left[H_{i,j_1}\xi_{i,j_1}\xi_{i,j_2}^\top H_{i,j_2}^\top \right]
=~& H_{i,j_1} \left( \Psi_{i,j_1}-\bar{\Psi}_{i,j_1} \right)
\underbrace{\mathbb{E}_{i,j_2}^+\left[ R_{i,j_1}\,R_{i,j_2} \left( \Psi_{i,j_2}-\bar{\Psi}_{i,j_2}\right)^\top \right]}_{(I)} H_{i,j_2}^\top
\end{align*}
Note that:
\begin{align*}
R_{i,j_1} = R_{i,j_2} - \sum_{j'=j_1}^{j_2-1}\Phi_{i,j'}^\top\theta_{j'}^*
\end{align*}
Thus:
\begin{align*}
(I) =~& \mathbb{E}_{i,j_2}^+\left[ R_{i,j_2}^2 \left( \Psi_{i,j_2}-\bar{\Psi}_{i,j_2}\right)^\top \right] - \sum_{j'=j_1}^{j_2-1}\Phi_{i,j'}^\top\theta_{j'}^* \mathbb{E}_{i,j_2}^+\left[ R_{i,j_2} \left( \Psi_{i,j_2}-\bar{\Psi}_{i,j_2}\right)^\top \right]\\
=~& \mathbb{E}_{i,j_2}^+\left[ R_{i,j_2}^2 \left( \Psi_{i,j_2}-\bar{\Psi}_{i,j_2}\right)^\top \right] \tag{by Lemma~\ref{lemma:identification_parameter}}\\
=~& \mathbb{E}_{i,j_2}^+\left[\mathbb{E}[ R_{i,j_2}^2 \mid T_{i,j_2}, \mathcal{F}_{i,j_2}^+]\left( \Psi_{i,j_2}-\bar{\Psi}_{i,j_2}\right)^\top \right] \\
=~& \mathbb{E}_{i,j_2}^+\left[\mathbb{E}_{i,j_2}^+[ R_{i,j_2}^2]\left( \Psi_{i,j_2}-\bar{\Psi}_{i,j_2}\right)^\top \right] \tag{by \eqref{eq:residual_independence}}\\
=~& \mathbb{E}_{i,j_2}^+\left[ R_{i,j_2}^2\right]\mathbb{E}_{i,j_2}^+\left[\left( \Psi_{i,j_2}-\bar{\Psi}_{i,j_2}\right)^\top \right] = \mathbf{0}.
\end{align*}
Therefore, we have:
\begin{equation*}
\mathrm{Cov}_i(H_{i,j_1}\xi_{i,j_1},H_{i,j_2} \xi_{i,j_2})=\mathbb{E}_i[H_{i,j_1}\xi_{i,j_1} \xi_{i,j_2}^\top H_{i,j_2}^\top] = \mathbb{E}_i\left[\mathbb{E}_{i,j_2}^+[H_{i,j_1}\xi_{i,j_1} \xi_{i,j_2}^\top H_{i,j_2}^\top]\right]= \mathbf{0}.
\end{equation*}
\paragraph{Part (b).} By the definition of $\xi_{i,j}$ and the fact that the weights $H_{i,j}$ are measurable in $\mathcal{F}_{i,j}^+$:
\begin{align*}
\mathbb{E}_{i,j}^+\left[H_{i,j_1} \xi_{i,j}\xi_{i,j}^\top H_{i,j_2}^\top\right] =~& H_{i,j}\, \mathbb{E}_{i,j}^+\left[R_{i,j}^2 \left(\Psi_{i,j} - \bar{\Psi}_{i,j}\right)\, \left(\Psi_{i,j} - \bar{\Psi}_{i,j}\right)^\top\right]\, H_{i,j}^\top\\
=~& H_{i,j}\, \mathbb{E}_{i,j}^+\left[\mathbb{E}[R_{i,j}^2\mid T_{i,j},\mathcal{F}_{i,j}^+] \left(\Psi_{i,j} - \bar{\Psi}_{i,j}\right)\, \left(\Psi_{i,j} - \bar{\Psi}_{i,j}\right)^\top\right]\, H_{i,j}^\top\\
=~& H_{i,j}\, \mathbb{E}_{i,j}^+\left[\mathbb{E}_{i,j}^+[R_{i,j}^2] \left(\Psi_{i,j} - \bar{\Psi}_{i,j}\right)\, \left(\Psi_{i,j} - \bar{\Psi}_{i,j}\right)^\top\right]\, H_{i,j}^\top \tag{by \eqref{eq:residual_independence}}\\
=~& H_{i,j}\, \mathbb{E}_{i,j}^+\left[R_{i,j}^2 \right]\, \mathbb{E}_{i,j}^+ \left[\left(\Psi_{i,j} - \bar{\Psi}_{i,j}\right)\, \left(\Psi_{i,j} - \bar{\Psi}_{i,j}\right)^\top\right]\, H_{i,j}^\top \\
=~& \mathbb{E}_{i,j}^+\left[R_{i,j}^2 \right] H_{i,j} \,\mathrm{Var}_{i,j}^+\left(\Psi_{i,j}\right)\, H_{i,j}^\top
\end{align*}
Therefore,
\begin{align*}
\mathrm{Var}_{i}(H_{i,j}\xi_{i,j})=\mathbb{E}_i[H_{i,j}\xi_{i,j}\xi_{i,j}^\top H_{i,j}^\top]=\mathbb{E}_i\left[ \mathbb{E}_{i,j}^+[H_{i,j}\xi_{i,j}\xi_{i,j}^\top H_{i,j}^\top]\right] = \mathbb{E}_i\left[\mathbb{E}^+_{i,j}\left[ R_{i,j}^2 \right] \cdot H_{i,j}\,\mathrm{Var}^+_{i,j}(\Psi_{i,j})\,
H_{i,j}^\top\right]j}^\top}.
\end{align*}
\section{Proof of Theorem \ref{thm:consistency}}
\label{appendix:consistency}
We hereby show that GMM estimator $\hat{\theta}_n\in\argmin_{\theta\in\Theta}\mathcal{L}_n(\theta)$ converges to the true parameter $\theta^*\in\argmin_{\theta\in\Theta}\mathcal{L}(\theta)$, where
\begin{align*}
\mathcal{L}_n(\theta) &= \left\|m_n(\theta)
\right\|_A^2, \quad \mbox{for}\quad m_n(\theta):=\frac{1}{n}\sum_{i=1}^n H_i\, \left(\beta_i+J_i\theta\right),\\
\mathcal{L}(\theta)&=\left\|\bar{m}_n(\theta)
\right\|_A^2, \quad \mbox{for}\quad \bar{m}_n(\theta):=\frac{1}{n}\sum_{i=1}^n H_i\left( \bar{\beta}_i+\bar{J}_i\theta\right).
\end{align*}
In Appendix \ref{appendix:thm_1_a}, we show that uniformly across $\theta\in\Theta$, we have:
\begin{equation}
\label{eq:empirical_loss_uniform_convergence}
\mathbb{E}\left[I^2\right] = O(L^2n^{\alpha_2-1}),\quad \mbox{for}\quad I:=
\max_{\theta\in\Theta}\|m_n(\theta)-\bar{m}_n(\theta)\|_A.
\end{equation}
Then for $\hat\theta_n$ that minimizes the loss $\mathcal{L}_n(\theta) = \|m_n(\theta)\|_A^2$, we have:
\begin{align*}
\|\bar{m}_n(\hat\theta_n)\|_A &= \|m_n(\hat\theta_n)\|_A + \|\bar{m}_n(\hat\theta_n)\|_A - \|m_n(\hat\theta_n)\|_A \\
&\leq \|m_n(\hat\theta_n)\|_A + \|\bar{m}_n(\hat\theta_n)-m_n(\hat\theta_n)\|_A \tag{by triangular inequality}\\
&\leq \|m_n(\theta^*)\|_A + \|\bar{m}_n(\hat\theta_n)-m_n(\hat\theta_n)\|_A \tag{$\hat\theta_n$ minimizes $\mathcal{L}_n(\theta)$}\\
&\leq \|\bar m_n(\theta^*)\|_A +\|\bar{m}_n(\theta^*)-m_n(\theta^*)\|_A+ \|\bar{m}_n(\hat\theta_n)-m_n(\hat\theta_n)\|_A \tag{by triangular inequality}\\
&=\|\bar{m}_n(\theta^*)-m_n(\theta^*)\|_A+ \|\bar{m}_n(\hat\theta_n)-m_n(\hat\theta_n)\|_A \tag{by moment equation $\bar m_n(\theta^*)=0$}\\
&\leq 2I.
\end{align*}
Therefore,
\begin{equation}
\label{eq:moment_convergence}
\mathbb{E}[ \|\bar{m}_n(\hat\theta_n)\|_A^2]\leq 4\mathbb{E}[I^2] \stackrel{\mbox{by \eqref{eq:empirical_loss_uniform_convergence}}}{=} O(L^2n^{\alpha_2-1}).
\end{equation}
On the other hand, we have:
\begin{align}
\label{eq:normalized_error_rate}
\mathbb{E}\left[
\left\|\frac{1}{n}\sum_{i=1}^n H_i \bar{J}_i(\hat{\theta}_n-\theta^*)\right\|_A^2
\right]& = \mathbb{E}\left[
\|\bar{m}_n(\hat\theta_n) - \bar{m}_n(\theta^*) \|_A^2
\right]\\
&= \mathbb{E}\left[
\|\bar{m}_n(\hat\theta_n) \|_A^2
\right]\tag{by moment equation $\bar m_n(\theta^*)=0$}\nonumber\\
& = O(L^2n^{\alpha_2-1})\tag{by \eqref{eq:empirical_loss_uniform_convergence}}.\nonumber
\end{align}
With \eqref{eq:normalized_error_rate}, Appendix \ref{appendix:thm_1_b} applies the induction method to invert $\frac{1}{n}\sum_{i=1}^n H_i \bar{J}_i$ and show that for each $j\in[0,L]$, it holds that:
\begin{align*}
\mathbb{E}\left[\left\|\hat{\theta}_{n,j}-\theta_{j}^*\right\|_2\right] =O\left(n^{\frac{(L-j+1)(\alpha_1+\alpha_2)-1}{2}}\right)\quad \mbox{and}\quad \mathbb{E}\left[\left\|\hat{\theta}_{n,j}-\theta_{j}^*\right\|^2_2\right] =O\left(n^{(2L-2j+1)\alpha_1+\alpha_2-1}\right)
\end{align*}
Concluding the proof.
\subsection{Uniform Convergence}
\label{appendix:thm_1_a}
We now show:
\begin{equation}
\mathbb{E}\left[I^2\right] = O(L^2n^{\alpha_2-1}),\quad \mbox{for}\quad I:=
\max_{\theta\in\Theta}\|m_n(\theta)-\bar{m}_n(\theta)\|_A.\tag{\ref{eq:empirical_loss_uniform_convergence}}
\end{equation}
We have:
\begin{align*}
I&=
\max_{\theta\in\Theta}\|m_n(\theta)-\bar{m}_n(\theta)\|_A\\
&=\max_{\theta\in\Theta}\left\|
\frac{1}{n}\sum_{i=1}^nH_i(\beta_i-\bar \beta_i) + \frac{1}{n}\sum_{i=1}^nH_i(J_i-\bar J_i) \theta
\right\|_A\\
&\lesssim\max_{\theta\in\Theta}\left\|
\frac{1}{n}\sum_{i=1}^nH_i(\beta_i-\bar \beta_i) + \frac{1}{n}\sum_{i=1}^nH_i(J_i-\bar J_i) \theta
\right\|_2 \tag{$\lambda_{\max}(A)=O(1)$}\\
&\leq \max_{\theta\in\Theta} \left\|
\frac{1}{n}\sum_{i=1}^nH_i(\beta_i-\bar \beta_i)\right\|_2 + \left\|\frac{1}{n}\sum_{i=1}^nH_i(J_i-\bar J_i) \theta
\right\|_2\tag{triangular inequality}\\
&\lesssim \underbrace{\left\|
\frac{1}{n}\sum_{i=1}^nH_i(\beta_i-\bar \beta_i)\right\|_2}_{C} + \underbrace{\left\|\frac{1}{n}\sum_{i=1}^nH_i(J_i-\bar J_i)
\right\|_2}_{D},
\end{align*}
where for a symmatric matrix $M$, we use $\lambda_{\max}(M)$ to denote its largest eigenvalue and use $\|M\|_2$ to denote its operator norm induced by vector 2-norm.
Appendix \ref{appendix:term_c} shows that $\mathbb{E}[C^2]=O(Ln^{\alpha_2-1})$, and Appendix \ref{appendix:term_d} shows that $\mathbb{E}[D^2]=O(L^2n^{\alpha_2-1})$. Therefore,
\[
\mathbb{E}[I^2]\lesssim \mathbb{E}[C^2]+ \mathbb{E}[D^2] = O(L^2n^{\alpha_2-1}),
\]
proving \eqref{eq:empirical_loss_uniform_convergence}.
\subsubsection{Bounding term C.}
\label{appendix:term_c}
We show:
\begin{equation}
\label{eq:beta_convergence}
\mathbb{E}\left[C^2\right]=\mathbb{E}\left[ \left\|\frac{1}{n}\sum_{i=1}^n H_i\, \beta_i - \frac{1}{n}\sum_{i=1}^n H_i\bar{\beta}_i\right\|_2^2\right] = O(Ln^{\alpha_2-1}).
\end{equation}
Note that $\{H_i(\beta_i-\bar \beta_i)\}$ is a martingale difference sequence adapted to $\{\mathcal{F}_i\}$.
This is because each $j$-th component
\[
\mathbb{E}_i[H_{i,j}(\beta_{i,j}-\bar \beta_{i,j})] = \mathbb{E}_i \mathbb{E}_{i,j}^+[H_{i,j}(\beta_{i,j}-\bar \beta_{i,j})] = \mathbb{E}_i \left[H_{i,j}\mathbb{E}_{i,j}^+[(\beta_{i,j}-\bar \beta_{i,j})]\right]=0,
\]
by the definition $\bar \beta_{i,j}=\mathbb{E}_{i,j}^+[\beta_{i,j}]$ and $H_{i,j}$ adapted to $\mathcal{F}_{i,j}^+$.
We have
\begin{align*}
\mathbb{E}\left[ \left\|\frac{1}{n}\sum_{i=1}^n H_i\, \beta_i - \frac{1}{n}\sum_{i=1}^n H_i\bar{\beta}_i\right\|_2^2 \right] &\stackrel{(i)}{=} \frac{1}{n^2}\sum_{i=1}^n\mathbb{E}\left[ (\beta_i - \bar{\beta}_i)^\top H_i^\top H_i(\beta_i - \bar{\beta}_i)\right]\\
& = \sum_{j=0}^L\frac{1}{n^2}\sum_{i=1}^n \mathbb{E}\left[ (\beta_{i,j} - \bar{\beta}_{i,j})^\top H_{i,j}^\top H_{i,j}(\beta_{i,j} - \bar{\beta}_{i,j})\right]\\
& = \sum_{j=0}^L\frac{1}{n^2}\sum_{i=1}^n \mathrm{Tr}\left(\mathbb{E}\left[ H_{i,j}(\beta_{i,j} - \bar{\beta}_{i,j})(\beta_{i,j} - \bar{\beta}_{i,j})^\top H_{i,j}^\top \right]\right),
\end{align*}
where (i) is due to that $\{H_i(\beta_i-\bar \beta_i)\}$ is a martingale difference sequence.
Now consider each diagonal entry $j$. We have:
\begin{align*}
&\frac{1}{n^2}\sum_{i=1}^n\mathrm{Tr}\left(\mathbb{E}\left[ H_{i,j}(\beta_{i,j} - \bar{\beta}_{i,j})(\beta_{i,j} - \bar{\beta}_{i,j})^\top H_{i,j}^\top\right]\right)\\
= & \frac{1}{n^2}\sum_{i=1}^n\mathrm{Tr}\left(\mathbb{E}\left[ H_{i,j}\mathrm{Var}_{i,j}^+(Y_i(\Psi_{i,j}-\bar{\Psi}_{i,j})) H_{i,j}^\top\right] \right) = O(n^{-2}) \sum_{i=1}^n\mathrm{Tr}\left(\mathbb{E}\left[ H_{i,j}\mathrm{Var}_{i,j}^+(\Psi_{i,j}) H_{i,j}^\top\right]\right) = O(n^{\alpha_2-1}),
\end{align*}
where the last equality is by Property \ref{property:weight_regularizing}(c). We thus have \eqref{eq:beta_convergence} hold.
\subsubsection{Bounding term D.}
\label{appendix:term_d}
We next show:
\begin{equation}
\label{eq:coeff_convergence}
\mathbb{E}\left[ \left\| \frac{1}{n}\sum_{i=1}^n H_i\, J_i - \frac{1}{n}\sum_{i=1}^n H_i\bar{ J}_i\right\|_2^2 \right]= O(L^2n^{\alpha_2-1}).
\end{equation}
Similarly, $\{H_i(J_i-\bar J_i)\}$ is a martingale difference sequence adapted to $\{\mathcal{F}_i\}$.
For any given $\theta\in\Theta$ with $\Theta$ being bounded, we have
\begin{align}
&\mathbb{E}\left[ \left\|\frac{1}{n}\sum_{i=1}^n H_i\, (J_i -\bar{J}_i)\theta\right\|_2^2 \right] \leq \mathbb{E}\left[ \left\|\frac{1}{n}\sum_{i=1}^n H_i\, (J_i -\bar{J}_i)\right\|_{2}^2 \cdot \|\theta^*\|_2^2 \right]\nonumber\\
\stackrel{(i)}{\lesssim} & \mathbb{E}\left[ \left\|\frac{1}{n}\sum_{i=1}^n H_i\, (J_i -\bar{J}_i)\right\|_{Frob}^2 \right] =\mathbb{E}\left[\mathrm{Tr}\left(\left(\frac{1}{n}\sum_{i=1}^n H_i\, (J_i -\bar{J}_i)\right) \cdot \left(\frac{1}{n}\sum_{i=1}^n H_i\, (J_i -\bar{J}_i)\right)^\top\right)_i)\right)^\top}\right]\nonumber\\
= &\mathrm{Tr}\left(\mathbb{E}\left[\left(\frac{1}{n}\sum_{i=1}^n H_i\, (J_i -\bar{J}_i)\right) \cdot \left(\frac{1}{n}\sum_{i=1}^n H_i\, (J_i -\bar{J}_i)\right)^\top\right]\right)ht)^\top\right]}\nonumber\\
\stackrel{(ii)}{=}& \frac{1}{n^2}\sum_{i=1}^n \mathrm{Tr}\left(\mathbb{E}\left[H_i\, (J_i -\bar{J}_i)(J_i -\bar{J}_i)^\top H_i^\top \right]\right)\nonumber \\
\stackrel{(iii)}{=}&\frac{1}{n^2}\sum_{i=1}^n \sum_{j=0}^L \sum_{k\geq j}\mathrm{Tr}\left(\mathbb{E}\left[H_{i,j}\, (J_{i,j,k} -\bar{J}_{i,j,k})(J_{i,j,k} -\bar{J}_{i,j,k})^\top H_{i,j}^\top \right]\right):=\frac{1}{n^2}\sum_{i=1}^n \sum_{j=0}^L \sum_{k\geq j}\mathbb{E}\left[\Lambda_{i,j,k}\right]\label{eq:jacobian_jk}
\end{align}
where $\Lambda_{i,j,k}:=\mathrm{Tr}\left(\mathbb{E}_{i,j}^+\left[H_{i,j}\, (J_{i,j,k} -\bar{J}_{i,j,k})(J_{i,j,k} -\bar{J}_{i,j,k})^\top H_{i,j}^\top \right]\right)$; (i) uses the boundedness of $\theta$ and that a matrix $\ell_2$ norm is bounded by its Frobenius norm; (ii) uses that $\{H_i\, (J_i -\bar{J}_i)\}_{i=1}^n$ is a martingale difference sequence; (iii) is due to that $H_i$ is a block-diagonal matrix and $J_i$ is a block upper-triangular matrix.
Now let's look into \eqref{eq:jacobian_jk} and analyze $\Lambda_{i,j,k}$.
\begin{align*}
\Lambda_{i,j,k}:=&\mathrm{Tr}\left(\mathbb{E}_{i,j}^+\left[H_{i,j}\, (J_{i,j,k} -\bar{J}_{i,j,k})(J_{i,j,k} -\bar{J}_{i,j,k})^\top H_{i,j}^\top \right]\right)\\
= & \mathrm{Tr}( H_{i,j}\mathbb{E}_{i,j}^+[ ((\Psi_{i,j} -\bar{\Psi}_{i,j})\Phi_{i,k}^\top - \mathbb{E}_{i,j}^+[(\Psi_{i,j} -\bar{\Psi}_{i,j})\Phi_{i,k}^\top])((\Psi_{i,j} -\bar{\Psi}_{i,j})\Phi_{i,k}^\top - \mathbb{E}_{i,j}^+[(\Psi_{i,j} -\bar{\Psi}_{i,j})\Phi_{i,k}^\top])^\top ]H_{i,j}^\top)\\
\stackrel{(i)}{\lesssim} & \mathrm{Tr}(H_{i,j}\mathrm{Var}_{i,j}^+(\Psi_{i,j} -\bar{\Psi}_{i,j})H_{i,j}^\top),
\end{align*}
where (i) is due to Lemma \ref{lemma:var_inequality}, boundedness of $\|\Phi_{i,k}\|_2$, and that if positive semi-definite matrices $A\preceq B$, then $\mathrm{Tr}(CAC^\top)\leq \mathrm{Tr}(CBC^\top)$ for any matrix $C$.
Continuing \eqref{eq:jacobian_jk}, we have:
\begin{align*}
\mathbb{E}\left[ \left\|\frac{1}{n}\sum_{i=1}^n H_i\, (J_i -\bar{J}_i)\theta\right\|_2^2 \right] &\lesssim \frac{1}{n^2}\sum_{i=1}^n \sum_{j=0}^L \sum_{k\geq j}\mathbb{E}\left[\mathrm{Tr}(H_{i,j}\mathrm{Var}_{i,j}^+(\Psi_{i,j} -\bar{\Psi}_{i,j})H_{i,j}^\top)\right]\stackrel{(i)}{=}O(L^2n^{\alpha_2-1}),
\end{align*}
where (i) is by Property \ref{property:weight_regularizing}(c). Thus we have \eqref{eq:coeff_convergence} hold.
\subsection{Induction Step}
\label{appendix:thm_1_b}
We now bound the estimation error $(\hat\theta_n - \theta^*)$. Equation
\eqref{eq:normalized_error_rate} shows that the normalized error satisfies
$ \mathbb{E}[\|\bar B_n (\hat{\theta}_n-\theta^*)
\|_{2}^2]=O(n^{\alpha_2-1})$,
for $ \Bar B_{n} := -n^{-1}\sum_{i=1}^nH_i \bar{J}_i$.
Note that $\Bar{B}_n$ is an upper block-triangular matrix. We will show that, for each diagonal block $\Bar{B}_{n,j,j}$, the matrix $\Bar{B}_{n,j,j}^\top \Bar{B}_{n,j,j}$ is invertible with high probability. This allows us to invert $\Bar{B}_{n,j,j}^\top \Bar{B}_{n,j,j}$ to bound the estimation error $(\hat{\theta}_{n,j} - \theta_j^*)$ in a backwards manner.
In particular, we apply induction method to prove the result. We start by introducing a few notations. With $\Bar B_{n} = -n^{-1}\sum_{i=1}^nH_i \bar{J}_i$, define $\Bar B_{n,j,j'}=n^{-1}\sum_{i=1}^n H_{i,j}\mathbb{E}_{i,j}^+\left[(\Psi_{i,j}-\bar{\Psi}_{i,j})\Phi_{i,j'}^\top\right]$ as its $(j,j')$ block.
Similarly, define $ B^0_{n,j,j'}=n^{-1}\sum_{i=1}^n \mathbb{E}_{i,j}\left[H_{i,j}(\Psi_{i,j}-\bar{\Psi}_{i,j})\Phi_{i,j'}^\top\right]$. Define their difference $\delta_{n,j,j'}:=B^0_{n,j,j'}-\Bar B_{n,j,j'}$. Appendix \ref{appendix:regularity_Bn} shows that:
\begin{equation}
\label{eq:bn_regularity}
\|\Bar B_{n,j,j'}\|^2_{Frob}=O\left(n^{\alpha_1}\right),\quad \|B^0_{n,j,j'}\|^2_{Frob}=O\left(n^{\alpha_2}\right), \quad \mathbb{E}\left[\left\|\delta_{n,j,j'}\right\|_{Frob}^2\right]=O\left(n^{\alpha_2-1}\right)
\end{equation}
By Property \ref{property:weight_regularizing}(a), we have $\forall l\in [0:L]$, $(B^0_{n,j,j})^\top B^0_{n,j,j} \succeq c_1^2 n^{-\alpha_1}I$, where $c_1$ is a constant introduced in Property \ref{property:weight_regularizing}. The regularity of $\bar B_{n,j,j}$ inherits from that of $B^0_{n,j,j}$.
Define the event $\mathcal{E}_n$ as:
\begin{align*}
\mathcal{E}_n = \left\{\Bar B_{n,j,j}^\top \Bar B_{n,j,j} \succeq \frac{c_1^2}{4} n^{-\alpha_1}I, \ \forall l\in\{0,\dots, L\}\right\},
\end{align*}
Appendix \ref{appendix:regularity_Bn} shows that $\mathcal{E}_n $ happens with high probability:
\begin{align}
\label{eq:en_high_probability}
&1-\mathbb{P}(\mathcal{E}_n)=O(Ln^{\alpha_1+\alpha_2-1}).
\end{align}
We now apply the induction step to bound the estimation error $\hat\theta_n - \theta^*$.
\paragraph{Base case: $j=L$.} We have:
\begin{align*}
\mathbb{E}\left[\|\hat{\theta}_{n,L}-\theta^*_L\|_2^2\right]& \lesssim n^{\alpha_1}\mathbb{E}\left[\left\|\Bar B_{n,L,L}(\hat{\theta}_{n,L}-\theta_L^*)\right\|^2_2\mathbf{1}(\mathcal{E}_n) \right] + \mathbb{P}(\mathcal{E}_n^c)\\
&\leq n^{\alpha_1}\mathbb{E}\left[\left\|\Bar B_{n,L,L}(\hat{\theta}_{n,L}-\theta_L^*)\right\|^2_2 \right] + \mathbb{P}(\mathcal{E}_n^c)\\
&\stackrel{(i)}{\leq} n^{\alpha_1}\mathbb{E}\left[\left\|\Bar B_{n}(\hat{\theta}_{n}-\theta^*)\right\|^2_2\right] + \mathbb{P}(\mathcal{E}_n^c)\\
&= n^{\alpha_1} \mathbb{E}\left[\left\|\frac{1}{n}\sum_{i=1}^n H_i\bar{J}_i \left(\hat{\theta}_n-\theta^*\right)
\right\|_{2}^2\right] + \mathbb{P}(\mathcal{E}_n^c)\stackrel{(ii)}{=} O(n^{\alpha_1+\alpha_2-1}),
\end{align*}
where (i) is because $B_n$ is a block upper triangular matrix and
(ii) is by \eqref{eq:normalized_error_rate} and \eqref{eq:en_high_probability}. This result also yield the $L^1$-norm convergence:
\begin{equation*}
\mathbb{E}\left[\|\hat{\theta}_{n,L}-\theta^*_L\|_2\right]\leq \mathbb{E}\left[\|\hat{\theta}_{n,L}-\theta^*_L\|_2^2\right]^{1/2} = O(n^{\frac{\alpha_1+\alpha_2-1}{2}}).
\end{equation*}
\paragraph{Induction step.} Now consider $j<L$ recursively. Assume the induction hypothesis holds that for $j'=j+1, \dots, L$, such that:
\begin{align}
\label{eq:induction}
\mathbb{E}\left[\left\|\hat{\theta}_{n,j'}-\theta_{j'}^*\right\|_2\right] =O\left(n^{\frac{(L-j'+1)(\alpha_1+\alpha_2)-1}{2}}\right)\quad \mbox{and}\quad \mathbb{E}\left[\left\|\hat{\theta}_{n,j'}-\theta_{j'}^*\right\|^2_2\right] =O\left(n^{(2L-2j'+1)\alpha_1+\alpha_2-1}\right).
\end{align}
Let's first prove the $L^2$-norm convergence. We have
\begin{align*}
& \mathbb{E}\left[\|\hat{\theta}_{n,j}-\theta^*_j\|_2^2\right] \lesssim n^{\alpha_1}\mathbb{E}\left[\left\|\Bar B_{n,j,j}(\hat{\theta}_{n,j}-\theta^*_j)\right\|^2_2\mathbf{1}(\mathcal{E}_n) \right] + \mathbb{P}(\mathcal{E}_n^c)\leq n^{\alpha_1}\mathbb{E}\left[\left\|\Bar B_{n,j,j}(\hat{\theta}_{n,j}-\theta^*_j)\right\|^2_2\right] + \mathbb{P}(\mathcal{E}_n^c)\\
&= n^{\alpha_1}\mathbb{E}\left[\left\|\Bar B_{n,j,j}(\hat{\theta}_{n,j}-\theta^*_j) + \sum_{j'>j}\Bar B_{n,j,j'}(\hat{\theta}_{n,j'}-\theta^*_{j'})
- \sum_{j'>j}\Bar B_{n,j,j'}(\hat{\theta}_{n,j'}-\theta^*_{j'})
\right\|^2_2 \right] + \mathbb{P}(\mathcal{E}_n^c)\\
&\stackrel{(i)}{\leq} 2L n^{\alpha_1}\mathbb{E}\left[\left\|\Bar B_{n,j,j}(\hat{\theta}_{n,j}-\theta^*_j) + \sum_{j'>j}\Bar B_{n,j,j'}(\hat{\theta}_{n,j'}-\theta^*_{j'}) \right\|^2_2 \right] + 2L
n^{\alpha_1}
\sum_{j'>j} \mathbb{E}\left[\left\|\Bar B_{n,j,j'}(\hat{\theta}_{n,j'}-\theta^*_{j'})\right\|^2_2 \right]+ \mathbb{P}(\mathcal{E}_n^c)\\
&\stackrel{(ii)}{\leq} 2L n^{\alpha_1}\mathbb{E}\left[\left\|\Bar B_{n}(\hat{\theta}_{n}-\theta^*)\right\|^2_2 \right]+
2L
n^{\alpha_1}
\sum_{j'>j} \mathbb{E}\left[\left\|\Bar B_{n,j,j'}(\hat{\theta}_{n,j'}-\theta^*_{j'})\right\|^2_2 \right]+ \mathbb{P}(\mathcal{E}_n^c)\\
&= 2L n^{\alpha_1} \mathbb{E}\left[\left\|\frac{1}{n}\sum_{i=1}^n H_i\bar{J}_i \left(\hat{\theta}_n-\theta^*\right)
\right\|_{2}^2\right] +
2L
n^{\alpha_1}
\sum_{j'>j} \mathbb{E}\left[\left\|\Bar B_{n,j,j'}(\hat{\theta}_{n,j'}-\theta^*_{j'})\right\|^2_2 \right]+ \mathbb{P}(\mathcal{E}_n^c)\\
&\leq 2L n^{\alpha_1} \mathbb{E}\left[\left\|\frac{1}{n}\sum_{i=1}^n H_i\bar{J}_i \left(\hat{\theta}_n-\theta^*\right)
\right\|_{2}^2\right] +
2L
n^{\alpha_1}
\sum_{j'>j} \mathbb{E}\left[\left\|\Bar B_{n,j,j'}\right\|_{Frob}^2\left\|(\hat{\theta}_{n,j'}-\theta^*_{j'})\right\|^2_2 \right]+ \mathbb{P}(\mathcal{E}_n^c)\\
&\stackrel{(iii)}{=} 2LO(n^{\alpha_2-1}) + 2L\sum_{j'>j} O\left(n^{\alpha_1+\alpha_1+(2L-2j'+1)\alpha_1+\alpha_2-1}\right) + O(n^{\alpha_1+\alpha_2-1}) = O\left(n^{(2L-2j+1)\alpha_1+\alpha_2-1}\right),
\end{align*}
where (i) is by triangular inequality and Cauchy-Swarchz inequality, (ii) is because $\Bar B_n$ is a block upper triangular matrix and (iii) is by \eqref{eq:normalized_error_rate}, \eqref{eq:bn_regularity}, \eqref{eq:induction}, and \eqref{eq:en_high_probability}.
We then look into the $L^1$-norm convergence by having:
\begin{align*}
& \mathbb{E}\left[\|\hat{\theta}_{n,j}-\theta^*_j\|_2\right] \lesssim n^{\frac{\alpha_1}{2}}\mathbb{E}\left[\left\|\Bar B_{n,j,j}(\hat{\theta}_{n,j}-\theta^*_j)\right\|_2\mathbf{1}(\mathcal{E}_n) \right] + \mathbb{P}(\mathcal{E}_n^c)\leq n^{\frac{\alpha_1}{2}}\mathbb{E}\left[\left\|\Bar B_{n,j,j}(\hat{\theta}_{n,j}-\theta^*_j)\right\|_2\right] + \mathbb{P}(\mathcal{E}_n^c)\\
&\leq n^{\frac{\alpha_1}{2}}\mathbb{E}[\|\Bar B_{n,j,j}(\hat{\theta}_{n,j}-\theta^*_j) + \sum_{j'>j}\Bar B_{n,j,j'}(\hat{\theta}_{n,j'}-\theta^*_{j'}) \|_2 ] +
n^{\frac{\alpha_1}{2}}
\sum_{j'>j} \mathbb{E}\left[\left\|\Bar B_{n,j,j'}(\hat{\theta}_{n,j'}-\theta^*_{j'})\right\|_2 \right]+ \mathbb{P}(\mathcal{E}_n^c)\\
&\stackrel{(i)}{\leq} n^{\frac{\alpha_1}{2}}\mathbb{E}\left[\left\|\Bar B_{n}(\hat{\theta}_{n}-\theta^*)\right\|_2 \right]+
n^{\frac{\alpha_1}{2}}
\sum_{j'>j} \mathbb{E}\left[\left\|\Bar B_{n,j,j'}(\hat{\theta}_{n,j'}-\theta^*_{j'})\right\|_2 \right]+ \mathbb{P}(\mathcal{E}_n^c)\\
&= n^{\frac{\alpha_1}{2}} \mathbb{E}\left[\left\|\frac{1}{n}\sum_{i=1}^n H_i\bar{J}_i \left(\hat{\theta}_n-\theta^*\right)
\right\|_{2}\right] +
n^{\frac{\alpha_1}{2}}
\sum_{j'>j} \mathbb{E}\left[\left\|\Bar B_{n,j,j'}(\hat{\theta}_{n,j'}-\theta^*_{j'})\right\|^2_2 \right]+ \mathbb{P}(\mathcal{E}_n^c)\\
&\stackrel{(ii)}{\leq} O(n^{\frac{\alpha_1+\alpha_2-1}{2}}) +
n^{\frac{\alpha_1}{2}}
\sum_{j'>j} \mathbb{E}\left[\left\|B^0_{n,j,j'}(\hat{\theta}_{n,j'}-\theta^*_{j'})\right\|_2 \right]+
n^{\frac{\alpha_1}{2}}
\sum_{j'>j} \mathbb{E}\left[\left\|\delta_{n,j,j'}(\hat{\theta}_{n,j'}-\theta^*_{j'})\right\|_2 \right] + O(n^{\alpha_1+\alpha_2-1}) \\
&\stackrel{(iii)}{\leq} O(n^{\frac{\alpha_1+\alpha_2-1}{2}}) +
n^{\frac{\alpha_1}{2}}
\sum_{j'>j} \mathbb{E}\left[\left\|B^0_{n,j,j'}\right\|_{Frob}\left\|(\hat{\theta}_{n,j'}-\theta^*_{j'})\right\|_2 \right]+
n^{\frac{\alpha_1}{2}}
\sum_{j'>j} \mathbb{E}\left[\left\|\delta_{n,j,j'}\right\|^2_{Frob}\right]^{\frac{1}{2}}\mathbb{E}\left[\left\|(\hat{\theta}_{n,j'}-\theta^*_{j'})\right\|^2_2 \right]^{\frac{1}{2}} \\
&\stackrel{(iv)}{=} O(n^{\frac{\alpha_1+\alpha_2-1}{2}}) + n^{\frac{\alpha_1}{2}}\sum_{j'>j} O(n^{\frac{\alpha_2}{2}+\frac{(L-j'+1)(\alpha_1+\alpha_2)-1}{2}}) +
n^{\frac{\alpha_1}{2}}
\sum_{j'>j} O(n^{\frac{\alpha_2-1}{2} + \frac{(2L-2j'+1)\alpha_1+\alpha_2-1}{2}})\\
&= O(n^{\frac{\alpha_1+\alpha_2-1}{2}}) + O(n^{\frac{(L-j+1)(\alpha_1+\alpha_2)-1}{2}}) + O(n^{\frac{(2L-2j)\alpha_1+2\alpha_2-2}{2}})
\\
&\stackrel{(v)}{=}O(n^{\frac{(L-j+1)(\alpha_1+\alpha_2)-1}{2}}).
\end{align*}
where (i) uses that $\bar B_n$ is a block upper triangular matrix, (ii) uses \eqref{eq:normalized_error_rate}, triangular inequality, and \eqref{eq:en_high_probability}, (iii) uses Cauchy-Schwartz inequality and that a matrix $\ell_2$ norm is bounded by its Frobenius norm, (iv) uses the induction assumption \eqref{eq:induction}, and (v) uses that $\alpha_2\leq \alpha_1\leq 1/L$ as specified in Property \ref{property:weight_regularizing}.
\subsubsection{Asymptotic regularity of $\Bar{B}_n$.}
\label{appendix:regularity_Bn}
Note that $\Bar{B}_n$ is block upper triangular, where the column sizes of blocks correspond to the decomposition of $\theta=(\theta_0, \theta_1, \dots, \theta_L)$. For any $j\leq j'$ we have:
\begin{align*}
\Bar B_{n,j,j'}=n^{-1}\sum_{i=1}^n H_{i,j}\mathbb{E}_{i,j}^+\left[(\Psi_{i,j}-\bar{\Psi}_{i,j})\Phi_{i,j'}^\top\right] .
\end{align*}
Define the matrix $B_n^0$ as follows: for each of its $(j,j')$-th block with $j\leq j'$, define
\begin{align*}
B^0_{n,j,j'}=n^{-1}\sum_{i=1}^n \mathbb{E}_{i,j}\left[H_{i,j}(\Psi_{i,j}-\bar{\Psi}_{i,j})\Phi_{i,j'}^\top\right] ,\mbox{ and }\delta_{n,j,j'}=B^0_{n,j,j'}-\Bar B_{n,j,j'}.
\end{align*}
First, notice that the singular values of the $j$-th diagonal block $B^0_{n,j,j}$ are lower-bounded when Property \ref{property:weight_regularizing}(a) holds.
For the $0$-th diagonal block corresponding to $\theta_0^*$, we have $B^0_{n,0,0}= n^{-1}\sum_{i=1}^n \mathbb{E}_{i,0}[H_{i,0}] =1\geq c_1 n^{-\frac{\alpha_1}{2}}$ for $H_{i,0}=1$ and large $n$.
Next we shall show that for $j\leq j'$, we have
\begin{align*}
\mathbb{E}\left[\left\|\delta_{n,j,j'}\right\|_{Frob}^2\right]=O\left(n^{\alpha_2-1}\right), \, \|B^0_{n,j,j'}\|^2_{Frob}=O\left(n^{\alpha_2}\right), \, \|\Bar B_{n,j,j'}\|^2_{Frob}=O\left(n^{\alpha_1}\right),
\end{align*}
and that the event $\mathcal{E}_n$ happens with high probability:
\begin{align*}
1 - \mathbb{P}(\mathcal{E}_n)=O(Ln^{\alpha_1+\alpha_2-1}).
\end{align*}
\paragraph{Asymptotic neglibility of $\delta_{n,j,j'}$.} For a matrix $A$, we use $\mathrm{Tr}(A)$ to denote its trace, i.e. the sum of its diagonal entries. With $j\leq j'$, by the definition of the Frobenius norm, we have
\begin{align*}
\mathbb{E}[\|\delta_{n,j,j'}\|_{Frob}^2 ]=~& \mathbb{E}\left[\mathrm{Tr}\left(\delta_{n,j,j'}\, \delta_{n,j,j'}^\top\right)\right]\\
=~& n^{-2}\mathbb{E}\Big[\mathrm{Tr}\Big(\left\{\sum_{i=1}^n \left(H_{i,j}\mathbb{E}_{i,j}^+[(\Psi_{i,j}-\bar{\Psi}_{i,j})\Phi_{i,j'}^\top] -\mathbb{E}_{i,j}[H_{i,j}(\Psi_{i,j}-\bar{\Psi}_{i,j})\Phi_{i,j'}^\top ]\right)\right\}\\
&\quad\quad \times\left\{\sum_{i=1}^n \left(H_{i,j}\mathbb{E}_{i,j}^+[(\Psi_{i,j}-\bar{\Psi}_{i,j})\Phi_{i,j'}^\top] -\mathbb{E}_{i,j}[H_{i,j}(\Psi_{i,j}-\bar{\Psi}_{i,j})\Phi_{i,j'}^\top ]\right)\right\}^\top \Big)\Big]\\
\stackrel{(i)}{=}&
n^{-2}\sum_{i=1}^n\mathrm{Tr}\Big( \mathbb{E}\Big[\left(H_{i,j}\mathbb{E}_{i,j}^+[(\Psi_{i,j}-\bar{\Psi}_{i,j})\Phi_{i,j'}^\top] -\mathbb{E}_{i,j}[H_{i,j}(\Psi_{i,j}-\bar{\Psi}_{i,j})\Phi_{i,j'}^\top ]\right)\\
&\quad \quad \times \left(H_{i,j}\mathbb{E}_{i,j}^+[(\Psi_{i,j}-\bar{\Psi}_{i,j})\Phi_{i,j'}^\top] -\mathbb{E}_{i,j}[H_{i,j}(\Psi_{i,j}-\bar{\Psi}_{i,j})\Phi_{i,j'}^\top ]\right)^\top
\Big]\Big)\\
\stackrel{(ii)}{\leq} & n^{-2}\sum_{i=1}^n\mathrm{Tr} \left(
\mathbb{E}\left[ H_{i,j} \mathbb{E}^+_{i,j}[(\Psi_{i,j}-\bar{\Psi}_{i,j})\Phi_{i,j'}^\top] \mathbb{E}^+_{i,j}[\Phi_{i,j'}(\Psi_{i,j}-\bar{\Psi}_{i,j})^\top] H_{i,j}^\top \right]\right)\\
\stackrel{(iii)}{\lesssim} & n^{-2}d\sum_{i=1}^n\mathrm{Tr} \left(
\mathbb{E}\left[ H_{i,j} \mathbb{E}^+_{i,j}\left[(\Psi_{i,j}-\bar{\Psi}_{i,j})(\Psi_{i,j}-\bar{\Psi}_{i,j})^\top \right] H_{i,j}^\top \right]\right)^\top }
} \stackrel{(iv)}{=} O(n^{\alpha_2-1})
\end{align*}
where (i) uses the fact that $\delta_{n,j,j'}$ is a sum of martingale difference sequence, (ii) uses Lemmas \ref{lemma:trace_inequality_3} \& \ref{lemma:trace_inequality_4}, (iii) uses that $\|\Phi_{i,j'}\|_2$ is bounded, Lemma \ref{lemma:matrix_inequality} and Lemma \ref{lemma:trace_inequality_3}, and (iv) uses Property \ref{property:weight_regularizing}(c).
\paragraph{Uniform bound of $B^0_{n,j,j'}$.} With $j\leq j'$, we have
\begin{align*}
\|B^0_{n,j,j'}\|_{Frob}^2=&n^{-2}\left\|
\sum_{i=1}^n \mathbb{E}_{i,j}[H_{i,j}(\Psi_{i,j}-\bar{\Psi}_{i,j})\Phi_{i,j'}^\top ]\right\|^2_{Frob}
\\
\stackrel{(i)}{\leq} &n^{-2}
\left(\sum_{i=1}^n \left\|\mathbb{E}_{i,j}[H_{i,j}(\Psi_{i,j}-\bar{\Psi}_{i,j})\Phi_{i,j'}^\top ]\right\|_{Frob}\right)^2
\\
\stackrel{(ii)}{\leq} &n^{-1}
\sum_{i=1}^n \left\|\mathbb{E}_{i,j}[H_{i,j}(\Psi_{i,j}-\bar{\Psi}_{i,j})\Phi_{i,j'}^\top ]\right\|^2_{Frob}
\\
= &n^{-1}
\sum_{i=1}^n \mathrm{Tr}\left(\mathbb{E}_{i,j}[H_{i,j}(\Psi_{i,j}-\bar{\Psi}_{i,j})\Phi_{i,j'}^\top ]\mathbb{E}_{i,j}[\Phi_{i,j'}(\Psi_{i,j}-\bar{\Psi}_{i,j})^\top H_{i,j}^\top] \right)
\\
\stackrel{(iii)}{\lesssim} &n^{-1}
\sum_{i=1}^n \mathrm{Tr}\left(\mathbb{E}_{i,j}[H_{i,j}(\Psi_{i,j}-\bar{\Psi}_{i,j})(\Psi_{i,j}-\bar{\Psi}_{i,j})^\top H_{i,j}^\top] \right)
\\
= &
n^{-1}\sum_{i=1}^n \mathrm{Tr}\left(\mathbb{E}_{i,j}[H_{i,j}\mathrm{Var}_{i,j}^+(\Psi_{i,j})
H_{i,j}^\top] \right) = O(n^{\alpha_2}).
\end{align*}
where (i) is by triangular inequality, (ii) is by Cauchy-Schwartz inequality, (iii) is by the assumption that $\|\Phi_{i,j'}\|_2$ is bounded, Lemma \ref{lemma:matrix_inequality} and Lemma \ref{lemma:trace_inequality_2}.
\paragraph{Uniform bound of $\Bar B_{n,j,j'}$.} With $j\leq j'$, we have
\begin{align*}
\|\Bar B_{n,j,j'}\|_{Frob}^2=&n^{-2}\left\|
\sum_{i=1}^n H_{i,j}\mathbb{E}_{i,j}^+[(\Psi_{i,j}-\bar{\Psi}_{i,j})\Phi_{i,j'}^\top] \right\|^2_{Frob}
\stackrel{(i)}{\leq} n^{-1}
\sum_{i=1}^n \left\|H_{i,j}\mathbb{E}_{i,j}^+[(\Psi_{i,j}-\bar{\Psi}_{i,j})\Phi_{i,j'}^\top] \right\|^2_{Frob}\\
\leq &n^{-1}
\sum_{i=1}^n \mathrm{Tr}\left(H_{i,j}\mathbb{E}_{i,j}^+[(\Psi_{i,j}-\bar{\Psi}_{i,j})\Phi_{i,j'}^\top]\mathbb{E}_{i,j}^+[ \Phi_{i,j'}(\Psi_{i,j}-\bar{\Psi}_{i,j})^\top] H_{i,j}^\top \right)
\\
\stackrel{(ii)}{\lesssim} &n^{-1}
\sum_{i=1}^n \mathrm{Tr}\left(H_{i,j}H_{i,j}^\top \right) = O(n^{\alpha_1}).
\end{align*}
where (i) is by triangular inequality and Cauchy-Schwartz inequality, (ii) is due to the fact that $\Phi_{i,j'}, \Psi_{i,j},\bar{\Psi}_{i,j}$ are bounded such that
\[
H_{i,j}\mathbb{E}_{i,j}^+[(\Psi_{i,j}-\bar{\Psi}_{i,j})\Phi_{i,j'}^\top] \mathbb{E}_{i,j}^+[\Phi_{i,j'}(\Psi_{i,j}-\bar{\Psi}_{i,j})^\top] H_{i,j}^\top\precsim H_{i,j}H_{i,j}^\top;
\]
then with Lemma \ref{lemma:trace_inequality_2} we have (ii).
\paragraph{High probability event $\mathcal{E}_n$.}
We finally show that $\mathbb{P}(\mathcal{E}^c_n)=O(Ln^{\alpha_1+\alpha_2-1})$. When $\mathcal{E}_n$ does not happen, there exists a $j\in\{0,1,\dots, L\}$ and an
eigenvector $x$ of $\Bar B_{n,j,j}^\top \Bar B_{n,j,j}$, with $\|x\|_2=1$ such that $\Bar B_{n,j,j}^\top \Bar B_{n,j,j} x = \lambda x$ and $\lambda<c_1^2n^{-\alpha_1}/4$. In that case:
\begin{equation}
\label{eq:ct_ex}
x^\top \Bar B_{n,j,j}^\top \Bar B_{n,j,j} x < \frac{c_1^2 n^{-\alpha_1}}{4}.
\end{equation}
Since $\delta_{n,j,j}=\Bar B_{n,j,j}-B_{n,j,j}^0$, we have that:
\begin{align}
\frac{c^2n^{-\alpha_1}}{4} > x^\top \Bar B_{n,j,j}^\top \Bar B_{n,j,j} x \stackrel{(i)}{\geq} \frac{1}{2} x^\top (B_{n,j,j}^0)^\top B_{n,j,j}^0 x - 2\,\|\delta_{n,j,j}\|_{Frob}^2 \stackrel{(ii)}{\geq} \frac{c_1^2n^{-\alpha_1}}{2} -2\,\|\delta_{n,j,j}\|_{Frob}^2
\end{align}
where (i) is by Lemma \ref{lem:lower-eigenvalue}, and (ii) is because the minimum eigenvalue value of $(B_{n,j,j}^0)^\top B_{n,j,j}^0$ is at least $c_1^2n^{-\alpha}$ by Property \ref{property:weight_regularizing}(a) and the vector $x$ is unit-norm. Rearranging yields:
\begin{align}
\|\delta_{n,j,j}\|_{Frob}^2 \geq \frac{c_1^2n^{-\alpha_1}}{8}
\end{align}
Thus we have,
by Markov's inequality:
\begin{align*}
\mathbb{P}(\mathcal{E}^c_n)&\leq \mathbb{P}\left( \exists j, \|\delta_{n,j,j}\|_{Frob}^2 \geq \frac{c_1^2n^{-\alpha_1}}{8}\right)\leq \sum_{j=0}^L\mathbb{P}\left( \|\delta_{n,j,j}\|_{Frob}^2 \geq \frac{c_1^2n^{-\alpha_1}}{8}\right)\\
&\leq \sum_{j=0}^L\frac{8n^{\alpha_1}}{c_1^2}\mathbb{E}\left[\|\delta_{n,j,j}\|_{Frob}^2\right]=O(Ln^{\alpha_1+\alpha_2-1}).
\end{align*}
\section{Proof of Theorem \ref{thm:be_feasible_2}}
\label{appendix:be_feasible}
For notation convenience, we write
$F_{i,j}:=f_i(S_{i,1:j}, T_{i,1:j-1})$ and $\hat{F}_{i,j}:=\hat{f}_{i,j}(S_{i,1:j}, T_{i,1:j-1})$. Also, we let
$F_i:=(F_{i,1},\dots, F_{i,L})$, and $\hat{F}_i=(\hat{F}_{i,1}, \dots, \hat{F}_{i,L})$.
Define the variance estimate: $\widehat{\Xi}_n=n^{-1}\sum_{i=1}^n \mbox{diag}\left\{W_{i,j}\mathrm{Var}^+_{i,j}(\Psi_{i,j})W_{i,j}^\top\right\}$.
Let $\hat{\theta}_n$ be the AW-GMM estimator defined in \eqref{eq:gmm_estimator_linear}. With assumptions in Theorem \ref{thm:be_feasible_2}, the weights $H_{i,j}$ satisfy Property \ref{property:weight_regularizing} with constants $(\alpha_1, 0)$.
We can then invoke Theorem \ref{thm:consistency} and get:
\begin{equation}
\mathbb{E}[\|\hat{\theta}_{n,j}-\theta^*_j\|_2]=O\left(n^{\frac{(L-j+1)\alpha_1-1}{2}}\right), \forall j\in[0:L].
\end{equation}
Thus by Markov inequality, we have for any positive $\delta$:
\begin{equation*}
\mathbb{P}\left(\|\hat{\theta}_{n,j}-\theta^*_j\|_2\geq \delta\right)\leq O\left(n^{\frac{(L-j+1)\alpha_1-1}{2}}\right)\Big/\delta.
\end{equation*}
With $\theta^*$ lying in the interior of $\Theta$ and $\Theta$ being bounded, for $n$ larger than a constant, with probability at least $1-O(n^{\frac{(L+1)\alpha_1-1}{2}})$, $\hat{\theta}_n$ lies in the interior of $\Theta$. Define $\mathcal{E}_{int}$ as the event that $\hat{\theta}_n$ lies in the interior of $\Theta$. Then for $n$ larger than a constant, $\mathbb{P}(\mathcal{E}_{int}^c)=O(n^{\frac{(L+1)\alpha_1-1}{2}})$.
Now let's establish the strong Gaussian approximation when event $\mathcal{E}_{int}$ happens.
With $\hat{\theta}_n$ being in the interior of $\Theta$, we can invoke the first-order condition of the GMM solution:
\begin{align}
\tag{\ref{eq:link_mds_error_2}}
B_n^\top A\left(\frac{1}{n}\sum_{i=1}^n H_i \xi_i\right)
=& B_n^\top A B_n(\hat\theta_n-\theta^*),\quad \mbox{where}\quad B_n:=-\frac{1}{n}\sum_{i=1}^n H_i J_i.
\end{align}
Appendix \ref{appendix:uniform_mds} shows that the sum of MDS $\frac{1}{\sqrt{n}}\sum_{i=1}^n H_i\xi_{i}$ can be approximated uniformly by a Gaussian random variable $Z_\xi\sim \mathcal{N}\left(0, \widehat{\Xi}_n\right)$, such that
\begin{align}
\label{eq:high_prob_gaussian_approximation}
\sup_{C\in\mathcal{C}}\left| \mathbb{P}\left(n^{-\frac{1}{2}}\sum_{i=1}^n H_i\xi_{i}\in C\right)-\mathbb{P}\left( Z_\xi \in C\right)\right|
= O \left( L^2\cdot n^{-\min(\frac{1-\alpha_1}{12}, \frac{\alpha_3}{5}, \frac{\gamma}{5})}\right),
\end{align}
where $\mathcal{C}$ is the set of all convex subsets in $\mathbb{R}^{1+dL}$.
Therefore for any convex set $\bar C$,
\begin{equation}
\mathbb{P}\left( \left\{\sqrt{n}B_n^\top A B_n (\hat{\theta}_n -\theta)\in \bar C\right\}\cap \mathcal{E}_{int}\right)=\mathbb{P}\left(\left\{n^{-\frac{1}{2}}\sum_{i=1}^n B_n^\top A H_i\xi_{i}\in \bar C\right\}\cap \mathcal{E}_{int}\right) \label{eq:nor_eq_1}.
\end{equation}
Now for the convex set $\bar{C}$, define the induced set $\tilde{C}$ as follows:
\begin{equation}
\tilde{C} = \{x:B_n^\top A x \in \bar C\}. \label{eq:nor_eq_2}
\end{equation}
This definition yields:
\begin{align}
\mathbb{P}\left(n^{-\frac{1}{2}}\sum_{i=1}^n H_i \xi_{i}\in \tilde{C}\right) = \mathbb{P}\left(n^{-\frac{1}{2}}\sum_{i=1}^n B_n^\top A H_i\xi_{i}\in \bar C\right) \quad \mbox{and}\quad
\mathbb{P}\left(Z_\xi\in \tilde{C}\right) = \mathbb{P}\left( B_n^\top A Z_\xi\in \bar C\right). \label{eq:nor_eq_3}
\end{align}
Meanwhile, it's straightforward to see that $\tilde{C}$ is also a convex set; thus by the strong Gaussian approximation in \eqref{eq:high_prob_gaussian_approximation}, we have
\begin{align}
\left| \mathbb{P}\left(n^{-\frac{1}{2}}\sum_{i=1}^n H_i \xi_{i}\in \tilde{C}\right)-\mathbb{P}\left( Z_\xi\in \tilde{C}\right)\right|
= O \left( L^2\cdot n^{-\min(\frac{1-\alpha_1}{12}, \frac{\alpha_3}{5}, \frac{\gamma}{5})}\right) \label{eq:nor_eq_4}
\end{align}
Combining \eqref{eq:nor_eq_1}, \eqref{eq:nor_eq_3}, \eqref{eq:nor_eq_4}, we have for any convex set $\bar{C}\in\mathcal{C}$:
\begin{align}
&\left| \mathbb{P}\left( \sqrt{n}B_n^\top A B_n (\hat{\theta}_n -\theta)\in \bar C \right)-\mathbb{P}\left( B_n^\top A Z_\xi\in \bar C\right)\right| \label{eq:normalized_error_gaussian_approximation}\\
\leq &\left| \mathbb{P}\left( \left\{\sqrt{n}B_n^\top A B_n (\hat{\theta}_n -\theta)\in \bar C\right\}\cap\mathcal{E}_{int}\right)-\mathbb{P}\left(\left\{ B_n^\top A Z_\xi\in \bar C\right\} \cap\mathcal{E}_{int}\right)\right| + 2\mathbb{P}(\mathcal{E}_{int}^c) \nonumber\\
=& \left| \mathbb{P}\left( \left\{n^{-\frac{1}{2}}\sum_{i=1}^n B_n^\top A H_i\xi_i \in \bar C\right\}\cap\mathcal{E}_{int}\right)-\mathbb{P}\left(\left\{ B_n^\top AZ_\xi\in \bar C\right\} \cap\mathcal{E}_{int}\right)\right|+ 2\mathbb{P}(\mathcal{E}_{int}^c)\nonumber\\
\leq & \left| \mathbb{P}\left( n^{-\frac{1}{2}}\sum_{i=1}^n B_n^\top A H_i\xi_i \in \bar C\right)-\mathbb{P}\left( B_n^\top A Z_\xi\in \bar C\right)\right|+ 4\mathbb{P}(\mathcal{E}_{int}^c)\nonumber\\
\leq & \left| \mathbb{P}\left( n^{-\frac{1}{2}}\sum_{i=1}^n H_i\xi_i \in \tilde C\right)-\mathbb{P}\left( Z_\xi\in \tilde C\right)\right|+ 4\mathbb{P}(\mathcal{E}_{int}^c)
\tag{by definition of $\tilde C$ in \eqref{eq:nor_eq_2}}
\nonumber\\
&= O \left( L^2\cdot n^{-\min(\frac{1-\alpha_1}{12}, \frac{\alpha_3}{5}, \frac{\gamma}{5})}\right)+ O\left(n^{\frac{(L+1)\alpha_1-1}{2}}\right) = O \left( L^2\cdot n^{-\min(\frac{1-\alpha_1}{12}, \frac{1-(L+1)\alpha_1}{2}, \frac{\alpha_3}{5}, \frac{\gamma}{5})}\right) .\nonumber
\end{align}
We finally show that the normalizing matrix $B_n^\top AB_n$
is invertible with high probability.
Define event $\mathcal{E}^M_n$ for $B_n=-n^{-1}\sum_{i=1}^n H_iJ_i$:
\begin{align*}
\mathcal{E}^M_n = \left\{B_{n,j,j}^\top B_{n,j,j} \succeq \frac{c_1^2}{4M^2} n^{-\alpha_1} I, \ \forall l\in[0: L]\right\}.
\end{align*}
Appendix \ref{appendix:regularity_B_normal} shows that $\mathcal{E}^M_n$ happens with probability at least $1-O(Ln^{\alpha_1-1})$. When $\mathcal{E}^M_n$ happens, each $B_{n,j,j}$ has full column rank. Since $B_n$ is a block upper triangular matrix, $B_n$ also has full column rank. As a result, $B_n^\top A B_n$, which has the same rank as $B_n$ (see Lemma \ref{lemma:col_rank_matrix}), is thus invertible.
Now for any convex set $C\in\mathcal{C}$, define the induced set
\begin{equation}
\label{eq:induce_c_bar}
\bar C:=\{x: (B_n^\top A B_n)^\dagger x\in C\},
\end{equation}
which is also a convex set by construction. We have:
\begin{align*}
&\left| \mathbb{P}\left( \sqrt{n}(\hat{\theta}_n -\theta)\in C \right)-\mathbb{P}\left((B_n^\top A B_n)^\dagger B_n^\top A Z_\xi\in C\right)\right| \\
\leq &\left| \mathbb{P}\left( \left\{\sqrt{n}(\hat{\theta}_n -\theta)\in C\right\}\cap \mathcal{E}^M_n\right)-\mathbb{P}\left(\left\{(B_n^\top A B_n)^\dagger B_n^\top A Z_\xi\in C\right\} \cap \mathcal{E}^M_n\right)\right| + 2\mathbb{P}((\mathcal{E}^M_n)^c) \\
= &\left| \mathbb{P}\left( \left\{\sqrt{n}(B_n^\top A B_n)^{-1}B_n^\top A B_n(\hat{\theta}_n -\theta)\in C\right\}\cap \mathcal{E}^M_n\right)-\mathbb{P}\left(\left\{(B_n^\top A B_n)^\dagger B_n^\top A Z_\xi\in C\right\} \cap \mathcal{E}^M_n\right)\right|\\
&\quad \quad + 2\mathbb{P}((\mathcal{E}^M_n)^c)
\tag{When $\mathcal{E}^M_n$ happens, $B_n^\top A B_n$ is invertible. }
\\
= &\left| \mathbb{P}\left( \left\{\sqrt{n}(B_n^\top A B_n)^{\dagger}B_n^\top A B_n(\hat{\theta}_n -\theta)\in C\right\}\cap \mathcal{E}^M_n\right)-\mathbb{P}\left(\left\{(B_n^\top A B_n)^{\dagger}B_n^\top A Z_\xi\in C\right\} \cap \mathcal{E}^M_n\right)\right| + 2\mathbb{P}((\mathcal{E}^M_n)^c)
\tag{When $\mathcal{E}^M_n$ happens, $B_n^\top A B_n$ is invertible and $(B_n^\top A B_n)^{-1}=(B_n^\top A B_n)^\dagger$. }\\
= & \left| \mathbb{P}\left( \left\{\sqrt{n}B_n^\top A B_n(\hat{\theta}_n -\theta)\in\bar C\right\}\cap \mathcal{E}^M_n\right)-\mathbb{P}\left(\left\{ B_n^\top A Z_\xi\in \bar C\right\} \cap \mathcal{E}^M_n\right)\right| + 2\mathbb{P}((\mathcal{E}^M_n)^c)
\tag{by definition of $\bar C$ in \eqref{eq:induce_c_bar}}
\\
= & \left| \mathbb{P}\left( \sqrt{n}B_n^\top A B_n(\hat{\theta}_n -\theta)\in\bar C\right)-\mathbb{P}\left(B_n^\top A Z_\xi\in \bar C\right)\right| + 4\mathbb{P}((\mathcal{E}^M_n)^c) \\
= & O \left( L^2\cdot n^{-\min(\frac{1-\alpha_1}{12}, \frac{1-(L+1)\alpha_1}{2}, \frac{\alpha_3}{5}, \frac{\gamma}{5})}\right) + O(Ln^{\alpha_1-1})\tag{by \eqref{eq:normalized_error_gaussian_approximation}}\\
=& O \left( L^2\cdot n^{-\min(\frac{1-\alpha_1}{12}, \frac{1-(L+1)\alpha_1}{2}, \frac{\alpha_3}{5}, \frac{\gamma}{5})}\right).
\end{align*}
This concludes our proof for Theorem \ref{thm:be_feasible_2}.
\subsection{Strong Gaussian Approximation of $n^{-\frac{1}{2}}\sum_{i=1}^n H_i\xi_i$}
\label{appendix:uniform_mds}
We show the strong Gaussian approximation results of $n^{-\frac{1}{2}}\sum_{i=1}^n H_i \xi_i$. To do it, we follow three steps:
\begin{enumerate}
\item Connect $n^{-\frac{1}{2}}\sum_{i=1}^n H_i\xi_i$ with $\mathcal{N}(0, \Sigma_n)$, where recall that
\[
\Sigma_n = \mbox{diag}\left\{\Sigma_{n,0}, \Sigma_{n,1}, \dots, \Sigma_{n,L}\right\}, \quad \mbox{where} \quad \Sigma_{n,j}=n^{-1}\sum_{i=1}^n\mathbb{E}\left[\frac{F_{i,j}}{\hat{F}_{i,j}}
W_{i,j}\mathrm{Var}^+_{i,j}(\Psi_{i,j})W_{i,j}^\top
\right].
\]
\item Connect $\mathcal{N}(0, \Sigma_n)$ with $\mathcal{N}(0, \Xi_n)$, where recall that
\[
\Xi_n=\mbox{diag}\left\{
\Xi_{n,0}, \Xi_{n,1},\dots \Xi_{n,L}
\right\},\quad \mbox{where}\quad \Xi_{n,j}= n^{-1}\mathbb{E}\left[\sum_{i=1}^n W_{i,j}\mathrm{Var}^+_{i,j}(\Psi_{i,j})W_{i,j}^\top\right],
\]
and by Property \ref{property:weight_regularizing}(b), we have
\[
\Xi_n \succeq C_2 \cdot I.
\]
\item Connect $\mathcal{N}(0, \Xi_n)$ with $\mathcal{N}(0, \widehat{\Xi}_n)$.\footnote{For any almost surely positive semi-definite random matrix $\Sigma$, we use $\mathcal{N}(0,\Sigma)$
to define the distribution of $X:=\Sigma^{1/2}Z$ with $Z\sim \mathcal{N}(0,1)$ independent of $\Sigma$.} Here we define $\widehat{\Xi}_n$ as
\[
\widehat{\Xi}_n:=
\mbox{diag}\left\{
\widehat\Xi_{n,0}, \ \widehat\Xi_{n,1},\dots \widehat\Xi_{n,L}
\right\},\quad \mbox{where}\quad \widehat\Xi_{n,j}=
n^{-1}\mbox{diag}\left\{\sum_{i=1}^n W_{i,j}\mathrm{Var}^+_{i,j}(\Psi_{i,j})W_{i,j}^\top\right\}.
\]
\end{enumerate}
Collectively, we shall be able to connect $n^{-1}\sum_{i=1}^n H_i\xi_i$ with $\mathcal{N}(0, \widehat{\Xi}_n)$, which achieves the inferential result we target.
\medskip
\subsubsection{Step 1: Connect $n^{-\frac{1}{2}}\sum_{i=1}^n H_i\xi_i$ with $\mathcal{N}(0, \Sigma_n)$.}
Recall that $\{ H_i\xi_i\}_{i=1}^n$ is a martingale difference sequence. To characterize its asymptotic behavior, we shall leverage the following proposition.
\begin{restatable}[Strong Gaussian Approximation for Martingale Vectors, \cite{cattaneo2022yurinskii}]{proposition}{cattaneo}
\label{prop:uniform_clt}
Let $\{\varphi_t\}_{t=1}^T $ be $\mathbb{R}^d$-valued squared integrable martingale difference sequence adapted to $\{\mathcal{G}_t\}_{t=1}^T$. Define $v_t:=\mathrm{Var}_{t}(\varphi_t) = \mathbb{E}[\varphi_t \varphi_t^\top\mid \mathcal{G}_{t-1}]$. Define $S_T:=\sum_{t=1}^T \varphi_t$, and let $\Sigma_T:=\frac{1}{T}\sum_{t=1}^T\mathbb{E}[v_t]$ and $\Omega_T:=\sum_{t=1}^T (v_t-\mathbb{E}[v_t])$. Then there exists a $W_T \sim \mathcal{N}(0, \Sigma_T)$ such that
\begin{equation*}
\sup_{C\in\mathcal{C}}|\mathbb{P}(T^{-1/2}S_T\in C)-\mathbb{P}(W_T\in C)| \lesssim \inf_{\eta>0} \left\{
\left(\frac{\beta_2d}{\eta^3}\right)^{1/3} + \left(\frac{\sqrt{d\mathbb{E}[\|\Omega_T\|_2]}}{\eta}\right)^{2/3} +\Delta_2(\mathcal{C}, \eta)
\right\},
\end{equation*}
where:
\begin{itemize}
\item $\mathcal{C}$ is a class of measurable subsets of $\mathbb{R}^d$.
\item $\Delta_p(\mathcal{C}, \eta)$ defines the Gaussian perimetric (anti-concentration) quantity \[
\Delta_2(\mathcal{C}, \eta) = \sup_{C\in\mathcal{C}}\left\{\mathbb{P}(W_T\in C_2^\eta\setminus C)\vee \mathbb{P}(W_T\in C\setminus C_2^{-\eta})\right\},
\]
with $C_2^\eta=\left\{x\in\mathbb{R}^d:\|x-C\|_2\leq \eta\right\}, C^{-\eta}_2=\mathbb{R}^d\setminus(\mathbb{R}^d\setminus C)_2^\eta$ and $\|x-C\|_2=\inf_{x'\in C}\|x-x'\|_2$.
\item $\beta_2=\sum_{t=1}^T \mathbb{E}\left[
\|\varphi_t\|^3_2 + \|v_t^{1/2}Z_t\|_2^3
\right]$, with $\{Z_t\}_{t=1}^T$ being i.i.d.~standard Gaussian variables on $\mathbf{R}^d$ independent of $\mathcal{F}_T$.
\end{itemize}
When $\Sigma_T$ is invertible, and $\mathcal{C}$ represents the set of all convex measurable subsets of $\mathbf{R}^d$, we have:
\begin{equation}
\label{eq:strong_gaussian_cattaneo}
\sup_{C\in\mathcal{C}}|\mathbb{P}(T^{-1/2}S_T\in C)-\mathbb{P}(W_T\in C)| \lesssim \inf_{\eta>0} \left\{
\left(\frac{\beta_2d}{\eta^3}\right)^{1/3} + \left(\frac{\sqrt{d\mathbb{E}[\|\Omega_T\|_2]}}{\eta}\right)^{2/3} +\eta\sqrt{\|T^{-1}\Sigma^{-1}_T\|_F}
\right\}.
\end{equation}
\end{restatable}
\begin{remark}
Proposition \ref{prop:uniform_clt} is adapted from Proposition A.1 in \cite{cattaneo2022yurinskii} by setting the vector norm parameter $p=2$.
\end{remark}
We now write out the terms in Proposition \ref{prop:uniform_clt} regarding $\{ H_i\xi_i\}_{i=1}^{n}$:
\begin{itemize}
\item $ v_i = \mathbb{E}_{i}\left[ H_i\xi_i\xi_i^\top H_i^\top\right]$, whose $j$-th diagonal block is $\mathbb{E}_{i}\left[\frac{F_{i,j}}{\hat{F}_{i,j}}
W_{i,j}\mathrm{Var}_{i,j}^+(\Psi_{i,j})W_{i,j}^\top
\right]$.\footnote{Recall that by Lemma \ref{lemma:mds} this $v_i$ is block-diagonal.}
\item $\Sigma_{n} =n^{-1}\mathbb{E}\left[ \sum_{i=1}^{n} v_i \right]$, which is a block-diagonal matrix with its $j$-th block being $ n^{-1}\sum_{i=1}^n\mathbb{E}\left[\frac{F_{i,j}}{\hat{F}_{i,j}}
W_{i,j}\mathrm{Var}^+_{i,j}(\Psi_{i,j})W_{i,j}^\top\right]$.
\item $\Omega_{n}=\sum_{i=1}^{n}(v_i-\mathbb{E}[v_i])$, which is a block-diagonal matrix with its $j$-th block being
\[
\begin{split}
&\sum_{i=1}^{n}
\left\{\mathbb{E}_{i}\left[\frac{F_{i,j}}{\hat{F}_{i,j}}
W_{i,j}\mathrm{Var}_{i,j}^+(\Psi_{i,j})W_{i,j}^\top
\right] - \mathbb{E}\left[\frac{F_{i,j}}{\hat{F}_{i,j}}
W_{i,j}\mathrm{Var}_{i,j}^+(\Psi_{i,j})W_{i,j}^\top
\right]\right\}.
\end{split}
\]
\end{itemize}
\paragraph{Regularity of $\Sigma_{n}$.} We first lower bound the eigenvalues of $\Sigma_{n}$. Note that by construction, $\Sigma_{n}$ is a symmetric block-diagonal matrix, and thus we only need to look at its $j$-th diagonal block $\Sigma_{n,j,j}$. Define $\Xi_{n} := n^{-1}\mbox{diag}_j\left\{\mathbb{E}\left[\sum_{i=1}^n W_{i,j}\mathrm{Var}^+_{i,j}(\Psi_{i,j})W_{i,j}^\top\right]\right\}$, and so $\Xi_n$ is also a block-diagonal and symmetric matrix. We have
\begin{align*}
&\left\| \Sigma_{n,j,j} -\Xi_{n,j,j}\right\|_{Frob}^2 =
\left\|\frac{1}{n}\left(\sum_{i=1}^n
\mathbb{E}\left[\frac{F_{i,j}-\hat{F}_{i,j}}{\hat{F}_{i,j}} W_{i,j}\mathrm{Var}^+_{i,j}(\Psi_{i,j})W_{i,j}^\top\right] \right)\right\|_{Frob}^2
\\
&\leq Ld
\left( \frac{1}{n}\sum_{i=1}^n
\mathbb{E}\left[\left|\frac{F_{i,j}-\hat{F}_{i,j}}{\hat{F}_{i,j}}\right| \mathrm{Tr} \left\{W_{i,j}\mathrm{Var}^+_{i,j}(\Psi_{i,j})W_{i,j}^\top\right\}\right]\right)^2
\tag{by Lemma \ref{lemma:sum_scalar_trace}}\\
& \leq Ld \left( \frac{1}{n\sigma^2}\sum_{i=1}^n \mathbb{E}\left[
\left|F_{i,j}-\hat{F}_{i,j}\right|\mathrm{Tr}\left\{ W_{i,j}\mathrm{Var}^+_{i,j}(\Psi_{i,j})W_{i,j}^\top\right\}
\right]\right)^2
\tag{$\hat F_{i,j}\geq \sigma^2$ by condition \eqref{eq:f_convergence}}\\
&\stackrel{(iii)}{\lesssim} \left( \frac{1}{n\sigma^2}\sum_{i=1}^n \mathbb{E}\left[
\left|F_{i,j}-\hat{F}_{i,j}\right|
\right]\right)^2
\tag{
$\mathrm{Tr}(W_{i,j}\mathrm{Var}^+_{i,j}(\Psi_{i,j})W_{i,j})=O(1)$ by
Property \ref{property:weight_stabilizing}(c)
}
\\
& = O(n^{-2\gamma_F})\tag{by convergence condition in \eqref{eq:f_convergence}}.
\end{align*}
Thus we have
\begin{equation}
\label{eq:diff_sigman_xin}
\left\| \Sigma_{n} -\Xi_{n}\right\|_{Frob}^2 = O(Ln^{-2\gamma_F}).
\end{equation}
Together with Property \ref{property:weight_stabilizing}(c) that says that $\Xi_{n}\succeq c_6\cdot I$ and Lemma \ref{lem:lower-eigenvalue}, we have that
\begin{align}
\label{eq:invertability_of_sigman}
& \Sigma_{n}^2\succeq \left\{\frac{c_2^2}{2}-O(Ln^{-2\gamma_F})\right\}\cdot I.
\end{align}
Therefore, for large enough $n$, $\Sigma_{n} \succeq \frac{c_2}{2}I$ is invertible. Moreover, we have:
\begin{equation}
\label{eq:bounded_sigman_inverse_frob}
\|\Sigma_{n}^{-1}\|_{Frob}\leq \sqrt{Ld} \|\Sigma_{n}^{-1}\|_{2}\leq \frac{2\sqrt{Ld}}{c_2}=O(1).
\end{equation}
\paragraph{Magnitude of $\Omega_{n}$.}
Note that $\Omega_{n}$ by definition is block-diagonal and symmetric.
\begin{align}
& \frac{1}{n}\mathbb{E}\left[\left\|\Omega_{n}\right\|_2\right] \label{eq:omega_n_bound}\\&\leq~ \max_{j=1}^L \frac{1}{n}
\mathbb{E}\left[\left\|
\sum_{i=1}^{n}
\left\{\mathbb{E}_{i}\left[\frac{F_{i,j}}{\hat{F}_{i,j}}
W_{i,j}\mathrm{Var}_{i,j}^+(\Psi_{i,j})W_{i,j}^\top
\right] - \mathbb{E}\left[\frac{F_{i,j}}{\hat{F}_{i,j}}
W_{i,j}\mathrm{Var}_{i,j}^+(\Psi_{i,j})W_{i,j}^\top
\right]\right\}
\right\|_2\right]
\right\|_2} \nonumber\\
&\leq ~\max_{j=1}^L \frac{1}{n}
\mathbb{E}\left[\left\|
\sum_{i=1}^{n}\left\{\mathbb{E}_{i}\left[\frac{F_{i,j} - \hat{F}_{i,j}}{\hat{F}_{i,j}}W_{i,j}\mathrm{Var}_{i,j}^+(\Psi_{i,j})W_{i,j}^\top \right]\right\}
\right\|_2\right]ght\|_2}\nonumber\\
&\quad + \max_{j=1}^L \frac{1}{n}\left\|\sum_{i=1}^{n}
\mathbb{E}\left[\frac{F_{i,j} - \hat{F}_{i,j}}{\hat{F}_{i,j}}W_{i,j}\mathrm{Var}^+_{i,j}(\Psi_{i,j})W_{i,j}^\top\right]
\right\|_2 \nonumber\\
&\quad + \max_{j=1}^L \frac{1}{n} \mathbb{E}
\left\|
\sum_{i=1}^{n} \left\{\mathbb{E}_{i}[W_{i,j}\mathrm{Var}^+_{i,j}(\Psi_{i,j})W_{i,j}^\top] - \mathbb{E}[W_{i,j}\mathrm{Var}^+_{i,j}(\Psi_{i,j})W_{i,j}^\top]\right\}
\right\|_2 \nonumber\\
&\lesssim ~\max_{j=1}^L \frac{1}{n}
\mathbb{E}\left[
\sum_{i=1}^{n} \mathbb{E}_{i}\left[\left|\frac{F_{i,j} - \hat{F}_{i,j}}{\hat{F}_{i,j}}\right| \mathrm{Tr}\left\{W_{i,j}\mathrm{Var}^+_{i,j}(\Psi_{i,j})W_{i,j}^\top \right\}\right]
\right]t]
}\nonumber\\
&\quad+ \max_{j=1}^L \frac{1}{n}\sum_{i=1}^{n}
\mathbb{E}\left[\left|\frac{F_{i,j} - \hat{F}_{i,j}}{\hat{F}_{i,j}}\right| \mathrm{Tr}\left\{W_{i,j}\mathrm{Var}^+_{i,j}(\Psi_{i,j})W_{i,j}^\top\right\}\right]
\tag{by Lemma \ref{lemma:sum_scalar_trace}}
\nonumber \\
&\quad + \max_{j=1}^L \frac{1}{n} \mathbb{E}
\left\|
\sum_{i=1}^{n} \left\{\mathbb{E}_{i}[W_{i,j}\mathrm{Var}^+_{i,j}(\Psi_{i,j})W_{i,j}^\top] - \mathbb{E}[W_{i,j}\mathrm{Var}^+_{i,j}(\Psi_{i,j})W_{i,j}^\top]\right\}
\right\|_2 \nonumber\\
&\lesssim ~\max_{j=1}^L \frac{2}{n}
\mathbb{E}\left[
\sum_{i=1}^{n} \mathbb{E}_{i}\left[\left|F_{i,j} - \hat{F}_{i,j}\right|\right]
\right]t]
}
\tag{$\mathrm{Tr}\left\{W_{i,j}\mathrm{Var}^+_{i,j}(\Psi_{i,j})W_{i,j}^\top\right\}=O(1)$ by Property \ref{property:weight_stabilizing}(c) and $\hat F_{i,j}\geq \sigma^2$ by construction}
\nonumber\\
&\quad + \max_{j=1}^L \frac{1}{n} \mathbb{E}
\left\|
\sum_{i=1}^{n} \left\{\mathbb{E}_{i}[W_{i,j}\mathrm{Var}^+_{i,j}(\Psi_{i,j})W_{i,j}^\top] - \mathbb{E}[W_{i,j}\mathrm{Var}^+_{i,j}(\Psi_{i,j})W_{i,j}^\top]\right\}
\right\|_2
\nonumber
\\
& \leq O(n^{-\gamma_F}) +
\max_{j=1}^L \frac{1}{n}
\sum_{i=1}^{n}
\mathbb{E}
\left\|\mathbb{E}_{i}[W_{i,j}\mathrm{Var}^+_{i,j}(\Psi_{i,j})W_{i,j}^\top] - \mathbb{E}[W_{i,j}\mathrm{Var}^+_{i,j}(\Psi_{i,j})W_{i,j}^\top]
\right\|_2
\tag{by convergence of $\hat F_{i,j}$ and triangular inequality}\\
& \leq O(n^{-\gamma_F}) +
\max_{j=1}^L \frac{1}{n}
\sum_{i=1}^{n}
\mathbb{E}
\left\|\mathbb{E}_{i,j}[W_{i,j}\mathrm{Var}^+_{i,j}(\Psi_{i,j})W_{i,j}^\top] - \mathbb{E}[W_{i,j}\mathrm{Var}^+_{i,j}(\Psi_{i,j})W_{i,j}^\top]
\right\|_2
\tag{by Jensen's inequality}\\
&=
~O(n^{-\min(\gamma_F, \alpha_3)})
\tag{by Property \ref{property:weight_stabilizing}(c)}.
\nonumber
\end{align}
\paragraph{Regularity of $\beta_2$.} We move onto discussing $\beta_2$ defined in Proposition \ref{prop:uniform_clt} and have that
\begin{align}
\beta_2 =~& \sum_{i=1}^{n} \mathbb{E}\left[\| H_i\xi_i\|^3_2 + \|v_i^{\frac{1}{2}}Z_i\|_2^3\right]\nonumber\\
\leq~& \sum_{i=1}^{n} \sqrt{\mathbb{E}\left[\| H_i\xi_i\|^4_2\right]\,\mathbb{E}\left[\| H_i\xi_i\|^2_2\right]} + \mathbb{E}\left[\left\|\mathbb{E}_{i}\left[\mbox{diag}_j\left\{\frac{F_{i,j}}{\hat{F}_{i,j}}
W_{i,j}\mathrm{Var}_{i,j}^+(\Psi_{i,j})W_{i,j}^\top\right\}
\right]^{\frac{1}{2}}Z_{i}\right\|_2^3\right]t\|_2^3}\nonumber\\
\stackrel{(i)}{\lesssim}~& \sum_{i=1}^{n} \sqrt{Ld'\,\mathbb{E}[\| H_i\xi_i\|^4_2]} + \mathbb{E}\left[\|Z_{i}\|_2^3\right]\nonumber\\
\leq~& \sum_{i=1}^{n} \sqrt{Ld'\, \mathbb{E}\left[\| H_i\xi_i\|^4_2\right]} + \sqrt{Ld'\mathbb{E}\left[\| Z_{i}\|^4_2\right]}\nonumber\\
\stackrel{(ii)}{\leq}~& \sum_{i=1}^{n} (Ld')^{\frac{1}{2}}\sqrt{\mathbb{E}
\left[ \| H_i\xi_i\|_2^4 \right]} + (Ld')^{\frac{3}{2}} \sqrt{\kappa_4}\label{eq:beta_2}
\end{align}
where $\kappa_4$ is the fourth-moment of a standard normal random variable, (i) uses that
\[
\frac{F_{i,j}}{\hat{F}_{i,j}}
W_{i,j}\mathrm{Var}_{i,j}^+(\Psi_{i,j})W_{i,j}^\top\preceq\frac{M^2}{\sigma^2}\mathrm{Tr}(W_{i,j}\mathrm{Var}_{i,j}^+(\Psi_{i,j})W_{i,j}^\top)\cdot I \precsim I,
\]
and that
\[
\mathbb{E}\left[\| H_i\xi_{i}\|^2\right]=\mathrm{Tr}(\mathbb{E}[v_i])=\mathrm{Tr}\left(\mathbb{E}\left[\mbox{diag}_j\left\{\frac{F_{i,j}}{\hat{F}_{i,j}}
W_{i,j}\mathrm{Var}_{i,j}^+(\Psi_{i,j})W_{i,j}^\top\right\}\right]\right)=O(Ld');
\]
and (ii) uses Jensen's inequality.\footnote{For any sequence $a_1,\ldots, a_K$: $(\sum_{i=1}^K a_i^2)^2 = K^2 \left(\frac{1}{K}\sum_{i=1}^K a_i^2\right)^2 \leq K\sum_{i=1}^K\left( a_i^2\right)^2$.}
Now we consider $\mathbb{E}\left[\| H_i\xi_i\|_2^4\right]$.
\begin{align}
&\mathbb{E}[\| H_i\xi_i\|_2^4]=\mathbb{E}\left[\mathbb{E}_{i}\left[ \xi_{i}^\top H_i^\top H_i \xi_{i}\xi_{i}^\top H_i^\top H_i \xi_{i}
\right] \right]right] } \nonumber\\
& \leq (L+1)\mathbb{E}\bigg[\sum_{j=0}^L\frac{ \mathbb{E}_{i}\left[R_{i,j}^4\right]\mathbb{E}_{i}\left[ \{(\Psi_{i,j}-\bar{\Psi}_{i,j})^\top W_{i,j}^\top W_{i,j}(\Psi_{i,j}-\bar{\Psi}_{i,j})\}^2 \right]}{\hat{F}_{i,j}^2}\bigg],\label{eq:fourth_moment_feasible}
\end{align}
where the inequality is by Cauchy-Schwartz inequality.
Note that $\hat{F}_{i,j}\geq \sigma^2$, and both $R_{i,j}^4$ and $\|\Psi_{i,j}-\bar{\Psi}_{i,j}\|_2^2$ are bounded by some finite constant, thus we have $(\Psi_{i,j}-\bar{\Psi}_{i,j})(\Psi_{i,j}-\bar{\Psi}_{i,j})^\top \precsim I_d$. Hence,
\begin{align}
\eqref{eq:fourth_moment_feasible}
\lesssim~& L\cdot \mathbb{E}\left[ \sum_{j=0}^L\mathbb{E}_{i}\left[ \mathrm{Tr}\left(\left(\Psi_{i,j}-\bar{\Psi}_{i,j}\right)^\top\,W_{i,j}^\top W_{i,j}W_{i,j}^\top W_{i,j}\, \left(\Psi_{i,j}-\bar{\Psi}_{i,j}\right)\right)i\right]_{i,j}\right)}}\right]\right]}\nonumber\\
=~& L\cdot \mathbb{E}\left[ \sum_{j=0}^L \mathrm{Tr}\left(W_{i,j}^\top W_{i,j}W_{i,j}^\top W_{i,j}\mathrm{Var}_{i,j}^+(\Psi_{i,j})\right)\right]=L\cdot \mathbb{E}\left[ \sum_{j=0}^L\mathrm{Tr}\left(W_{i,j}W_{i,j}^\top W_{i,j}\mathrm{Var}^+_{i,j}(\Psi_{i,j})W_{i,j}^\top\right)\right]\nonumber\\
\stackrel{(i)}{\leq}~& L\cdot \mathbb{E}\left[\sum_{j=0}^L \mathrm{Tr}\left(W_{i,j}W_{i,j}^\top\right)\mathrm{Tr}\{W_{i,j}\mathrm{Var}^+_{i,j}(\Psi_{i,j})W_{i,j}^\top\}\right] \stackrel{(ii)}{\lesssim} Ld'\cdot\mathbb{E}\left[ \sum_{j=0}^L\mathrm{Tr}\left(W_{i,j}W_{i,j}^\top\right)\right] . \label{eq:fourth_moment_continue}
\end{align}
where (i) uses that if $A,B$ are positive semi-definite, we have $\mathrm{Tr}(AB)\leq\mathrm{Tr}(A)\mathrm{Tr}(B)$, and (ii) uses Property \ref{property:weight_stabilizing}(c).
Combining \eqref{eq:beta_2} and \eqref{eq:fourth_moment_continue}, we have that
\begin{align}
\beta_2 &\lesssim ~ \sum_{i=1}^{n} d'L\sqrt{ \mathbb{E}\left[\sum_{j=0}^L \mathrm{Tr}\left(W_{i,j}^\top W_{i,j}\right)\right]} + (Ld')^{\frac{3}{2}} \sqrt{\kappa_4} \nonumber\\
&\leq d'Ln^{\frac{1}{2}}\sqrt{\sum_{i=1}^{n}\mathbb{E}\left[\sum_{j=0}^L \mathrm{Tr}\left(W_{i,j}^\top W_{i,j}\right)\right]} + (Ld')^{\frac{3}{2}}n\sqrt{\kappa_4}\nonumber\\
& = d'L n^{\frac{1}{2}}\sqrt{\sum_{j=0}^L\sum_{i=1}^{n}\mathbb{E}\left[ \mathrm{Tr}\left(W_{i,j}^\top W_{i,j}\right)\right]} + (Ld')^{\frac{3}{2}}n\sqrt{\kappa_4}\nonumber\\
&= O\left((Ld')^{\frac{3}{2}}n^{\frac{\alpha_1}{2}+1}\right),
\label{eq:beta_2_bound}
\end{align}
where the last equality is by Property \ref{property:weight_stabilizing}(b).
\paragraph{Applying Proposition \ref{prop:uniform_clt}.}
Combining \eqref{eq:bounded_sigman_inverse_frob}, \eqref{eq:omega_n_bound}, \eqref{eq:beta_2_bound},
We have
\begin{align}
&\sup_{C\in\mathcal{C}}\left|\mathbb{P}\left(n^{-1/2}\sum_{i=1}^{n} H_i\xi_i\in C\right)-\mathbb{P}\left(\mathcal{N}(0, \Sigma_{n})\in C\right)\right|\nonumber\\
\lesssim~&\inf_{\eta>0} \left\{
\left(\frac{\beta_2dL}{\eta^3}\right)^{\frac{1}{3}} + \left(\frac{\sqrt{dL\mathbb{E}[\|\Omega_{n}\|_2]}}{\eta}\right)^{\frac{2}{3}} +\eta\sqrt{\|n^{-1}\Sigma_{n}^{-1}\|_F}
\right\} \nonumber\\
\lesssim~& \inf_{\eta>0} \left\{
\frac{L^{\frac{2}{3}}n^{\frac{1}{3}+\frac{\alpha_1}{6}}}{\eta}
+ \frac{L^{\frac{1}{3}}n^{\frac{1-\min(\gamma_F, \alpha_3)}{3}} }{\eta^{\frac{2}{3}}}
+\eta n^{-\frac{1}{2}}L^{\frac{1}{4}}
\right\}. \label{eq:gaussian_error}
\end{align}
Choose $\eta= n^{\frac{5-2\gamma'}{10}}$ where $\gamma'\in(0,\min(\gamma_F, \alpha_3)]$. Continuing the above, we have
\begin{align*}
\eqref{eq:gaussian_error}
\leq ~& L^{\frac{2}{3}}\,\min_{\gamma'\in(0, \min(\gamma_F, \alpha_3)]}\max\left( n^{\frac{\alpha_1-1}{6} + \frac{\gamma'}{5}}, n^{-\frac{\gamma'}{5}} \right)\\
=~ &
L^{\frac{2}{3}}\,\min_{\gamma'\in(0, \min(\gamma_F, \alpha_3)]} n^{-\min\left(\frac{1-\alpha_1}{6} - \frac{\gamma'}{5},\frac{\gamma'}{5}\right)}\\
=~& L^{\frac{2}{3}}\, n^{-\min\left(\frac{1-\alpha_1}{12}, \frac{\alpha_3}{5}, \frac{\gamma_F}{5}\right)}.\tag{choose $\gamma'=\min\left(\frac{5(1-\alpha_1)}{12}, \gamma_F, \alpha_3\right)$}\nonumber
\end{align*}
Collectively, we have
\begin{equation}
\sup_{C\in\mathcal{C}}\left|\mathbb{P}\left(n^{-1/2}\sum_{i=1}^n H_i\xi_i\in C\right)-\mathbb{P}\left(\mathcal{N}(0, \Sigma_{n})\in C\right)\right| = O\left(
L^{\frac{2}{3}} n^{-\min(\frac{1-\alpha_1}{12}, \frac{\alpha_3}{5}, \frac{\gamma_F}{5})}
\right)
\label{eq:step_1}
\end{equation}
\medskip
For the remaining steps $2-3$, we need to invoke the following lemma.
\begin{lemma}[see \cite{barsov1987estimates,devroye2018total}]\label{lem:tvdist} Let $\Sigma_1, \Sigma_2$ be two positive definite covariance matrices and let $\lambda_1, \ldots, \lambda_d$ be the eigenvalues of $\Sigma_1^{-1}\Sigma_2-I$. Then $\text{TV}(\mathcal{N}(0,\Sigma_1), \mathcal{N}(0, \Sigma_2))\leq 2\min(\sqrt{\sum_i \lambda_i^2},1)$.
\end{lemma}
\subsubsection{Step 2: Connect $\mathcal{N}(0, \Sigma_n)$ with $\mathcal{N}(0, \Xi_n)$.} Note that by \eqref{eq:diff_sigman_xin}, we have that $\left\|\Sigma_n - \Xi_n\right\|_{Frob}^2=O(Ln^{-2\gamma_F})$. Moreover, under Property \ref{property:weight_stabilizing}(c), $\Xi_n \succeq c_2 I $ and thus is invertible; also $\Sigma_n\succeq \frac{c_2}{2}I$ for large $n$ by \eqref{eq:invertability_of_sigman}.
Consider
\begin{align*}
\Xi_n^{-1}\Sigma_n-I = \Xi_n^{-1} \cdot \left(\Sigma_n - \Xi_n\right).
\end{align*}
Let $\{\lambda_i\}$ be the eigenvalues of $\Xi_n^{-1}\Sigma_n-I $, we have
\begin{align*}
\sum_{i}\lambda_i^2 &\lesssim\|\Xi_n^{-1}\Sigma_n-I\|_{Frob}^2 \tag{by Lemma \ref{lemma:eigenval_frob}}\\
&= \mathrm{Tr}\left(\Xi_n^{-1} \cdot \left(\Sigma_n - \Xi_n\right) \cdot \left(\Sigma_n - \Xi_n\right) \Xi_n^{-1}\right)ght) \Xi_n^{-1}}\\
& = \mathrm{Tr}\left(\Xi_n^{-2} \cdot \left(\Sigma_n - \Xi_n\right)^2\right)ight)^2}\\
& \leq \mathrm{Tr}\left(\Xi_n^{-2}\right)\mathrm{Tr}\left(\left(\Sigma_n - \Xi_n\right)^2 \right)ght)^2 }\tag{by Lemma \ref{lemma:trace_inequality}} \\
&\leq \mathrm{Tr}\left(\Xi_n^{-1}\right)^2 \left\| \left(\Sigma_n - \Xi_n\right)\right\|_{Frob}^2 \tag{since $\Xi_n^{-1}\succ 0$}\\
&= O(L^3n^{-2\gamma_F}).
\end{align*}
Applying Lemma \ref{lem:tvdist}, we thus have
\begin{equation}
\text{TV}(\mathcal{N}(0,\Sigma_n), \mathcal{N}(0, \Xi_n))=O(L^{\frac{3}{2}}n^{-\gamma_F}).
\label{eq:step_2}
\end{equation}
\medskip
\subsubsection{Step 3: Connect $\mathcal{N}(0,\Xi_n)$ with $\mathcal{N}(0, \widehat{\Xi}_n)$.}
We first show that event $\widehat{\mathcal{E}}_n=\left\{ \widehat{\Xi}_n\succeq \frac{c_6}{2} I\right\}$ happens with high probability.
Define $\overline{\Xi}_n$ as a block-diagonal matrix with its $j$-th diagonal block as,
\[
\frac{1}{n}\sum_{i=1}^{n} \mathbb{E}_{i}[W_{i,j}\mathrm{Var}^+_{i,j}(\Psi_{i,j})W_{i,j}^\top].
\]
Note that $\widehat{\Xi}_n-\overline{\Xi}_n$ is the sum of martingale difference sequence. We have
\begin{align*}
&\mathbb{E}\left[
\left\|\widehat{\Xi}_n-\overline{\Xi}_n\right\|_{Frob}^2
\right]= \mathbb{E}\left[
\mathrm{Tr}\left(\left(\widehat{\Xi}_n-\overline{\Xi}_n\right)\cdot \left(\widehat{\Xi}_n-\overline{\Xi}_n\right)\right)e{\Xi}_n\right)}
\right] = \mathbb{E}\left[
\mathrm{Tr}\left(\left(\widehat{\Xi}_n-\overline{\Xi}_n\right)^2\right)ight)^2}
\right]\\
=&\frac{1}{n^2}\sum_{i=1}^{n}\mathbb{E} \Big[ \mathrm{Tr}\big(\left\{W_{i,j}\mathrm{Var}^+_{i,j}(\Psi_{i,j})W_{i,j}^\top - \mathbb{E}_{i}[W_{i,j}\mathrm{Var}^+_{i,j}(\Psi_{i,j})W_{i,j}^\top]\right\}^2\big) \Big]\\
\leq & \frac{1}{n^2}\sum_{i=1}^{n} \mathbb{E}\left[\mathrm{Tr}\left(
\left\{W_{i,j}\mathrm{Var}^+_{i,j}(\Psi_{i,j})W_{i,j}^\top \right\}^2
\right)\right]\tag{by Lemma \ref{lemma:trace_inequality_4}} \\
=& \frac{1}{n^2}\sum_{i=1}^{n} \mathbb{E}\left[\mathrm{Tr}\left(
\left\{W_{i,j}\mathrm{Var}^+_{i,j}(\Psi_{i,j})W_{i,j}^\top \right\}
\right)^2\right] \\
=& O(Ln^{-1}).\tag{Property \ref{property:weight_stabilizing}(c)}
\end{align*}
Thus
\[
\mathbb{E}\left[
\left\|\widehat{\Xi}_n-\overline{\Xi}_n\right\|_{2}
\right]\leq \mathbb{E}\left[
\left\|\widehat{\Xi}_n-\overline{\Xi}_n\right\|_{Frob}
\right] \leq \mathbb{E}\left[
\left\|\widehat{\Xi}_n-\overline{\Xi}_n\right\|_{Frob}^2
\right]^{1/2}=O\left(L^{\frac{1}{2}}n^{-\frac{1}{2}}\right).
\]
Therefore,
\[
\mathbb{E}\left[
\left\|\widehat{\Xi}_n - {\Xi}_n\right\|_{2}
\right]\leq \mathbb{E}\left[
\left\|\widehat{\Xi}_n-\overline{\Xi}_n\right\|_{2}
\right] + \mathbb{E}\left[
\left\|\overline{\Xi}_n - {\Xi}_n\right\|_{2}
\right] = O(L^{\frac{1}{2}}n^{-\min(\frac{1}{2}, \alpha_3)}),
\]
where we use the result that $ \mathbb{E}[\|\overline{\Xi}_n-\Xi_n\|_2] = O(n^{-\alpha_3})$, which holds as proven below:
\begin{align*}
&\mathbb{E}[\|\overline{\Xi}_n-\Xi_n\|_2]\\
\leq ~&\max_{j}\mathbb{E}\left[\left\|\frac{1}{n}\sum_{i=1}^{n} \mathbb{E}_{i}[W_{i,j}\mathrm{Var}^+_{i,j}(\Psi_{i,j})W_{i,j}^\top] -\frac{1}{n}\sum_{i=1}^{n} \mathbb{E}[W_{i,j}\mathrm{Var}^+_{i,j}(\Psi_{i,j})W_{i,j}^\top] \right\|_2\right]\\
\leq~& \max_{j}\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}\left[\left\| \mathbb{E}_{i}[W_{i,j}\mathrm{Var}^+_{i,j}(\Psi_{i,j})W_{i,j}^\top] - \mathbb{E}[W_{i,j}\mathrm{Var}^+_{i,j}(\Psi_{i,j})W_{i,j}^\top] \right\|_2\right]\tag{triangular inequality}\\
\leq~& \max_{j}\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}\left[\left\| \mathbb{E}_{i,j}[W_{i,j}\mathrm{Var}^+_{i,j}(\Psi_{i,j})W_{i,j}^\top] - \mathbb{E}[W_{i,j}\mathrm{Var}^+_{i,j}(\Psi_{i,j})W_{i,j}^\top] \right\|_2\right] \tag{Jensen's inequality}\\
=~& O(n^{-\alpha_3}).\tag{by Property \ref{property:weight_stabilizing}(c)}
\end{align*}
On the other hand, Property \ref{property:weight_stabilizing}(c) says that $\Xi_n\succeq c_2\cdot I$. Applying Lemma \ref{lem:high_probability_invertible}, we have that:
\[
\mathbb{P}(\widehat{\mathcal{E}}_n^c) = O(L^{\frac{1}{2}}n^{-\min(\alpha_3,\frac{1}{2})}).
\]
Under Property \ref{property:weight_stabilizing}(c), we have that $\Xi_n \succeq c_2 I $ and thus is invertible. Now conditioning on event $\widehat{\mathcal{E}}_n$ (and thus $\widehat{\Xi}_n$ is invertible), we have
\begin{align*}
\Xi_n^{-1}\widehat{\Xi}_n-I = \Xi_n^{-1} \cdot \left(\widehat{\Xi}_n - \Xi_n\right).
\end{align*}
We overload notations and let $\{\lambda_i\}$ be the eigenvalues of $\Xi_n^{-1}\widehat{\Xi}_n-I $, we have
\begin{align*}
& \sum_{i}\lambda_i^2 \lesssim \|\Xi_n^{-1}\widehat{\Xi}_n-I\|_{Frob}^2 \tag{by Lemma \ref{lemma:eigenval_frob}}\\
&= \mathrm{Tr}\left(\Xi_n^{-1} \cdot \left(\widehat{\Xi}_n - \Xi_n\right) \cdot \left(\widehat{\Xi}_n - \Xi_n\right) \Xi_n^{-1}\right)ght) \Xi_n^{-1}}\\
& \leq \mathrm{Tr}\left(\Xi_n^{-2}\right)\mathrm{Tr}\left(\left(\widehat{\Xi}_n - \Xi_n\right)^2 \right)ght)^2 } =\mathrm{Tr}\left(\Xi_n^{-1}\right)^2 \left\| \widehat{\Xi}_n - \Xi_n\right\|_{Frob}^2= O(L^3) \left\| \widehat{\Xi}_n - \Xi_n\right\|_{2}^2.
\end{align*}
Applying Lemma \ref{lem:tvdist},
we have that, conditioning on event $\widehat{\mathcal{E}}_n $ happens,
\begin{equation}
\label{eq:xin_hat_xin_tv}
\text{TV}(\mathcal{N}(0,\Xi_n), \mathcal{N}(0, \widehat{\Xi}_n))=O(L^{\frac{3}{2}})\left\| \widehat{\Xi}_n - \Xi_n\right\|_{2},
\end{equation}
where both sides are functions of the random matrix $\widehat{\Xi}_n$, and the total variation distance on the left-hand-side is defined conditional on $\widehat{\Xi}_n$, which serves as the covariance matrix of the Gaussian distribution $\mathcal{N}(0, \widehat{\Xi}_n)$.
Then for any $C\in\mathcal{C}$, we have
\begin{align*}
&\left| \mathbb{P}\left(\mathcal{N}(0,\Xi_n)\in C \right) - \mathbb{P}\left(\mathcal{N}(0,{\widehat{\Xi}}_n)\in A\right)\right|\\
\leq & \left| \mathbb{P}\left(\{\mathcal{N}(0,\Xi_n)\in C\} \cap \widehat{\mathcal{E}}_n \right) - \mathbb{P}\left(\{\mathcal{N}(0,{\widehat{\Xi}}_n)\in C\} \cap \widehat{\mathcal{E}}_n\right)\right| + 2 \mathbb{P}(\widehat{\mathcal{E}}_n^c)\\
= & \left| \mathbb{P}\left(\mathcal{N}(0,\Xi_n)\in C \mid \widehat{\mathcal{E}}_n\right) - \mathbb{P}\left(\mathcal{N}(0,{\widehat{\Xi}}_n)\in C \mid\widehat{\mathcal{E}}_n\right)\right|\mathbb{P}(\widehat{\mathcal{E}}_n) + 2 \mathbb{P}(\widehat{\mathcal{E}}_n^c)\\
= & \left|
\mathbb{E}\left[\mathbf{1}
\left\{\mathcal{N}(0,\Xi_n)\in C\right\} -
\mathbf{1}
\left\{\mathcal{N}(0,\widehat{\Xi}_n)\in C\right\}
\mid \widehat{\mathcal{E}}_n
\right]
\right|
\mathbb{P}(\widehat{\mathcal{E}}_n)
+ 2 \mathbb{P}(\widehat{\mathcal{E}}_n^c)\\
= & \left|
\mathbb{E}\left[\mathbb{E}\left[\mathbf{1}
\left\{\mathcal{N}(0,\Xi_n)\in C\right\} -
\mathbf{1}
\left\{\mathcal{N}(0,\widehat{\Xi}_n)\in C\right\}
\mid \widehat{\Xi}_n, \widehat{\mathcal{E}}_n\right]\mid \widehat{\mathcal{E}}_n
\right]\right|}_n
}
}
\mathbb{P}(\widehat{\mathcal{E}}_n)
+ 2 \mathbb{P}(\widehat{\mathcal{E}}_n^c)\\
= & \left|
\mathbb{E}\left[\mathbb{P}\left(\mathcal{N}(0,\Xi_n)\in C\mid \widehat{\Xi}_n, \widehat{\mathcal{E}}_n\right) -\mathbb{P}\left(\mathcal{N}(0,\widehat{\Xi}_n)\in C\mid \widehat{\Xi}_n, \widehat{\mathcal{E}}_n\right)
\mid
\widehat{\mathcal{E}}_n
\right]\right|
\mathbb{P}(\widehat{\mathcal{E}}_n)
+ 2 \mathbb{P}(\widehat{\mathcal{E}}_n^c)\\
\leq & \mathbb{E}\left[ \text{TV}(\mathcal{N}(0,\Xi_n), \mathcal{N}(0, \widehat{\Xi}_n))\mid \widehat{\mathcal{E}}_n\right] \mathbb{P}(\widehat{\mathcal{E}}_n) + 2 \mathbb{P}(\widehat{\mathcal{E}}_n^c)\\
= &\mathbb{E}\left[O(L^{3/2})\left\| \widehat{\Xi}_n - \Xi_n\right\|_{2}\mathbf{1}[\widehat{\mathcal{E}}_n]\right] + 2 \mathbb{P}(\widehat{\mathcal{E}}_n^c)
\tag{applying \eqref{eq:xin_hat_xin_tv}}\\
\leq & \mathbb{E}\left[O(L^{3/2})\left\|\widehat{\Xi}_n - \Xi_n\right\|_{2}\right] + 2 \mathbb{P}(\widehat{\mathcal{E}}_n^c)
= O(L^2 n^{-\min(\alpha_3,1/2)}).
\end{align*}
Thus
\begin{align}
&\sup_{C\in \mathcal{C}}\left| \mathbb{P}\left(\mathcal{N}(0,\Xi_n)\in C\right) - \mathbb{P}\left(\mathcal{N}(0,\widehat{\Xi}_n)\in C\right)\right| = O(L^2 n^{-\min(\alpha_3,1/2)}) \label{eq:step_3}
\end{align}
\medskip
\subsubsection{Put everything together.} Combining \eqref{eq:step_1}, \eqref{eq:step_2}, \eqref{eq:step_3}, we have
\begin{align*}
&\sup_{C\in \mathcal{C}}\left|\mathbb{P}\left(n^{-1/2}\sum_{i=1}^n H_i\xi_i\in C\right)-\mathbb{P}\left(\mathcal{N}(0, \widehat{\Xi}_n)\in C\right)\right|
= O \left( L^2\cdot n^{-\min(\frac{1-\alpha_1}{12}, \frac{\alpha_3}{5}, \frac{\gamma}{5})}\right).
\end{align*}
\subsection{Asymptotic regularity of $B_n$}
\label{appendix:regularity_B_normal}
Note that $B_n$ is block upper triangular, where the block sizes are corresponding to the decomposition of $\theta=( \theta_1, \dots, \theta_L)$.
Recall that by construction, $\hat{f}\in[\sigma^2, M^2]$.
For any $j\leq j'$ we have:
\begin{align*}
B_{n,j,j'}=n^{-1}\sum_{i=1}^n H_{i,j}(\Psi_{i,j}-\bar{\Psi}_{i,j})\Phi_{i,j'}^\top .
\end{align*}
Define
\begin{align*}
B^0_{n,j,j'}=n^{-1}\sum_{i=1}^n\mathbb{E}_{i,j}[ H_{i,j}(\Psi_{i,j}-\bar{\Psi}_{i,j})\Phi_{i,j'}^\top ],\mbox{ and }\delta_{n,j,j'}=B^0_{n,j,j'}-B_{n,j,j'}.
\end{align*}
For $j\in\{1,\dots,L\}$, we have:
\begin{align*}
\lambda_{\min} ((B^0_{n,j,j})^\top B^0_{n,j,j})\geq~&\frac{1}{M} \lambda_{\min}\left(
\left\{ \frac{1}{n}\sum_{i=1}^n \mathbb{E}_{i,j}[W_{i,j}\mathrm{Cov}^+_{i,j}(\Psi_{i,j},\Phi_{i,j}) ]\right\}^\top
\left\{ \frac{1}{n}\sum_{i=1}^n \mathbb{E}_{i,j}[W_{i,j}\mathrm{Cov}^+_{i,j}(\Psi_{i,j},\Phi_{i,j}) ]\right\}
\right) \\
\geq ~&\frac{c_2^2}{M}n^{-\alpha_1}.\tag{by Property \ref{property:weight_stabilizing}(a)}
\end{align*}
Next we shall show that for $j\leq j'$, we have
\begin{align*}
\mathbb{E}[\| \delta_{n,j,j'}\|_{Frob}^2]=O(n^{-1}),
\end{align*}
and that the event $\mathcal{E}^M_n$ happens with high probability:
\begin{align*}
\mathbb{P}(\left(\mathcal{E}^M_n\right)^c)=O(Ln^{\alpha_1-1}).
\end{align*}
\subsubsection{Asymptotic negligibility of $\delta_{n,j,j'}$.} With $j\leq j'$, we have
\begin{align*}
\mathbb{E}[\|\delta_{n,j,j'}\|_{Frob}^2 ]
=~& n^{-2}\mathbb{E}\Big[\mathrm{Tr}\Big(\left\{\sum_{i=1}^n H_{i,j}\left((\Psi_{i,j}-\bar{\Psi}_{i,j})\Phi_{i,j'}^\top -\mathbb{E}_{i,j}[(\Psi_{i,j}-\bar{\Psi}_{i,j})\Phi_{i,j'}^\top ]\right)\right\}\\
&\quad\quad \times\left\{\sum_{i=1}^n H_{i,j}\left((\Psi_{i,j}-\bar{\Psi}_{i,j})\Phi_{i,j'}^\top -\mathbb{E}_{i,j}[(\Psi_{i,j}-\bar{\Psi}_{i,j})\Phi_{i,j'}^\top ]\right)\right\}^\top \Big)\Big]\\
\stackrel{(i)}{=}~& n^{-2}\sum_{i=1}^n\mathrm{Tr}\Big( \mathbb{E}\Big[H_{i,j}\left((\Psi_{i,j}-\bar{\Psi}_{i,j})\Phi_{i,j'}^\top -\mathbb{E}_{i,j}[(\Psi_{i,j}-\bar{\Psi}_{i,j})\Phi_{i,j'}^\top ]\right)\\
&\quad \quad \times \left((\Psi_{i,j}-\bar{\Psi}_{i,j})\Phi_{i,j'}^\top -\mathbb{E}_{i,j}[(\Psi_{i,j}-\bar{\Psi}_{i,j})\Phi_{i,j'}^\top ]\right)^\top H_{i,j}^\top
\Big]\Big)\\
\stackrel{(ii)}{\leq}~& n^{-2}\sum_{i=1}^n\mathrm{Tr} \left(
\mathbb{E}\left[ H_{i,j} \mathbb{E}^+_{i,j}\left[(\Psi_{i,j}-\bar{\Psi}_{i,j})\Phi_{i,j'}^\top \Phi_{i,j'}(\Psi_{i,j}-\bar{\Psi}_{i,j})^\top \right] H_{i,j}^\top \right]\right)^\top }
}\\
\stackrel{(iii)}{\lesssim}~& n^{-2}\sum_{i=1}^n\mathrm{Tr} \left(
\mathbb{E}\left[ W_{i,j} \mathbb{E}^+_{i,j}\left[(\Psi_{i,j}-\bar{\Psi}_{i,j})(\Psi_{i,j}-\bar{\Psi}_{i,j})^\top \right] W_{i,j}^\top \right]\right)^\top }
} \stackrel{(iv)}{=} O(n^{-1}).
\end{align*}
where (i) uses the fact that $\delta_{n,j,j'}$ is a sum of martingale difference sequence, (ii) uses Lemmas \ref{lemma:trace_inequality_3} \& \ref{lemma:trace_inequality_4}, (iii) uses that $\Phi_{i,j'}$ is bounded and Lemma \ref{lemma:trace_inequality_3}, and (iv) uses Property \ref{property:weight_stabilizing}(c).
\subsubsection{High probability event $\mathcal{E}^M_n$.}
We finally show that $\mathbb{P}(\left(\mathcal{E}^M_n\right)^c)=O(Ln^{\alpha_1-1})$. When $\mathcal{E}^M_n$ does not happen, by Lemma~\ref{lemma:eigenvalue_tri}, there exists a $j_0\in\{0,1,,\dots, l\}$ and $x$ with $\|x\|_2=1$ such that
\begin{equation}
\label{eq:ct_ex_normal}
x^\top B_{n,j_0,j_0}^\top B_{n,j_0,j_0} x < \frac{c_2^2n^{-\alpha_1}}{4M}.
\end{equation}
Since $\delta_{n,j_0,j_0}=B_{n,j_0,j_0}-B_{n,j_0,j_0}^0$, by Lemma~\ref{lem:lower-eigenvalue}, we have the following:
\begin{align*}
\frac{c_1^2n^{-\alpha_1}}{4M^2} > x^\top B_{n,j_0,j_0}^\top B_{n,j_0,j_0} x \geq \frac{1}{2} x^\top (B_{n,j_0,j_0}^0)^\top B_{n,j_0,j_0}^0 x - 2\|\delta_{n,j_0,j_0}\|_{Frob}^2 \geq \frac{n^{-\alpha_1}c_2^2}{2M} - 2\|\delta_{n,j_0,j_0}\|_{Frob}^2.
\end{align*}
Rearranging yields the following results:
\begin{align}
\|\delta_{n,j_0,j_0}\|_{Frob}^2 > \frac{n^{-\alpha_1}c_2^2}{8M}
\end{align}
Thus by a union bound, we have
\begin{align*}
\mathbb{P}(\left(\mathcal{E}^M_n\right)^c)
\leq~& \mathbb{P}\left(\exists j: \|\delta_{n,j,j}\|_{Frob}^2 > \frac{n^{-\alpha_1}c_2^2}{8M}\right)\\
\leq~& \sum_{j=0}^L\mathbb{P}\left( \|\delta_{n,j,j}\|_{Frob}^2 > \frac{n^{-\alpha_1}c_2^2}{8M} \right) \\
\leq~& \sum_{j=0}^L\frac{8M \mathbb{E}[\|\delta_{n,j,j}\|_{Frob}^2]}{c_2^2n^{-\alpha_1}} =O(Ln^{\alpha_1-1}) .
\end{align*}
\section{Proof of Corollaries}
\subsection{Proof of Corollary \ref{cor:consistency}}
\subsubsection{Identity weights.}
Consider $H_{i,j}$ being the identity matrix, we verify Property \ref{property:weight_regularizing}.
\emph{Property \ref{property:weight_regularizing}(a)}. We have
\begin{align*}
\frac{1}{n}\sum_{i=1}^n \mathbb{E}_{i,j}\left[H_{i,j} \mathrm{Var}^+_{i,j}(\Phi_{i,j}) \right]= \frac{1}{n}\sum_{i=1}^n \mathbb{E}_{i,j}\left[(I_{d_\tau}\otimes X_{i,j})\mathrm{Var}^+_{i,j}(T_{i,j}) (I_{d_\tau}\otimes X_{i,j})^\top\right]\\
\succeq c \frac{1}{n}\sum_{i=1}^n t_{i,j}^{-\alpha}\mathbb{E}_{i,j}\left[ (I_{d_\tau}\otimes X_{i,j}) (I_{d_\tau}\otimes X_{i,j})^\top\right] = \frac{1}{n}\sum_{i=1}^n t_{i,j}^{-\alpha}I_{d_\tau}\otimes\mathbb{E}_{i,j}\left[ X_{i,j}X_{i,j}^\top\right]\gtrsim n^{-\alpha} I.
\end{align*}
Thus:
\[
\left\{\frac{1}{n}\sum_{i=1}^n \mathbb{E}_{i,j}\left[H_{i,j} \mathrm{Var}^+_{i,j}(\Phi_{i,j}) \right]\right\}^\top \left\{\frac{1}{n}\sum_{i=1}^n \mathbb{E}_{i,j}\left[H_{i,j} \mathrm{Var}^+_{i,j}(\Phi_{i,j}) \right]\right\}\gtrsim n^{-2\alpha} I.
\]
\emph{Property \ref{property:weight_regularizing}(b)}. We have
\begin{align*}
\frac{1}{n}\sum_{i=1}^n \mathrm{Tr}\left(H_{i,j}H_{i,j}^\top\right) = \frac{1}{n}\sum_{i=1}^n \mathrm{Tr}(I)= O(1).
\end{align*}
\emph{Property \ref{property:weight_regularizing}(c)}. We have
\begin{align*}
\frac{1}{n} \sum_{i=1}^n \mathrm{Tr}\left(\mathbb{E}_{i,j}[H_{i,j} \mathrm{Var}_{i,j}^+(\Phi_{i,j}) H_{i,j}^\top]\right)= \frac{1}{n} \sum_{i=1}^n \mathrm{Tr}\left(\mathbb{E}_{i,j}[\mathrm{Var}_{i,j}^+(\Phi_{i,j})]\right)= O(1).
\end{align*}
\subsubsection{Consistency weights.} Consider $H_{i,j} = (I_{d_\tau}\otimes X_{i,j})\mathrm{Var}^+_{i,j}(T_{i,j})^{-1/2}
(I_{d_\tau}\otimes X_{i,j})^\top $.
\emph{Property \ref{property:weight_regularizing}(a)}. We have
\begin{align*}
&\frac{1}{n}\sum_{i=1}^n \mathbb{E}_{i,j}\left[H_{i,j} \mathrm{Var}^+_{i,j}(\Phi_{i,j}) \right]\\
&= \frac{1}{n}\sum_{i=1}^n \mathbb{E}_{i,j}\left[
(I_{d_\tau}\otimes X_{i,j})\mathrm{Var}^+_{i,j}(T_{i,j})^{-1/2}
(I_{d_\tau}\otimes X_{i,j})^\top
(I_{d_\tau}\otimes X_{i,j})\mathrm{Var}^+_{i,j}(T_{i,j}) (I_{d_\tau}\otimes X_{i,j})^\top
\right]\\
&= \frac{1}{n}\sum_{i=1}^n \mathbb{E}_{i,j}\left[
\|X_{i,j}\|_2^2(I_{d_\tau}\otimes X_{i,j})\mathrm{Var}^+_{i,j}(T_{i,j})^{-1/2}
\mathrm{Var}^+_{i,j}(T_{i,j})
(I_{d_\tau}\otimes X_{i,j})^\top)
\right]\\
&\succeq \frac{c^2}{n}\sum_{i=1}^n t_{i,j}^{-\alpha/2}\mathbb{E}_{i,j}\left[
(I_{d_\tau}\otimes X_{i,j})
(I_{d_\tau}\otimes X_{i,j})^\top
\right] =\frac{c^2}{n}\sum_{i=1}^n t_{i,j}^{-\alpha/2} \mathbb{E}_{i,j}\left[
I_{d_\tau}\otimes( X_{i,j} X_{i,j}^\top)
\right] \gtrsim n^{-\alpha/2} I.
\end{align*}
Thus:
\[
\left\{\frac{1}{n}\sum_{i=1}^n \mathbb{E}_{i,j}\left[H_{i,j} \mathrm{Var}^+_{i,j}(\Phi_{i,j}) \right]\right\}^\top \left\{\frac{1}{n}\sum_{i=1}^n \mathbb{E}_{i,j}\left[H_{i,j} \mathrm{Var}^+_{i,j}(\Phi_{i,j}) \right]\right\}\gtrsim n^{-\alpha} I.
\]
\emph{Property \ref{property:weight_regularizing}(b)}. We have
\begin{align*}
&\frac{1}{n}\sum_{i=1}^n \mathrm{Tr}\left(H_{i,j}H_{i,j}^\top\right) \\
&= \frac{1}{n}\sum_{i=1}^n \mathrm{Tr}( (I_{d_\tau}\otimes X_{i,j})\mathrm{Var}^+_{i,j}(T_{i,j})^{-1/2}
(I_{d_\tau}\otimes X_{i,j})^\top
(I_{d_\tau}\otimes X_{i,j})
\mathrm{Var}^+_{i,j}(T_{i,j})^{-1/2}
(I_{d_\tau}\otimes X_{i,j})^\top
)\\
&= \frac{1}{n}\sum_{i=1}^n\| X_{i,j}\|_2^{2} \mathrm{Tr}( (I_{d_\tau}\otimes X_{i,j})\mathrm{Var}^+_{i,j}(T_{i,j})^{-1/2}
\mathrm{Var}^+_{i,j}(T_{i,j})^{-1/2}
(I_{d_\tau}\otimes X_{i,j})^\top
)\\
& \leq \frac{c^{-1}}{n}\sum_{i=1}^nt_{i,j}^\alpha\| X_{i,j}\|_2^{2} \mathrm{Tr}( (I_{d_\tau}\otimes X_{i,j})
(I_{d_\tau}\otimes X_{i,j})^\top
)\\
& = \frac{c^{-1}}{n}\sum_{i=1}^nt_{i,j}^\alpha\| X_{i,j}\|_2^{2} \| X_{i,j}\|_2^{2} \mathrm{Tr}( I
) = O(n^{\alpha}).
\end{align*}
\emph{Property \ref{property:weight_regularizing}(c)}. We have
\begin{align*}
&\frac{1}{n} \sum_{i=1}^n \mathrm{Tr}\left(\mathbb{E}_{i,j}[H_{i,j} \mathrm{Var}_{i,j}^+(\Phi_{i,j}) H_{i,j}^\top]\right)\\
=& \frac{1}{n}\sum_{i=1}^n \mathrm{Tr}\Big(\mathbb{E}_{i,j}\Big[
(I_{d_\tau}\otimes X_{i,j})\mathrm{Var}^+_{i,j}(T_{i,j})^{-1/2}
(I_{d_\tau}\otimes X_{i,j})^\top
(I_{d_\tau}\otimes X_{i,j})\mathrm{Var}^+_{i,j}(T_{i,j}) \\
&\quad \quad (I_{d_\tau}\otimes X_{i,j})^\top
(I_{d_\tau}\otimes X_{i,j})
\mathrm{Var}^+_{i,j}(T_{i,j})^{-1/2}
(I_{d_\tau}\otimes X_{i,j})^\top
\Big]\Big)\\
=& \frac{1}{n}\sum_{i=1}^n \mathrm{Tr}\left(\mathbb{E}_{i,j}\left[ \|X_{i,j}\|_2^4
(I_{d_\tau}\otimes X_{i,j})\mathrm{Var}^+_{i,j}(T_{i,j})^{-1/2}
\mathrm{Var}^+_{i,j}(T_{i,j}) \mathrm{Var}^+_{i,j}(T_{i,j})^{-1/2}
(I_{d_\tau}\otimes X_{i,j})^\top
\right]\right)\\
=& \frac{1}{n}\sum_{i=1}^n\mathrm{Tr}\left( \mathbb{E}_{i,j}\left[ \| X_{i,j}\|_2^4
(I_{d_\tau}\otimes X_{i,j})
(I_{d_\tau}\otimes X_{i,j})^\top
\right]\right) =\frac{1}{n}\sum_{i=1}^n \mathrm{Tr}\left(\mathbb{E}_{i,j}\left[ \| X_{i,j}\|_2^4
I_{d_\tau}\otimes( X_{i,j} X_{i,j}^\top)
\right]\right) \\
=& \frac{1}{n}\sum_{i=1}^n K\cdot\mathrm{Tr}\left(\mathbb{E}_{i,j}\left[ \| X_{i,j}\|_2^4
X_{i,j} X_{i,j}^\top
\right]\right)= \frac{1}{n}\sum_{i=1}^n K\cdot\left(\mathbb{E}_{i,j}\left[ \| X_{i,j}\|_2^6
\right]\right) = O(K c^3)= O(1).
\end{align*}
\subsection{Proof of Corollary \ref{cor:uniform_ci}}
\label{appendix:cor_ci}
Let $a$ be a given confidence level. Define set
\[
\mathcal{A}(\ell) =\{x\in\mathbb{R}^{1+dL}: |\ell^\top x|\leq q_{\ell,\frac{a}{2}}\},
\].
By definition of $q_{\ell,\frac{\alpha}{2}}$, we have
\begin{align*}
\mathbb{P} \left(
(B_n^\top A B_n)^\dagger B_n^\top A Z_\xi\in \mathcal{A}(\ell)\right)=&
\mathbb{P}\left(
|\ell^\top (B_{n}^\top A B_{n})^{-1} B_n^\top A Z_\xi|\leq q_{\ell,\frac{\alpha}{2}}
\right)\nonumber\\
=&
\mathbb{P}\left(
|\mathcal{N}(0, n^{-1}
\ell^\top (B_{n}^\top A B_{n})^{-1} B_n^\top A
\widehat{\Xi}_n
A B_n (B_{n}^\top A B_{n})^{-1} \ell
)
|\leq q_{\ell,\frac{a}{2}}
\right)\\
=& 1-a.
\end{align*}
The Strong Gaussian Approximation results says:
\begin{equation*}
\left| \mathbb{P}\left( \sqrt{n} (\hat{\theta}_n -\theta^*)\in \mathcal{A}(\ell)\right)-\mathbb{P}\left((B_n^\top A B_n)^\dagger B_n^\top A Z_\xi\in \mathcal{A}(\ell)\right)\right|
= O \left( L^2\cdot n^{-\min(\frac{1-\alpha_1}{12},\frac{1-(L+1)\alpha_1}{2}, \frac{\alpha_3}{5}, \frac{\gamma}{5})}\right).
\end{equation*}
We thus have:
\begin{align*}
\mathbb{P}\left( \ell^\top \theta^* \in \left[
\ell^\top \hat{\theta}\pm n^{-\frac{1}{2}}q_{\ell,\frac{a}{2}}\right]\right) =& \mathbb{P}\left(
\sqrt{n}(\hat{\theta}_n-\theta^*)\in \mathcal{A}(\ell)
\right)\\
=& \mathbb{P} \left(
(B_n^\top A B_n)^\dagger B_n^\top A Z_\xi\in \mathcal{A}(\ell)\right) + O \left( L^2\cdot n^{-\min(\frac{1-\alpha_1}{12},\frac{1-(L+1)\alpha_1}{2}, \frac{\alpha_3}{5}, \frac{\gamma}{5})}\right)\\
=& (1-a) + O \left( L^2\cdot n^{-\min(\frac{1-\alpha_1}{12},\frac{1-(L+1)\alpha_1}{2}, \frac{\alpha_3}{5}, \frac{\gamma}{5})}\right).
\end{align*}
\subsection{Proof of Corollary \ref{cor:uniform_cb}}
\label{appendix:cor_cb}
Let $a$ be a given confidence level. Define set
\[
\mathcal{A}(a) = \left\{
x\in\mathbb{R}^{1+dL}: \left\|(\hat{D}_n^{\dagger})^{1/2} x\right\|_\infty\leq q_{\infty,\, 1-a}
\right\}.
\]
Then we have
\begin{align*}
\mathbb{P} \left(\theta^*\in \left[\hat{\theta}_n \pm n^{-\frac{1}{2}} q_{\infty,1-a}\cdot \sqrt{\hat{d}_n} \right]\right)
& = \mathbb{P} \left(n^{\frac{1}{2}}(\theta^*-\hat{\theta}_n)\in \left[ \pm q_{\infty,1-a}\cdot \sqrt{\hat{d}_n} \right]\right)\\
&= \mathbb{P} \left((\hat{D}_n^\dagger)^{1/2}n^{\frac{1}{2}}(\theta^*-\hat{\theta}_n)\in \left[ \pm q_{\infty,1-a}\cdot\mathbf{1}\left\{\hat d_n\neq \mathbf{0}\right\} \right]\right) \\
& = \mathbb{P} \left(\left\|(\hat{D}_n^\dagger)^{1/2}n^{\frac{1}{2}}(\theta^*-\hat{\theta}_n)\right\|_\infty\leq q_{\infty,1-a} \right)\\
&=\mathbb{P} \left( n^{\frac{1}{2}}(\theta^*-\hat{\theta}_n)\in \mathcal{A}(a) \right)
\end{align*}
On the other hand, we have:
\begin{align*}
\mathbb{P}\left(\left\{B_{n}^\top A B_{n}\right\}^\dagger B_n^\top A Z_\xi\in A(a)\right) & = \mathbb{P}\left(
\left\|(\hat{D}_n^{\dagger})^{1/2}
\left\{B_{n}^\top A B_{n}\right\}^\dagger B_n^\top A Z_\xi
\right\|_\infty\leq q_{\infty,\, 1-a}
\right)\\
& = \mathbb{P}\left(
\|\mathcal{N}(0,(\hat{D}_n^{\dagger})^{\frac{1}{2}}
\left\{B_{n}^\top A B_{n}\right\}^\dagger B_n^\top A\widehat{\Xi}_n AB_n \left\{B_{n}^\top A B_{n}\right\}^\dagger
(\hat{D}_n^{\dagger})^{\frac{1}{2}})\|_\infty
\leq q_{\infty,\, 1-a}\right)\\
&=1-a \tag{by definition of $q_{\infty,1-a}$}.
\end{align*}
Moreover, by Strong Gaussian Approximation, we have
\begin{equation*}
\left| \mathbb{P}\left(n^{\frac{1}{2}}(\theta^*-\hat{\theta}_n)\in A(a)\right)-\mathbb{P}\left(\left\{B_{n}^\top A B_{n}\right\}^\dagger B_n^\top A Z_\xi\in A(a)\right)\right| =O\left( L^2\cdot n^{-\min\left(\frac{1-\alpha_1}{12}, \frac{1-(L+1)\alpha_1}{2},\frac{\alpha_3}{5},\frac{\gamma}{5} \right)}\right)right)}}.
\end{equation*}
Collectively, we have
\begin{align*}
\mathbb{P}\left\{\theta^*\in \left[\hat{\theta}_n \pm n^{-\frac{1}{2}} q_{\infty,1-a}\cdot \sqrt{\hat{d}_n} \right]\right\}
&= \mathbb{P}\left\{\left\{n^{\frac{1}{2}}(\theta^*-\hat{\theta}_n)\in \mathcal{A}(a)\right\}\right\})\right\}}
\\
&= \mathbb{P}\left\{\left\{B_{n}^\top A B_{n}\right\}^\dagger B_n^\top A Z_\xi\in A(a)\right\}i\in A(a)} +O\left( L^2\cdot n^{-\min\left(\frac{1-\alpha_1}{12}, \frac{1-(L+1)\alpha_1}{2},\frac{\alpha_3}{5},\frac{\gamma}{5} \right)}\right)right)}}
\\
&=1-a + O\left( L^2\cdot n^{-\min\left(\frac{1-\alpha_1}{12}, \frac{1-(L+1)\alpha_1}{2}, \frac{\alpha_3}{5},\frac{\gamma}{5} \right)}\right)right)}}.
\end{align*}
\subsection{Proof of Corollary \ref{cor:normality}}
Recall that $V_{i,j}:=\mathbb{E}_{i,j}[X_{i,j}X_{i,j}^\top]$. Consider weighting scheme $H_{i,j} := \hat{f}_{i,j}^{-1/2}W_{i,j}$, where $\hat{f}_{i,j}$ satisfies Condition \eqref{eq:f_convergence} and
\[
W_{i,j}=(I_{d_{\tau}}\otimes\hat{V}_{i,j}^{-1/2}) (I_{d_\tau}\otimes X_{i,j})\mathrm{Var}^+_{i,j}(T_{i,j})^{-1/2}
(I_{d_\tau}\otimes X_{i,j})^\top\cdot\|X_{i,j}\|_2^{-2},
\]
for $ \hat{V}_{i,j}$ adapted to $\mathcal{F}_{i,j}$ and satisfying
$c^{-1}\cdot I \succeq\hat{V}_{i,j}\succeq c\cdot I$ for some $c>0$ and
\[
\frac{1}{n}
\sum_{i=1}^n \mathbb{E}\left[\left\|\hat{V}_{i,j}-V_{i,j}
\right\|_{2}\right] = O(n^{-\gamma}).
\]
We verify $W_{i,j}$ satisfies Property \ref{property:weight_stabilizing} with $\alpha_1=\alpha$ and $\alpha_3=\gamma$.
\paragraph{Checking Property \ref{property:weight_stabilizing}(a).} We have
\begin{align*}
\frac{1}{n}\sum_{i=1}^n \mathbb{E}_{i,j}\left[W_{i,j} \mathrm{Var}^+_{i,j}(\Phi_{i,j}) \right]&= \frac{1}{n}\sum_{i=1}^n \mathbb{E}_{i,j}\Big[
(I_{d_{\tau}}\otimes\hat{V}_{i,j}^{-1/2})
(I_{d_\tau}\otimes X_{i,j})\mathrm{Var}^+_{i,j}(T_{i,j})^{-1/2}
(I_{d_\tau}\otimes X_{i,j})^\top\\
&\quad\quad \quad \quad\quad \quad\quad \times \|X_{i,j}\|_2^{-2}
(I_{d_\tau}\otimes X_{i,j})\mathrm{Var}^+_{i,j}(T_{i,j}) (I_{d_\tau}\otimes X_{i,j})^\top
\Big]\\
&= \frac{1}{n}\sum_{i=1}^n (I_{d_{\tau}}\otimes\hat{V}_{i,j}^{-1/2})\mathbb{E}_{i,j}\left[
(I_{d_\tau}\otimes X_{i,j})\mathrm{Var}^+_{i,j}(T_{i,j})^{-1/2}
\mathrm{Var}^+_{i,j}(T_{i,j})
(I_{d_\tau}\otimes X_{i,j})^\top
\right]\\
&= \frac{1}{n}\sum_{i=1}^n (I_{d_{\tau}}\otimes\hat{V}_{i,j}^{-1/2})\mathbb{E}_{i,j}\left[
(I_{d_\tau}\otimes X_{i,j})
\mathrm{Var}^+_{i,j}(T_{i,j})^{1/2}
(I_{d_\tau}\otimes X_{i,j})^\top
\right]
\end{align*}
Since $\hat{V}_{i,j} \preceq c^{-1}\cdot I$, we have
both $(I_{d_{\tau}}\otimes\hat{V}_{i,j}^{-1/2})$ and $\mathbb{E}_{i,j}\left[
(I_{d_\tau}\otimes X_{i,j})
\mathrm{Var}^+_{i,j}(T_{i,j})^{1/2}
(I_{d_\tau}\otimes X_{i,j})^\top
\right]$ being symmetric positive definite. Therefore $\frac{1}{n}\sum_{i=1}^n \mathbb{E}_{i,j}\left[W_{i,j} \mathrm{Var}^+_{i,j}(\Phi_{i,j}) \right]$ has eigenvalues being real number and positive. We have:
\begin{align*}
&\lambda_{\min}\left(\left\{\frac{1}{n}\sum_{i=1}^n \mathbb{E}_{i,j}\left[W_{i,j} \mathrm{Var}^+_{i,j}(\Phi_{i,j}) \right]\right\}^\top \left\{\frac{1}{n}\sum_{i=1}^n \mathbb{E}_{i,j}\left[W_{i,j} \mathrm{Var}^+_{i,j}(\Phi_{i,j}) \right]\right\}\right)\\
\geq &\lambda_{\min}\left(\frac{1}{n}\sum_{i=1}^n \mathbb{E}_{i,j}\left[W_{i,j} \mathrm{Var}^+_{i,j}(\Phi_{i,j}) \right]\right)^2\\
\geq &
c^{1/2}\lambda_{\min}\left(
\frac{1}{n}\sum_{i=1}^n \mathbb{E}_{i,j}\left[
(I_{d_\tau}\otimes X_{i,j})
\mathrm{Var}^+_{i,j}(T_{i,j})^{1/2}
(I_{d_\tau}\otimes X_{i,j})^\top
\right]\right)^2\\
\geq & c^{1/2}\lambda_{\min}\left(
\frac{1}{n}\sum_{i=1}^n t_{i,j}^{-\alpha/2}\mathbb{E}_{i,j}\left[
(I_{d_\tau}\otimes X_{i,j})
(I_{d_\tau}\otimes X_{i,j})^\top
\right]\right)^2 \\
=& c^{1/2}\lambda_{\min}\left(
\frac{1}{n}\sum_{i=1}^n t_{i,j}^{-\alpha/2}I_{d_\tau}\otimes\mathbb{E}_{i,j}\left[
X_{i,j}X_{i,j}^\top
\right]\right)^2 = \Omega(1) n^{-\alpha}.
\end{align*}
\paragraph{Checking Property \ref{property:weight_stabilizing}(b).}
\begin{align*}
&\frac{1}{n}\sum_{i=1}^n \mathrm{Tr}\left(W_{i,j}W_{i,j}^\top\right) \\
&= \frac{1}{n}\sum_{i=1}^n \mathrm{Tr}\big( (I_{d_{\tau}}\otimes\hat{V}_{i,j}^{-1/2}) (I_{d_\tau}\otimes X_{i,j})\mathrm{Var}^+_{i,j}(T_{i,j})^{-1/2}
(I_{d_\tau}\otimes X_{i,j})^\top\\
&\quad\quad \times \| X_{i,j}\|_2^{-2}
(I_{d_\tau}\otimes X_{i,j}\cdot\| X_{i,j}\|_2^{-2})\mathrm{Var}^+_{i,j}(T_{i,j})^{-1/2}
(I_{d_\tau}\otimes X_{i,j})^\top(I_{d_{\tau}}\otimes\hat{V}_{i,j}^{-1/2})
\big)\\
&= \frac{1}{n}\sum_{i=1}^n\| X_{i,j}\|_2^{-2} \mathrm{Tr}( (I_{d_{\tau}}\otimes\hat{V}_{i,j}^{-1}) (I_{d_\tau}\otimes X_{i,j})\mathrm{Var}^+_{i,j}(T_{i,j})^{-1/2}
\mathrm{Var}^+_{i,j}(T_{i,j})^{-1/2}
(I_{d_\tau}\otimes X_{i,j})^\top
)\\
&\lesssim \frac{c^{-1}}{n}\sum_{i=1}^n\| X_{i,j}\|_2^{-2} \mathrm{Tr}( (I_{d_\tau}\otimes X_{i,j})\mathrm{Var}^+_{i,j}(T_{i,j})^{-1/2}
\mathrm{Var}^+_{i,j}(T_{i,j})^{-1/2}
(I_{d_\tau}\otimes X_{i,j})^\top
)\\
& \leq \frac{c^{-2}}{n}\sum_{i=1}^nt_{i,j}^\alpha\| X_{i,j}\|_2^{-2} \mathrm{Tr}( (I_{d_\tau}\otimes X_{i,j})
(I_{d_\tau}\otimes X_{i,j})^\top
)\\
& = \frac{c^{-2}}{n}\sum_{i=1}^nt_{i,j}^\alpha\| X_{i,j}\|_2^{-2} \mathrm{Tr}( (I_{d_\tau}\otimes X_{i,j})^\top (I_{d_\tau}\otimes X_{i,j})
)\\
& = \frac{c^{-2}}{n}\sum_{i=1}^nt_{i,j}^\alpha\| X_{i,j}\|_2^{-2} \mathrm{Tr}( \|X_{i,j}\|_2^{2} I_{d_\tau}
) = \frac{c^{-2}}{n}\sum_{i=1}^nt_{i,j}^\alpha\| X_{i,j}\|_2^{-2} \|X_{i,j}\|_2^{2} \mathrm{Tr}( I_{d_\tau}
) = O(n^{\alpha}).
\end{align*}
\paragraph{Checking Property \ref{property:weight_stabilizing}(c).}
We have
\begin{align*}
& W_{i,j} \mathrm{Var}_{i,j}^+(\Phi_{i,j}) W_{i,j}^\top\\
=& (I_{d_{\tau}}\otimes\hat{V}_{i,j}^{-1/2})
(I_{d_\tau}\otimes X_{i,j})\mathrm{Var}^+_{i,j}(T_{i,j})^{-1/2}
\left((I_{d_\tau}\otimes X_{i,j})^\top\cdot\| X_{i,j}\|_2^{-2}\right)
(I_{d_\tau}\otimes X_{i,j})\mathrm{Var}^+_{i,j}(T_{i,j}) \\
&\quad \times (I_{d_\tau}\otimes X_{i,j})^\top
(I_{d_\tau}\otimes X_{i,j}\cdot\| X_{i,j}\|_2^{-2})\mathrm{Var}^+_{i,j}(T_{i,j})^{-1/2}
(I_{d_\tau}\otimes X_{i,j})^\top
(I_{d_{\tau}}\otimes\hat{V}_{i,j}^{-1/2}) \\
=& (I_{d_{\tau}}\otimes\hat{V}_{i,j}^{-1/2})
(I_{d_\tau}\otimes X_{i,j})\mathrm{Var}^+_{i,j}(T_{i,j})^{-1/2}
(\|X_{i,j}\|_2^{2}\cdot I_{d_\tau}\cdot\|X_{i,j}\|_2^{-2})
\mathrm{Var}^+_{i,j}(T_{i,j}) \\
&\quad \times (\| X_{i,j}\|_2^{2}\cdot I_{d_\tau}\cdot\| X_{i,j}\|_2^{-2})\mathrm{Var}^+_{i,j}(T_{i,j})^{-1/2}
(I_{d_\tau}\otimes X_{i,j})^\top
(I_{d_{\tau}}\otimes\hat{V}_{i,j}^{-1/2}) \\
=& (I_{d_{\tau}}\otimes\hat{V}_{i,j}^{-1/2})
(I_{d_\tau}\otimes X_{i,j})\mathrm{Var}^+_{i,j}(T_{i,j})^{-1/2}
\mathrm{Var}^+_{i,j}(T_{i,j}) \mathrm{Var}^+_{i,j}(T_{i,j})^{-1/2} \\
&\quad
(I_{d_\tau}\otimes X_{i,j})^\top
(I_{d_{\tau}}\otimes\hat{V}_{i,j}^{-1/2}) \\
=& (I_{d_{\tau}}\otimes\hat{V}_{i,j}^{-1/2})
(I_{d_\tau}\otimes X_{i,j})
(I_{d_\tau}\otimes X_{i,j})^\top
(I_{d_{\tau}}\otimes\hat{V}_{i,j}^{-1/2}) \\
=& (I_{d_{\tau}}\otimes\hat{V}_{i,j}^{-1/2})\left(I_{d_\tau}\otimes
(X_{i,j} X_{i,j}^\top)
\right)(I_{d_{\tau}}\otimes\hat{V}_{i,j}^{-1/2}) \\
=& I_{d_{\tau}}\otimes(\hat{V}_{i,j}^{-1/2}
X_{i,j} X_{i,j}^\top
\hat{V}_{i,j}^{-1/2}) \\
\end{align*}
With $c\cdot I \preceq \hat{V}_{i,j} \preceq c^{-1}I$, $c\cdot I \preceq V_{i,j} \preceq c^{-1}I$ and $\| X_{i,j}\|_2^2\leq c^{-1}$, we have
\begin{align*}
\mathrm{Tr}(W_{i,j} \mathrm{Var}_{i,j}^+(\Phi_{i,j}) W_{i,j}^\top) =O(c^{-2}),
\end{align*}
and
\begin{align*}
\frac{1}{n}\sum_{i=1}^n \mathbb{E}_{i,j}[W_{i,j} \mathrm{Var}_{i,j}^+(\Phi_{i,j}) W_{i,j}^\top] \succeq c^2 \cdot I.
\end{align*}
Finally, we have
\begin{align*}
&\frac{1}{n}\sum_{i=1}^n \mathbb{E}\left[\|\mathbb{E}_{i,j}[W_{i,j} \mathrm{Var}_{i,j}^+(\Phi_{i,j}) W_{i,j}^\top] - I\|_2\right]
= \frac{1}{n}\sum_{i=1}^n \mathbb{E}\left[\|\mathbb{E}_{i,j}[ I_{d_{\tau}}\otimes(\hat{V}_{i,j}^{-1/2}
X_{i,j} X_{i,j}^\top
\hat{V}_{i,j}^{-1/2})] - I\|_2\right]\\
= & \frac{1}{n}\sum_{i=1}^n \mathbb{E}\left[\|I_{d_{\tau}}\otimes(\hat{V}_{i,j}^{-1/2}V_{i,j}\hat{V}_{i,j}^{-1/2}-I) \|_2\right]
= \frac{1}{n}\sum_{i=1}^n \mathbb{E}\left[\|\hat{V}_{i,j}^{-1/2}V_{i,j}\hat{V}_{i,j}^{-1/2}-I \|_2\right] \\
\leq & \frac{1}{n}\sum_{i=1}^n \mathbb{E}\left[\|\hat{V}_{i,j}^{-\frac{1}{2}}\|^2_2\|V_{i,j}- \hat{V}_{i,j} \|_2\right]
\leq \frac{c^{-1}}{n}\sum_{i=1}^n \mathbb{E}\left[\|V_{i,j}- \hat{V}_{i,j} \|_2\right] = O(n^{-\gamma}),
\end{align*}
which also leads to $\frac{1}{n}\sum_{i=1}^n \|\mathbb{E}[W_{i,j} \mathrm{Var}_{i,j}^+(\Phi_{i,j}) W_{i,j}^\top] - I\|_2= O(n^{-\gamma})$ by Jensen's inequality.
Therefore,
\begin{align*}
& \frac{1}{n} \sum_{i=1}^n \mathbb{E}\left\|
\mathbb{E}_{i,j}[W_{i,j}\mathrm{Var}^+_{i,j}(\Phi_{i,j})W_{i,j}^\top] - \mathbb{E}[W_{i,j}\mathrm{Var}_{i,j}^+(\Phi_{i,j})W_{i,j}^\top]
\right\|_{2}\\
\leq ~& \frac{1}{n} \sum_{i=1}^n \mathbb{E}\left[\|
\mathbb{E}_{i,j}[W_{i,j}\mathrm{Var}^+_{i,j}(\Phi_{i,j})W_{i,j}^\top] -I\|_2\right]+ \|\mathbb{E}[W_{i,j}\mathrm{Var}_{i,j}^+(\Phi_{i,j})W_{i,j}^\top]
-I\|_{2} = O(n^{-\gamma}).
\end{align*}