Higher-Order Debiased Estimators for General Treatment Models

bibunit[plainnat] \begin{abstract} We have witnessed tremendous progress in developing the foundation for econometrics and causal inference in the past decades. The most popular paradigm in the current literature is the classical (first-order) semiparametric theory, in which a key building block is the (first-order) influence function. However, it is now well known that estimators based on influence functions can be sub-optimal in terms of convergence rates in various settings. To address this issue, higher-order influence functions (HOIF) are developed, generalizing the classical semiparametric theory. However, most existing results in this regard focus on treatment effect parameters defined in explicit forms, such as average treatment effects (ATE). In applications, economists are often confronted with tasks of inferring more complex parameters, such as quantile treatment effects (QTE) or effects of complicated treatment regimes/policy. These more complex parameters can often only be implicitly defined as the solution to nonlinear estimating equations, which correspond to $M$/$Z$-estimation problems. Our current understanding of these problems is mainly limited to the classical semiparametric theory. Given the foundational role of HOIF for estimating explicit parameters such as ATE, a modest step toward enriching the statistical foundation of econometrics and causal inference is to develop the corresponding higher-order estimators for those more complex parameters. To this end, we consider parameters of a class of non-separable structural models in the econometrics literature and develop a class of higher-order estimators for the target parameters. Statistical properties of these higher-order estimators are derived using recent advances in $U$-processes theory. Our proposed higher-order estimators relax complexity-reducing assumptions, quantified by H\"{o}lder smoothness, imposed on the nuisance parameters compared to existing alternative estimators for many important parameters in this class, including QTE and quantile dose-response functions, among others. Numerical experiments, including simulation studies and a real data analysis, are also conducted to corroborate our theoretical claims and illustrate how higher-order estimators can be used in practice. \end{abstract} Keywords: Bias Reduction, Causal Inference, General Treatment Models, Higher-Order Influence Functions, Quantile Treatment Effects, Two-Step Procedures \allowdisplaybreaks \section{Introduction} \subsection{Motivation} Recent decades have witnessed surged interests in the econometrics and causal inference literature on the construction of nearly rate-optimal estimators for parameters/estimands/(statistical) functionals, such as the average treatment effect (ATE). One particular fruitful econometric paradigm is through the lens of semiparametric theory schick1986asymptotically, newey1990semiparametric, which delivers, among other things, the celebrated doubly-robust Augmented Inverse Probability Weighting (AIPW) estimator of the ATE robins1994estimation, hahn1998role, chernozhukov2018double, based on the (first-order) influence function of the target parameter. The asymptotic variance of such estimators is known to achieve the semiparametric efficiency bound, provided that sufficient complexity-reducing assumptions are imposed on the nuisance parameters (the outcome regression and the propensity score in the case of ATE). \subsubsection*{Sub-optimality of first-order estimators} Despite the popularity of estimators based on influence functions (henceforth, first-order estimators), growing evidence suggests that they can be sub-optimal in certain settings. Specifically, even when the parameter of interest is estimable at the parametric $n^{-1/2}$ rate, first-order estimators may fail to achieve $n^{-1/2}$-consistency, letting alone attain the semiparametric efficiency bound. This sub-optimality can arise in a variety of contexts, including, e.g., when the nuisance parameters are assumed to belong to \text{H\"{o}lder}-type smoothness classes liu2017semiparametric, robins2023minimax, structure-agnostic classes for certain types of parameters robins1997toward, balakrishnan2026fundamental, liu2024assumption, jin2025sharp, or hybrid smoothness--structure-agnostic classes bonvini2024doubly; or when the nuisance parameters are in the many-covariate or small bandwidth settings, where a direct plug-in estimator based on first-order influence functions possibly may lead to a “leave-in” bias cattaneo2018kernel, cattaneo2019two. Using ATE and classical nonparametric nuisance models as an illustration, when both the outcome regression and the propensity score belong to \text{H\"{o}lder} smoothness classes with smoothness index $s$ and variable dimension $d$, standard\footnote{When we add the qualification “standard”, we mean first-order estimators where the nuisance parameter estimates are computed from an independent nuisance sample as in chernozhukov2018double or kallus2024localized.} first-order estimators of ATE are $n^{-1/2}$-consistent when $\frac{2 s / d}{1 + 2 s / d} > 0.5$, i.e., $s / d > 0.5$. However, it is well known that $s / d > 0.25$ is both sufficient and necessary (hence minimal) for ATE to be $n^{-1/2}$-estimable robins2009semiparametric, robins2023minimax. Indeed, estimators based on \emph{higher-order influence functions} (HOIFs) robins2008higher, robins2016technical, liu2017semiparametric, a generalization of the classical first-order influence functions, have been shown to be $n^{-1/2}$-consistent for ATE when $s / d > 0.25$ robins2023minimax, liu2017semiparametric. In the remainder of this paper, we refer to these estimators as the “higher-order estimators” for convenience. More importantly, to the best of our knowledge, no other simpler estimators are $n^{-1/2}$-consistent under this minimal \text{H\"{o}lder} smoothness condition, despite years of research efforts. These results also prompt recent advances in applying the HOIF framework to statistical problems similar to but more complex than ATE, such as conditional average treatment effects kennedy2024minimax and average dose-response functions colangelo2026double, bonvini2022fast. Our approach is also connected to the literature addressing “leave-in bias” in the many-covariate and small-bandwidth settings cattaneo2018kernel, cattaneo2019two; see Section (ref), Remark (ref) and Supplementary Material Section (ref) for an extended discussion. The present paper extends this line of work to a broader class of parameters central to modern econometric analysis, which will be made precise next. \subsubsection*{Challenges due to implicitly-defined parameters} The existing HOIF framework primarily concerns parameters that are explicitly defined, including but not limited to ATE and parameters used in conditional independence tests shah2020hardness as prominent examples. However, many causal parameters of interest in modern economics applications lack closed-form expressions and are instead defined as solutions to moment conditions or optimization problems. The quantile treatment effect (QTE) provides a canonical example. In policy evaluation, when outcomes are heterogeneous or heavy-tailed, QTE offers a more informative summary of treatment effects than the mean manski2004statistical, firpo2007efficient,chernozhukov2005iv,powell2020quantile. Yet, QTE can only be characterized as the solution to a variational problem ($M$-estimation) or to a moment equation ($Z$-estimation) (see Example (ref) later). Another prominent example in the recent econometrics literature is the $\alpha$-expected shortfall (ES) fan2025policy, which measures the average outcome in the worst $\alpha$ fraction of the distribution. This quantity is of particular interest in risk management and policy evaluation under tail risk, where policymakers are concerned with extreme rather than average outcomes. Borrowing a terminology from robins2016technical, such parameters are \emph{implicitly defined}. These implicitly defined parameters can pose significant challenges to estimation and inference. Another important challenge arises from the non-separability between the nuisance parameter and the parameter of interest su2019non, chernozhukov2025linear. Using QTE with binary treatment as an example, as shown in Example (ref) and Remark (ref) later, one of the nuisance parameters, the generalized outcome regression model, depends on the true but unknown QTE parameter. To address this problem, kallus2024localized develop the \emph{localized debiased machine learning} (LDML) estimator, which first produces an initial estimator of the QTE using only the propensity score independent of the QTE itself. With this initial estimator, LDML proceeds to estimate the generalized outcome regression model and then constructs a rate doubly-robust estimator based on the influence function of the QTE. However, as shown in kallus2024localized and for the sake of completeness in Section (ref), for the LDML estimator to be $n^{-1/2}$-consistent, it demands certain strong assumptions on the rate of convergence of the initial QTE estimator, hence also on the regularity of the propensity score; see Proposition (ref). ai2021unified bypass the two-step procedure and estimate the QTE by using only the propensity score. But unlike the LDML estimator, the QTE estimator of ai2021unified fails to be rate double-robust. It is then natural to raise the following question: \begin{quote} \emph{Can estimators of QTE or other more general implicitly-defined causal parameters be constructed that improve upon existing ones based on influence functions, in terms of convergence rates and the required assumptions on the initial estimators?} \end{quote} We provide a positive answer to the above inquiry by going beyond the classical first-order semiparametric theory, in a sense formalized later in the paper. For implicitly defined parameters, however, the literature beyond first-order estimators is rather scarce, with the only exception being Section 6 of robins2016technical, which touches on the issue without providing a treatment with further technical details. Also, the regression slope in a partially linear model was considered in robins2016technical as a leading example, which has an explicit form\footnote{Partial-linear semiparametric regression model reads as $Y = \beta A + b (X) + \varepsilon$, where $\varepsilon$ has mean zero and $\varepsilon \protect\mathpalette{\protect\independenT}{\perp} A, X$. Then the regression slope $\beta \equiv {\mathbb{E}} [\mathrm{cov} (A, Y | X)] / {\mathbb{E}} [\mathrm{var} (A | X)]$ is simply the ratio between the expected conditional covariance of $A, Y$ given $X$ and the expected conditional variance of $A$ given $X$, two well-studied examples of parameters, often referred to as doubly-robust functionals rotnitzky2021characterization.}. Our paper partially fills this gap. To strike a balance between the generality and practical utility of our results, we opt for considering a broad class of parameters defined via generic moment equations with possibly non-smooth, non-separable generalized residual functions (see (ref)). This formulation accommodates discrete, continuous, and more complicated treatment types in a unified way, and it aligns with a growing econometrics literature that represents causal targets through generic moments. For example, ao2021multivalued use a generic moment function applied to counterfactual distributions under multivalued treatments. Within this unified formulation, we cover a wide range of parameters encountered in applied causal inference, including ATE, QTE, and more generally the General Treatment Models (GTMs) recently coined in the econometrics literature ai2021unified, chen2025local. Extending the HOIF framework to parameters of GTMs is thus instrumental towards building a broadly applicable point and interval estimation scheme for parameters encountered in econometrics. \subsubsection*{Summary of main results and contributions} Compared to the HOIF theory for explicitly defined parameters, the following challenges remain to be addressed, which constitute the main contributions of our paper. \begin{enumerate}[label = (\arabic*)] • We derive HOIFs for causal parameters implicitly defined through GTMs, facilitating the construction of improved higher-order estimating equations, compared to state-of-the-art estimators based on the first-order theory, such as the LDML estimator. While we mainly analyze estimators solving second-order estimating equations, we refer to them as higher-order estimators to avoid introducing additional terminology and we also present the full $m$-th order construction in Supplementary Material Section (ref). Among the parameters of GTMs, the QTE with binary treatment is one of the most commonly encountered examples in practice. We obtain a $n^{-1/2}$-consistent estimator for QTE under the weakest \text{H\"{o}lder} smoothness assumptions on the nuisance parameters to date, almost matching the minimal smoothness assumptions for ATE to be $n^{-1/2}$-estimable. • State-of-the-art first-order estimators require an initial estimator of the target parameter to estimate the nuisance parameters. This extra step incurs an additional rate-condition on the initial estimator and, consequently, on certain nuisance parameter(s). As we will see in Section (ref), for the QTE with binary treatment, our proposed higher-order estimators can drastically relax this rate-condition, highlighting another advantage of higher-order estimators (see Table (ref) in Section (ref)). • The original formulation of GTMs by ai2021unified assumes a parametric specification for the target parameters. When the treatment takes finitely many discrete values, finite-dimensional GTMs are effectively nonparametric. However, in many modern applications, treatment variables are often continuous. It may then be necessary to consider infinite-dimensional GTMs to avoid model misspecification bias; or otherwise, the finite-dimensional GTMs can at best be interpreted as a projection of the true causal or structural parameters onto a parametric model. To our knowledge, colangelo2026double and bonvini2022fast are among the first to consider estimating the average dose-response function (see Example (ref) later) using higher-order estimators when the treatment is continuous. To incorporate such scenarios, we extend our higher-order estimators to nonparametric GTMs, thereby generalizing the results of colangelo2026double and bonvini2022fast to include, for example, the quantile dose-response function (see Example (ref) later). \end{enumerate} Before proceeding, we illustrate the value of incorporating HOIF theory for GTM parameters through a simulation study; full details are provided in Sections (ref) and (ref). For clarity of exposition, we present the discussion at a high level and defer some technical rigor to later sections. Figure (ref) reports the simulation results. The parameter of interest $\beta^{*}$ is the QTE at the $25\%$ quantile, with binary treatment and a one-dimensional covariate that suffices to control for confounding. The histograms display the simulated distribution of the centered estimator, $\widehat{\beta}-\beta^{*}$. In the left panel of Figure (ref), we adopt the aforementioned state-of-the-art LDML estimator kallus2024localized, with nuisance parameters (see Example (ref) later) estimated by neural networks chen2024causal. We design the simulation such that the nuisance parameters have low smoothness in \text{H\"{o}lder} sense (with smoothness close to $s = 0.4$). As a result, even LDML with neural networks fails to deliver an estimator that is close to being unbiased, which is not too surprising given the empirical evidence on the difficulty of fitting functions of low regularity by neural networks in practice xu2022deepmed. In contrast, the right panel of Figure (ref) shows that, after centering by the true QTE, our HOIF-based estimator is well centered around zero, though with slightly higher variance than the LDML estimator. This empirical result demonstrates that the HOIF framework can be a useful tool even when modern neural networks are employed for nonparametric first-step estimation, complementing state-of-the-art methods based on classical semiparametric theory. \begin{figure}[htbp] \caption{Comparison between the first-order (LDML) estimator and the higher-order estimator (HOE) of QTE. The histograms over 1000 Monte Carlo simulations are displayed, and the red curves are the corresponding normal probability density function with the same mean and variance as the Monte Carlo mean and variance.} \end{figure} \subsection{A survey of related works} Our work is related to several strands of literature in econometrics and statistics. We now give a brief (and very likely incomplete) survey. \begin{enumerate}[label = (\arabic*)] • \textbf{Two-step procedures.} The problem studied in this paper is related to classical and modern semiparametric two-step estimation in econometrics, where some first-stage nuisance components may be indexed by the target parameter. Classical semiparametric theory typically imposes conditions under which the first-stage bias is asymptotically negligible, so that the limiting distribution of the target estimator is not affected; see, among others, newey1990semiparametric, newey1994large,chen2003estimation and the references therein. More recent work develops estimation and inference procedures that remain valid when first-stage bias is non-negligible. For example, cattaneo2018kernel and cattaneo2019two study two-step estimators where the first-stage bias can enter the asymptotic distribution of the target parameter, using kernel-based and linear-regression-based nuisance estimators, respectively. Specifically, cattaneo2019two consider a setting with many covariates, where the first-stage nuisance estimator is obtained by linear regression. They show that a non-negligible leave-in bias arises when the same sample is used for both steps, and propose jackknife-type methods to consistently estimate and remove this bias, yielding a bias-corrected estimator. cattaneo2018kernel develop a small-bandwidth asymptotic framework for kernel-based first-step estimators, where the smoothing bias is not assumed negligible. They show that the standard nonparametric bootstrap automatically corrects for this bias, so that percentile confidence intervals achieve correct coverage. When the asymptotic bias cannot be consistently estimated, cavaliere2024bootstrap propose a distributional adjustment to restore bootstrap validity. Their method first obtains a bootstrap p-value, which is not asymptotically uniform due to the bias. They show that this p-value converges to a non‑uniform limit that does not depend on the bias. By applying a probability integral transformation estimated via a second bootstrap, they obtain a corrected $p$-value that is asymptotically uniform, thus restoring valid inference without directly estimating the bias. cattaneo2025higher develop higher-order distributional approximations for kernel-based estimators of the density-weighted average derivative. By employing Edgeworth expansions, they refine the large-sample approximation of the estimator's distribution, which leads to improved finite-sample inference, particularly enhancing the coverage accuracy of confidence intervals. Another broad class of approaches focuses on estimators with a small bias property (SBP) newey2004twicing, meaning that the bias of the target estimator vanishes faster than the bias of the first-stage nuisance estimators. newey2004twicing shows that a class of semiparametric estimators based on twicing kernels satisfy this property. Tackling a similar problem and motivated by the HOIF theory, chen2024method proposed moment-based $\sqrt{n}$-consistent estimators of a class of causal parameters even when the nuisance parameters, assumed to be generalized linear models, are not consistently estimable. • \textbf{Influence functions and semiparametric inference. } Our work builds on and extends semiparametric and nonparametric inference based on influence functions, which play a central role in semiparametric efficiency theory. Influence-function-based estimation is a classical topic in semiparametric theory and has been studied for decades. For example, schick1986asymptotically develops a general method for constructing asymptotically efficient estimators in semiparametric models. newey1990semiparametric provides an introduction to semiparametric efficiency bounds, discussing their nature, calculation, and application in constructing semiparametric estimators and deriving their limiting distributions. For a more detailed and comprehensive treatment of influence functions for general semiparametric models and statistical inference, see bickel1998efficient. As mentioned, our results can be viewed as a further development of the higher-order generalization of the classical influence-function-based approaches robins2008higher. For binary treatment, robins1994estimation and hahn1998role derived the efficient influence functions for ATE and ATT (average treatment effect on the treated). cattaneo2010efficient considered a general treatment effect model with possibly over-identified and non-smooth moment conditions under multi-valued treatments, covering a broad class of treatment-effect parameters including average and quantile treatment effects, and established the corresponding efficiency bounds. Building on these efficiency results, the relevant literature focuses on developing practical estimators for treatment effects; see, for example, hirano2003efficient and the references therein. Recent work also leverages influence functions to obtain valid inference with high-dimensional or machine-learned nuisance components chernozhukov2018double. colangelo2026double further extend DML-type ideas to kernel-based inference on the average dose--response function under continuous treatments. rotnitzky2021characterization and chernozhukov2022locally characterized a large class of explicitly-defined parameters, including ATE and ATT, for which DML-type estimators can be constructed based on their influence functions. newey2018cross developed refined DML-type estimators that achieve faster convergence rates under smoothness classes with different nuisance parameters estimated from different independent subsamples. Papers mentioned above mainly concern explicitly-defined parameters, including ATE and ATT. In terms of implicitly-defined parameters, firpo2007efficient studies efficient influence functions for QTE with binary treatment. ao2021multivalued develop a unified moment-based framework for multi-valued treatments, covering means, quantiles, and distribution functions, and propose an efficient influence function based estimator for multivalued treatment effects for the treated. Their approach also enables decomposition analysis, separating wage structure effects from composition effects. Applying their method to evaluate the Workforce Investment Act (WIA) program, they find that heterogeneity in participation levels is an important dimension for evaluating social programs, demonstrating the practical relevance of their methods for policy evaluation. Complementarily, ai2021unified provide a unified framework for continuous treatments and derive efficient influence functions for general treatment models. As mentioned, more recently, kallus2024localized propose the LDML estimator for the QTE, along the line of DML estimators. \end{enumerate} \subsection{Notation} Before proceeding, we collect some frequently used notation. Let $L_{2, 0}$ denote the space of mean zero $L_{2}$ functions with respect to a dominating measure $\nu$. Without loss of generality, we take $\nu$ to be the Lebesgue measure throughout the paper. We write $a_n \lesssim b_n$ (resp., $a_n \gtrsim b_n$) if $a_n \leq C b_n$ (resp., $a_{n} \geq C b_{n}$) for some universal constant $C > 0$ independent of $n$; $a_n \asymp b_n$ if $a_n \lesssim b_n$ and $b_n \gtrsim a_n$; and $a_n \gg b_n$ (resp., $a_n \ll b_n$) if $b_n/a_n \to 0$ (resp., $b_n/a_n \to \infty$). For $a,b \in {\mathbb{R}}$, let $a \vee b = \max\{a,b\}$ and $a \wedge b = \min\{a, b\}$. For a vector $\bar{a} \in {\mathbb{R}}^{q}$, let $\| \bar{a} \| \coloneqq (\bar{a}^{\top}\bar{a})^{1/2}$ denote its Euclidean norm and $\bar{a}^{\otimes 2} = \bar{a}\,\bar{a}^{\top} $ denote the tensor (Kronecker) product. For a matrix $A\in{\mathbb{R}}^{p\times q}$ we write its operator norm as $ \Vert A\Vert_{\mathrm{op}} \coloneqq \inf\{c>0:\,\Vert Av\Vert \le c\Vert v\Vert \ \text{for every}\ v\in{\mathbb{R}}^{q}\} $. Given any function $f: {\mathcal{O}} \to {\mathbb{R}}^{m}$ with fixed integer $m \geq 1$, we denote its $L_{\infty}$-norm as $\|f (\cdot)\|_{\infty} \coloneqq \sup_{o \in {\mathcal{O}}} \|f (o)\|$. The higher-order estimators introduced later involve $L_{2} (\mathbb{P})$ projections onto the span of a size-$k$ dictionary $\bar{\phi}_{k} \coloneqq (\phi_{1}, \ldots, \phi_{k})$, where each $\phi_{j}: {\mathcal{O}} \rightarrow {\mathbb{R}}$. For any $f: {\mathcal{O}} \rightarrow {\mathbb{R}}$, this projection is denoted by $\Pi (f \mid \bar{\phi}_{k}) (\cdot) \coloneqq {\mathbb{E}} \{f (O) \bar{\phi}_{k} (O)\}^{\top} \Sigma_k^{-1} \bar{\phi}_{k} (\cdot)$, where $\Sigma_{k} \coloneqq {\mathbb{E}} \{\bar{\phi}_{k} (O)^{\otimes 2}\}$ is the population Gram matrix. The orthocomplement of this projection is $\Pi^{\perp} (f \mid \bar{\phi}_{k}) (\cdot) \coloneqq f (\cdot) - \Pi (f \mid \bar{\phi}_{k}) (\cdot)$. For any measurable function $h: {\mathcal{O}}^{m} \rightarrow {\mathbb{R}}$, let ${\mathbb{U}}_{n, m}$ be the $m$-th order $U$-statistic operator, \begin{align*} {\mathbb{U}}_{n, m} h \equiv {\mathbb{U}}_{n, m} h (O_{1}, \cdots, O_{m}) = \frac{(n - m)!}{n!} \sum_{1 \leq i_{1} \neq \cdots \neq i_{m} \leq n} h (O_{i_{1}}, \cdots, O_{i_{m}}). \end{align*} ${\mathbb{U}}_{n, 1} \equiv {\mathbb{P}}_{n}$ is understood to be the empirical mean. We next review some basic concepts from empirical processes theory; see van2023weak for further details. For a probability measure $Q$ and a constant $p > 1$, let $\| \cdot \|_{Q,p}$ denote the $L_{p} (Q)$-seminorm, i.e., $\|f\|_{Q, p} \coloneqq \{\int |f (x)|^{p} {\mathrm d} Q (x)\}^{1/p}$ for finite $p$ while $\|f\|_{Q,\infty}$ denotes the essential supremum of $|f|$ with respect to $Q$. Let ${\mathcal{F}}$ be a class of symmetric measurable functions $f : {\mathcal{O}}^{m} \rightarrow {\mathbb{R}}$ with envelope $\mathrm{F}$, that is, $\sup_{f\in{\mathcal{F}}} |f| \leq \mathrm{F} $. If $\| \mathrm{F} \|_{Q,p} > 0$, we write $\mathsf{N} ( {\mathcal{F}},\|\cdot\|_{Q,p} , \epsilon\|\mathrm{F}\|_{Q,p})$ for the minimal number of $\|\cdot\|_{Q,p}$-balls of radius $\epsilon\|\mathrm{F}\|_{Q,p}$ needed to cover ${\mathcal{F}}$. We finally recall the following definition. \begin{definition}[VC type class] A function class ${\mathcal{F}}$ with envelope $\mathrm{F}$ is said to be of VC type with characteristic $(A, v)$ if $\sup_{Q} \mathsf{N} \left( {\mathcal{F}}, \| \cdot \|_{Q, 2}, \epsilon \| \mathrm{F} \|_{Q, 2} \right) \leq (A / \epsilon)^{v}$ for all $\epsilon \in (0, 1]$, where $\sup_{Q}$ denotes the supremum over all discrete finite probability measures. \end{definition} For vector-valued functions ${\mathcal{F}} \coloneqq \{f: \mathcal O^{m} \to {\mathbb{R}}^{q}\}$ (with fixed $q$), these definitions are interpreted with $|\cdot|$ replaced by the Euclidean norm $\|\cdot\|$ on ${\mathbb{R}}^{q}$. All mentioned entropy and maximal-inequality results remain valid (up to universal constants). \subsection{Organization} The remainder of the paper is organized as follows. Section (ref) introduces the observation scheme and the general treatment model (GTM) and reviews the state-of-the-art estimators for GTM based on (first-order) influence functions from the standard semiparametric theory. Section (ref) presents the main results: In Section (ref), using the QTE as an illustration, we motivate the importance of introducing higher-order estimators to reduce bias and provide the intuition of why HOIFs can be used for bias correction; Section (ref) introduces the HOIF-based estimator for finite-dimensional GTM parameters, and Section (ref) establishes its statistical properties; Section (ref) specializes the results obtained in Section (ref) to QTE. Section (ref) addresses the estimation of infinite-dimensional GTM parameters, allowing greater modeling flexibility. Section (ref) presents empirical studies of the HOIF estimator. In particular, Section (ref) investigates the finite-sample performance of the HOIF estimator in low-regularity settings using simulation experiments, followed by a real data analysis in Section (ref). Section (ref) concludes the paper with a discussion of caveats and future directions. Proofs and less essential technical results are collected in the Supplementary Material. \section{Problem Setup and Review of Existing Results} \subsection{Causal parameters defined via general treatment models} Let $T \in {\mathcal{T}}$ be a general treatment variable that can be discrete, continuous, or a mix of both. Let $Y (t) \in {\mathcal{Y}}\subset\mathbb{R}$ denote the real-valued potential outcome when $T$ is intervened/set to $t$, and let $Y$ be the observed outcome. To control for confounding, we observe a vector of compactly supported baseline covariates $X \in {\mathcal{X}} \equiv [0, 1]^{d} \subseteq {\mathbb{R}}^{d}$, for some integer $d \geq 1$. Let ${\mathbb{P}}$ denote the joint distribution of the observed data $(X, T, Y)$ and ${\mathbb{E}}$ the corresponding expectation. Let $\Gamma: {\mathcal{Y}} \times {\mathcal{T}} \times {\mathcal{B}} \rightarrow {\mathbb{R}}^{p}$ be a known generalized residual function, which can be nonsmooth and non-separable. We are interested in the parameter $\bm{\beta}^*=(\beta^*_0,\beta^*_1,...,\beta^*_{p-1})^{\top}\in {\mathcal{B}} \subset {\mathbb{R}}^{p}$, defined as the \emph{unique} solution to the following integral moment equations: \begin{equation} \int_{\mathcal{T}} {\mathbb{E}} \{\Gamma (Y (t), t, \bm{\beta}^{\ast})\} \omega (t) {\mathrm d} t = 0, \end{equation} where $\omega(\cdot): {\mathcal{T}} \rightarrow {\mathbb{R}}$ is a user-specified weighting function. The general formulation (ref) can incorporate most causal parameters in the existing literature with flexible specifications of $\Gamma(\cdot)$ and $\omega(\cdot)$; see Examples (ref)--(ref) later in this section. We refer to (ref) as the \emph{general treatment model} (GTM), a unified moment-based framework for binary, multi-valued, and continuous treatments that covers means, quantiles, distribution functions, and various treatment effects. Similar frameworks are also studied in cattaneo2010efficient, ao2021multivalued, ai2021unified. Our goal is to conduct statistical inference on $\bm{\beta}^{\ast}$ using $n$ independent and identically distributed (i.i.d.) observations $\{O_{i} \equiv (X_{i}, T_{i}, Y_{i}), i = 1, \cdots, n\} \overset{\mathrm{i.i.d.}}{\sim} \mathbb{P}$, under weak complexity-reducing assumptions on the nuisance parameters (e.g., the smoothness/rate conditions imposed in Assumption (ref)), to be defined soon. Since (ref) involves the potential outcome $Y (t)$, we need to impose the following identification assumptions to turn (ref) into a functional of the observed data distribution. Let $\mathsf{p}_{T\mid X}(t\mid x)$ denote the conditional density/mass function of $T$ given $X=x$, referred to as the \emph{generalized propensity score}. \begin{assumption} We assume the following: (i) (Consistency) $Y \equiv Y (T)$; (ii) (Positivity) The generalized propensity score satisfies $0 < c_1 \leq \mathsf{p}_{T|X}(t|x) \leq c_2 < \infty$ for all $t \in \mathcal{T}$ with $\omega(t) > 0$, for constants $c_1, c_2$; (iii) (Unconfoundedness) $Y(t) \protect\mathpalette{\protect\independenT}{\perp} T \mid X$ for all $t \in \mathcal{T}$. \end{assumption} Assumption (ref) comprises the standard identification conditions for causal effects in observational studies. Notably, the stated positivity assumption allows practitioners to choose weight tailored for regions of ${\mathcal{T}}$ where positivity holds. Under Assumption (ref), $\bm{\beta}^{\ast}$, which is defined through the potential-outcome moment equations (ref), is identified as the unique solution to the following weighted moment equations based on the observed data distribution: \begin{equation} {\mathbb{E}} \{\xi (X, T) \Gamma (Y, T, \bm{\beta}^*) \}= 0, \end{equation} where the function $\xi (x, t) \coloneqq {\omega (t)}/{\mathsf{p}_{T | X} ( t |x)}$ is also called the \textit{stabilized weight} in the epidemiology literature hernan2000marginal. In this paper, we focus on two types of $\omega$: (1) $\omega (\cdot) \equiv \mathsf{p}_{T} (\cdot)$, the marginal distribution of the treatment $T$, which yields a finite-dimensional causal parameter (see Section (ref) and Section (ref)); (2) $\omega (\cdot) \equiv \mathsf{p}_T(\cdot)\delta_t(\cdot)$, the point mass at a particular value $t$ of $T$, which yields a nonparametric causal curve as a function of $t$ (see Section (ref)). Our framework can be readily extended to the case where $\omega (\cdot)$ is any known function of $T$ or additionally depends on $X$ (see Example (ref)), but we focus on these two choices to avoid unnecessary complications. We present several commonly encountered examples of treatment effect parameters defined via GTM (ref), demonstrating the importance of developing new and improved statistical procedures for this general class. \begin{example}[Average Treatment Effect (ATE)] When $\Gamma (Y(t),t,\bm{\beta}) = ( (Y(t) - \beta_0)\cdot (1 - t), (Y(t) - \beta_1)\cdot t )^{\top}$, $\omega(t) = \mathsf{p}_{T}(t)$ and ${\mathcal{T}} = \{0, 1\}$, where $\bm{\beta} = (\beta_{0}, \beta_{1})^{\top}$, the solution to (ref) is $\beta^*_{0} = {\mathbb{E}} [Y (0)]$ and $\beta^*_{1} = {\mathbb{E}}[Y (1)]$. ATE has been widely studied in the econometric literature; see hahn1998role,hirano2003efficient,chernozhukov2018double, among others. \end{example} \begin{example}[Quantile Treatment Effect (QTE)] When $\Gamma(Y(t),t,\bm{\beta}) = ( (\tau - \mathbbm{1} \{Y(t) \leq \beta_0\})(1 - t), (\tau - \mathbbm{1} \{Y(t) \leq \beta_1\})t)^{\top}$, $\omega(t) = \mathsf{p}_{T}(t)$ and ${\mathcal{T}} = \{0, 1\}$, where again $\bm{\beta} =(\beta_0,\beta_{1})^{\top}$, the solution to (ref) is $\beta^*_0 = \inf\{q: \mathbb{P}(Y (0) \leq q) \geq \tau\}$ and $\beta^*_{1} = \inf\{q: \mathbb{P}(Y (1) \leq q) \geq \tau\}$, which are the $\tau^{th}$ quantiles of the potential outcomes. QTE has been studied in chernozhukov2005iv,firpo2007efficient,kallus2024localized, among others. We will revisit QTE in Sections (ref) and (ref) later as the most important concrete example of GTM parameters. \end{example} \begin{example}[Distributional Treatment Effect (DTE)] For any $y\in\mathbb{R}$, let $\Gamma (Y(t),t,\bm{\beta}) = ( (\mathbbm{1}\{Y(t)\leq y\} - \beta_0 )( 1- t) , (\mathbbm{1}\{Y(t)\leq y\} - \beta_1) t)^{\top}$, $\omega(t) = \mathsf{p}_{T}(t)$, and ${\mathcal{T}} = \{0, 1\}$, where $\bm{\beta} = (\beta_{0}, \beta_{1})^{\top}$. The solution to (ref) is $\beta^*_{0} = \mathbb{P} (Y (0) \leq y)$ and $\beta^*_{1} = \mathbb{P} (Y (1) \leq y)$. DTE has been studied in chernozhukov2013inference and donald2014estimation, among others. \end{example} \begin{example}[$\alpha$-Expected Shortfall (ES)] Fix $\alpha\in(0,1)$, let ${\mathcal{T}}=\{0,1\}$ and set $\omega(t)\equiv \mathsf{p}_T(t)$. Define \begin{equation*} \Gamma\big(Y(t),t,\bm{\beta}\big) = t \cdot \begin{pmatrix} \alpha-\mathbbm{1}\{Y(t)\le \beta_0\}\\left[3pt] \dfrac{\mathbbm{1}\{Y(t) \le \beta_0\}}{\alpha}\big(Y(t)-\beta_1\big) \end{pmatrix}, \end{equation*} with $\bm{\beta} =(\beta_0,\beta_{1})^{\top}$. Then the solution to (ref) is given by \begin{align*} \beta_0^* = \inf\{q: \mathbb{P}(Y (1) \leq q) \geq \alpha\}, \qquad \beta_1^* = \mathbb{E} (Y(1)\mid Y(1) \le \beta_0^*). \end{align*} Thus, $\beta_0^*$ is the $\alpha$-quantile of $Y(1)$ and $\beta_1^*$ is the associated expected shortfall. Expected shortfall and related tail-risk functionals are widely used in policy learning in econometrics and related areas; see, e.g., patton2019dynamic, fan2025policy. \end{example} Up to this point, all the examples have set $\omega (\cdot) \equiv \mathsf{p}_{T} (\cdot)$. In the next example, $\omega (\cdot)$ is taken as a known function of $T$. \begin{example}[Stochastic Intervention Effect] The weighting function $\omega (t)$ can be further generalized to depend on the covariates, i.e., $\omega (t, x)$. When $\omega (t, x)$ is a known, user-specified density function of $T$, $\bm{\beta}^{\ast}$ corresponds to the stochastic (dynamical) intervention effect kennedy2019nonparametric. \end{example} \begin{example}[Marginal Structural Models] In its most common formulation, marginal structural models (MSMs) hernan2000marginal postulate a parametric model for the marginal mean of potential outcomes: \begin{equation} {\mathbb{E}} Y (t) = g (t; \bm{\beta}^{\ast}). \end{equation} Model (ref) can be rewritten as Model (ref) by letting $\Gamma (Y (t), t, \bm{\beta}^{\ast}) = Y (t) - g (t; \bm{\beta}^{\ast})$. The weighting function $\omega$ needs to be arbitrary if (ref) is taken as part of the modeling assumptions. See Supplementary Material Section (ref) for further discussions. \end{example} If we set the weighting function $\omega (\cdot)$ in (ref) as the product of the marginal density function of $T$ and the Dirac delta function $\mathsf{p}_T(\cdot)\delta_t (\cdot)$ for each possible treatment value $t \in {\mathcal{T}}$, then the target parameter $\bm{\beta}^* \equiv \beta^{\ast} (\cdot)\in \mathbb{R}^{\infty}$ becomes a function of $t$, that is, an infinite-dimensional ($p=\infty$) parameter. In Section (ref), we further develop our approach for such infinite-dimensional GTM parameters. Notably, this framework encompasses two important infinite-dimensional parameters: the average and quantile dose–response functions (ADRF and QDRF). \begin{example}[Average Dose-Response Function (ADRF)] When $\Gamma (Y(t),t,\beta) = Y(t) - \beta$ and $\omega (\cdot) = \mathsf{p}_T(\cdot)\delta_t (\cdot)$, the solution to (ref) is $\beta^*(t) = {\mathbb{E}} [ Y(t)] $. The ADRF has been widely studied in the literature; see kennedy2017non, su2019non and colangelo2026double, among others. \end{example} \begin{example}[Quantile Dose-Response Function (QDRF)] When $\Gamma(Y(t),t,\beta) = \tau - \mathbbm{1} \{Y(t) \leq \beta\}$ and $\omega (\cdot) = \mathsf{p}_T(\cdot)\delta_t (\cdot)$, the solution to (ref) is $\beta^*(t) = \inf\{q : \mathbb{P} (Y(t)\leq q) \geq \tau\}$, the $\tau^{th}$ quantile of the potential outcome $Y(t)$. The QDRF has been studied by galvao2015uniformly and su2019non. \end{example} \begin{remark} Finally, though not our main focus, we note that other unknown weighting functions $\omega$ can also be considered. For example, when $T \in \mathcal{T} = \mathbb{R}$ is a continuous treatment variable with probability density function $\mathsf{p}_T(t)$, let $\Gamma(Y(t), t, \beta_0) = Y(t) + \beta_0 t$ and $\omega(t) = \partial_t \mathsf{p}_T(t)$, the derivative of $\mathsf{p}_T(t)$. Under Assumption (ref) and using integration by parts, the solution to (ref) becomes $\beta_0^* =\int_{{\mathcal{T}}} \partial_t {\mathbb{E}} \{\mathsf{p}_{T \mid X}^{-1}(T \mid X)Y \mid T=t\} \mathsf{p}_T^2(t) {\mathrm d} t$. This is the density-weighted average derivative effect (DWADE), which has been studied in recent works such as cattaneo2025higher under a randomized experiment, i.e., with a known propensity score function $\mathsf{p}_{T \mid X}$; see Remark (ref) for further discussion. Extending the analysis to observational studies is a topic worth pursuing in future work. \end{remark} \subsection{Review of existing results based on influence functions} As mentioned, we take $\omega (\cdot) \equiv \mathsf{p}_{T} (\cdot)$ here and in Section (ref), so that $\xi (x,t) = \mathsf{p}_{T}(t)/ \mathsf{p}_{T|X} (t|x)$. Before presenting our new results, we first review state-of-the-art estimators of $\bm{\beta}^{\ast}$ based on influence functions from the (first-order) semiparametric efficiency theory. For convenience, we define $b_{\bm{\beta}} (x, t) \coloneqq {\mathbb{E}} \{\Gamma (Y, T, \bm{\beta}) | X = x, T = t\}$. Since $\Gamma (Y, T, \bm{\beta})$ essentially plays the same role as the outcome $Y$ in ATE, we refer to $b_{\bm{\beta}} (\cdot,\cdot)$ as the \emph{generalized outcome regression model}. We denote by $\theta_{\bm{\beta}} \coloneqq (\xi, b_{\bm{\beta}})$ the collection of nuisance parameters for a fixed $\bm{\beta}$. It is worth noting that $\bm{\beta}$ may be non-separable from the nuisance function $b_{\bm{\beta}}$ and can enter in a nonlinear manner. From (ref) and by the law of iterated expectations, $\bm{\beta}^{\ast}$ solves the equation \begin{equation} \psi (\bm{\beta}) \equiv \psi (\theta_{\bm{\beta}}) \coloneqq {\mathbb{E}} \{\xi (X, T) \Gamma (Y, T, \bm{\beta})\} = {\mathbb{E}} \{\xi (X, T) b_{\bm{\beta}} (X, T)\} = 0. \end{equation} ai2021unified developed the semiparametric efficiency theory for $\bm{\beta}^*$. Since the rest of this paper relies on this preliminary result, we reproduce it below in Proposition (ref): \begin{proposition} Under Assumption (ref), the efficient influence function (EIF) of $\psi (\bm{\beta})$ is \begin{align*} \mathsf{IF}_{\psi}^{(1)}(\bm{\beta}) = & \ \xi (X, T) \left\{ \Gamma (Y, T, \bm{\beta}) - b_{\bm{\beta}} (X, T) \right\} + \left\{ \int_{ {\mathcal{T}}} b_{\bm{\beta}} (X, t) \mathsf{p}_{T} (t) {\mathrm d} t - {\mathbb{E}} \{\xi (X, T) b_{\bm{\beta}} (X, T)\} \right\} \notag \\ & + \left\{ \int_{{\mathcal{X}}} b_{\bm{\beta}} (x, T) \mathsf{p}_{X} (x) {\mathrm d} x - {\mathbb{E}} \{\xi (X, T) b_{\bm{\beta}} (X, T)\} \right\}. \end{align*} Furthermore, the EIF of $\bm{\beta}^{\ast}$ that solves $\psi (\bm{\beta}) \equiv 0$ is \begin{equation} \mathsf{IF}_{\bm{\beta}^*}^{(1)} = H_{\bm{\beta}^{\ast}}^{-1} \cdot \mathsf{IF}_{\psi}^{(1)}(\bm{\beta}^{\ast}), \text{ where } H_{\bm{\beta}} \coloneqq - \nabla_{\bm{\beta}} \psi (\bm{\beta}) = - \nabla_{\bm{\beta}} {\mathbb{E}} \{\xi (X, T) \Gamma (Y, T, \bm{\beta})\}. \end{equation} \end{proposition} For certain results in this paper, we impose the following \text{H\"{o}lder} smoothness assumptions on the nuisance parameters $\xi$ and $b_{\bm{\beta}}$. \begin{assumption} $\xi (\cdot, \cdot)$ and $b_{\bm{\beta}} (\cdot, \cdot)$ are assumed to be \text{H\"{o}lder} smooth functions with smoothness indices $s_{1}$ and $s_{2}$ respectively, for every $\bm{\beta} \in {\mathcal{B}}$. We let $s \coloneqq (s_{1} + s_{2}) / 2$ denote the average smoothness of the two nuisance parameters. The nuisance estimators $\widehat{\xi}$ and $\widehat{b}_{\bm{\beta}}$ converge to $\xi$ and $b_{\bm{\beta}}$ respectively at the minimax optimal rates ($n^{- \frac{s_{1}}{d + 2 s_{1}}}$ and $n^{- \frac{s_{2}}{d + 2 s_{2}}}$) in $L_{2} (\mathbb{P})$-distance. \end{assumption} The \text{H\"{o}lder} smoothness conditions in Assumption (ref) are standard for quantifying the complexity of nuisance functions in nonparametric and semiparametric theory. The stated minimax rates are attainable by classical series, sieve, and kernel estimators under appropriate choice of tuning parameters such as the bandwidth or the number of basis functions; see, e.g., stone1982optimal. Assumption (ref) is used to translate the high-level rate conditions on the nuisance estimators into primitive smoothness requirements; see the discussion following Proposition (ref), Remark (ref), and the discussion below Assumption (ref). In the binary QTE example in Section (ref), this condition is specialized as Assumption (ref) and is used in the proof of Theorem (ref). As noted in the Introduction, a major statistical challenge in applying $\mathsf{IF}_{\psi}^{(1)}$ to construct a first-order estimator of $\bm{\beta}$ is that the nuisance parameter $b_{\bm{\beta}}$ itself depends on the unknown parameter $\bm{\beta}$. Directly applying Proposition (ref) therefore requires estimating the entire nuisance function, meaning that $b_{\bm{\beta}}$ needs to be estimated for every candidate $\bm{\beta}$. One simple way to bypass this difficulty is to first obtain an initial estimator $\bm{\beta}_{\mathsf{init}}$ close to $\bm{\beta}^{\ast}$, for example, by solving the empirical version of (ref) and then substitute $b_{\bm{\beta}_{\mathsf{init}}}$ into the influence function $\mathsf{IF}_{\psi}^{(1)}$. This approach is referred to as localized debiased machine learning (LDML), recently developed in kallus2024localized. To formally construct this first-order estimator $\widehat{\bm{\beta}}^{(1)}$, we define \begin{align*} \widehat{\Xi} (O, O'; \bm{\beta}) \coloneqq \widehat{\xi} (X, T) \{\Gamma (Y ,T , \bm{\beta}) - \widehat{b}_{\bm{\beta}_{\mathsf{init}}} (X, T)\} + \widehat{b}_{\bm{\beta}_{\mathsf{init}}} (X, T') , \end{align*} where $O'=(Y',T',X')$ is an independent copy of $O = (Y,T,X)$. We construct the first-order estimator $\widehat{\psi}^{(1)}(\bm{\beta})$ of $\psi(\bm{\beta})$ by the following second-order $U$-statistic (similar construction also appeared in tchetgen2010doubly, bonvini2022sensitivity): \begin{align} \widehat{\psi}^{(1)}(\bm{\beta}) \coloneqq {\mathbb{U}}_{n, 2} \{\widehat{\Xi} (O_{1}, O_{2}; \bm{\beta})\} =\frac{1}{n(n-1)} \sum_{i = 1}^{n}\sum_{i \neq j, j= 1}^{n} \widehat{\Xi} (O_{i}, O_{j}; \bm{\beta}) , \end{align} and define $\widehat{\bm{\beta}}^{(1)}$ as the solution to $\widehat{\psi}^{(1)}(\bm{\beta}) \equiv 0$. For convenience, we assume that the nuisance estimators $\widehat{\xi}$, $\widehat{b}_{\bm{\beta}_{\mathsf{init}}}$ and the initial estimator $\bm{\beta}_{\mathsf{init}}$ of $\bm{\beta}^*$ are all computed using an independent separate sample of size $n$, referred to as the \emph{nuisance sample}. The $U$-statistic operator ${\mathbb{U}}_{n, m}$, for any $m \geq 2$, is taken only over the main sample throughout the paper. All results are stated conditional on the nuisance sample. For brevity, we often suppress this dependence in our notation, particularly the dependence on $\bm{\beta}_{\mathsf{init}}$. To improve efficiency, one can follow the standard cross-fitting procedure of chernozhukov2018double by flipping the roles of the main and nuisance samples to construct a cross-fit version of the estimator $\widehat{\psi}^{(1)} (\bm{\beta})$. However, to avoid notation clutter, we will stick to the estimator without cross-fitting. \begin{remark} In the case of a binary treatment $T \in \{0,1\}$, we have $\xi(X,T) = T/\mathbb{P}(T = 1 | X) + (1 - T)/\mathbb{P}(T = 0 | X)$, $\bm{\beta}=(\beta_0,\beta_1)^{\top}$ and $b_{\bm{\beta}} (X, T) = T \cdot {\mathbb{E}} \{\Gamma(Y,T,\beta_1) \mid X, T=1\} + (1 - T) \cdot {\mathbb{E}} \{\Gamma(Y,T,\beta_0) \mid X, T = 0\}$. The influence function of $\mathsf{IF}_{\psi}^{(1)}(\bm{\beta})$ then reduces to \begin{align*} \mathsf{IF}_{\psi}^{(1)}(\bm{\beta}) = \left( \begin{matrix} T \cdot \xi(X,T) \left\{ \Gamma(Y,1,\beta_1) - b_{\bm{\beta}} (X, T = 1) \right\} + b_{\bm{\beta}} (X, T=1) \\ (1-T)\cdot \xi(X,T) \left\{ \Gamma(Y,0,\beta_0) - b_{\bm{\beta}} (X, T=0) \right\} + b_{\bm{\beta}} (X, T = 0) \end{matrix} \right). \end{align*} This expression was used to construct the first-order LDML estimating equation for QTE, as proposed in Section 1.1 of kallus2024localized; see Section (ref) for a further discussion. For a continuous treatment $T$, the influence function leads to a second-order $U$-statistic estimator. \end{remark} The statistical properties of $\widehat{\bm{\beta}}^{(1)}$ are formally stated in Proposition (ref) below. This result is analogous to Theorem 1 of kallus2024localized. A proof sketch is provided in Section (ref), which also serves to motivate the introduction of our higher-order estimators. \begin{proposition} Suppose that Assumption (ref) holds and there exist three diminishing sequences $r_{n, \xi}, r_{n, b}, r_{n, \bm{\beta}_{\mathsf{init}}}\rightarrow 0$ as $n \rightarrow \infty$ such that \begin{equation*} \begin{split} \Vert \widehat{\xi} - \xi \Vert_{\mathbb{P},2} \lesssim r_{n, \xi},\ \Vert \widehat{b}_{\bm{\beta}_{\mathsf{init}}} - b_{\bm{\beta}_{\mathsf{init}}} \Vert_{\mathbb{P},2} \lesssim r_{n, b} \text{ and } \Vert \bm{\beta}_{\mathsf{init}} - \bm{\beta} \Vert \lesssim r_{n, \bm{\beta}_{\mathsf{init}}}. \end{split} \end{equation*} If we further assume that conditions (i)--(vi) of Theorem 3 in kallus2024localized hold after the substitutions $\theta\mapsto\bm{\beta},\; \mu\mapsto b,\; U+V\mapsto\psi$ and that $r_{n, \xi} ( r_{n, b} + r_{n, \bm{\beta}_{\mathsf{init}}}) = o (n^{- 1 / 2})$, then \begin{align*} \sqrt{n} (\widehat{\bm{\beta}}^{(1)} - \bm{\beta}^{\ast}) \overset{d}{\rightarrow} N (0, \sigma^{2}), \end{align*} where $\sigma^{2} = \mathrm{var} (\mathsf{IF}_{\bm{\beta}^*}^{(1)})$ and $\mathsf{IF}_{\bm{\beta}^*}^{(1)}$ is the EIF of $\bm{\beta}^*$ defined in (ref). \end{proposition} Proposition (ref) shows that, for $\widehat{\bm{\beta}}^{(1)}$ to be $n^{-1/2}$-consistent and asymptotically normal, the initial estimator $\bm{\beta}_{\mathsf{init}}$ must converge to $\bm{\beta}^{\ast}$ fast enough to accurately estimate the nuisance function $b_{\bm{\beta}}$. If $\bm{\beta}_{\mathsf{init}}$ is obtained by solving the empirical version of (ref) using estimated stabilized weights $\widehat{\xi}$ that converge to the true $\xi$ at rate $r_{n,\xi}$, then $r_{n,\bm{\beta}_{\mathsf{init}}} = O(r_{n,\xi})$. Thus, unless $r_{n,\bm{\beta}_{\mathsf{init}}}$ converges to zero at an even faster rate than $r_{n,\xi}$ (a highly unlikely scenario), we require both $r_{n,\xi}$ and $r_{n,b}$ to be $o_{\mathbb{P}}(n^{-1/4})$. Under Assumption (ref), these rate conditions imply that $\frac{2s_{1}/d}{1 + 2 s_{1}/d} > 0.5$ and $\frac{s_{1}/d}{1 + 2 s_{1}/d} + \frac{s_{2}/d}{1 + 2 s_{2}/d} > 0.5$. In the special case $s_{1}=s_{2}=s$, this further reduces to $s/d > 0.5$, under which both nuisance parameters even satisfy the Donsker condition. However, inspired by the known optimal rate for ATE estimation under \text{H\"{o}lder}-type nuisance assumptions robins2023minimax, it is reasonable to conjecture that $s/d \geq 0.25$ suffices to guarantee a $n^{-1/2}$‐consistent estimator for finite-dimensional parameters of GTM under reasonable regularity conditions on the generalized residual function $\Gamma$. The remainder of this section is devoted to constructing a higher-order estimator of $\bm{\beta}^{\ast}$ requiring much relaxed \text{H\"{o}lder}-type smoothness assumptions on the nuisance parameters and the initial $\bm{\beta}_{\mathsf{init}}$, by leveraging the HOIF framework. \section{Higher-Order Estimators: The Finite-Dimensional Case} In this section, we present our main methodological and theoretical results on new higher-order estimators for finite-dimensional GTM parameters. We first motivate the need for higher-order estimators beyond the first-order estimator $\widehat{\bm{\beta}}^{(1)}$ through a simple numerical experiment in Section (ref), accompanied by a brief and intuitive introduction to HOIF theory. Section (ref) then discusses how to construct higher-order estimators for GTM parameters, and Section (ref) studies their statistical properties. Finally, Section (ref) specializes these general results to the QTE. \subsection{Motivation for higher-order estimators: The QTE example} Before proceeding, we provide a more technically rigorous motivation for constructing higher-order estimators for GTMs, complementing the empirical motivation based on the simulation results in Figure (ref) in the Introduction. We focus on $\beta_{1}^{\ast}$, the $\tau$-quantile of the potential outcome $Y(1)$, which was shown to be a GTM parameter in Example (ref). As discussed in Remark (ref) and Proposition (ref) (following kallus2024localized), we can analyze the bias of the LDML estimator $\widehat{\beta}^{(1)}$ by examining the bias of the first-order moment equation $\widehat{\psi}^{(1)}(\beta_1)$ for $\psi(\beta_1)$ at a given $\beta_1$: \begin{equation*} \begin{split} \mathbb{E} \{\widehat{\psi}^{(1)} (\beta_{1}) - \psi (\beta_{1})\} & = \mathbb{E} \left[ \left\{ \frac{\mathbb{P} (T = 1 \mid X)}{\widehat{\mathbb{P}} (T = 1 \mid X)} - 1 \right\} \{b_{\beta_{1}} (X, T = 1) - \widehat{b}_{\beta_{1, \mathsf{init}}} (X, T = 1)\} \right] \\ & = \mathbb{E} \left[ \{\widehat{\xi} (X, T) - \xi (X, T)\} \{b_{\beta_{1}} (X, T) - \widehat{b}_{\beta_{1, \mathsf{init}}} (X, T)\} \right], \end{split} \end{equation*} where we recall that $\xi (x, t) = t / \mathbb{P} (T = 1 \mid X = x)$, $b_{\beta_{1}} (x, t) = t \cdot \mathbb{E} [\tau - \mathbbm{1} \{Y \leq \beta_{1}\} \mid X = x, T = t]$, and $\widehat{\xi}$ and $\widehat{b}_{\beta_{1}}$ are their respective estimates. Recall that $\widehat{\xi}$ and $\widehat{b}_{\beta_{1,\mathsf{init}}}$ are constructed on the nuisance sample, whereas the expectation is taken with respect to the main sample. Hence, conditional on the nuisance sample, $\widehat{\xi}$ and $\widehat{b}_{\beta_{1,\mathsf{init}}}$ can be treated as fixed functions, and the bias $\mathbb{E} [ \{\widehat{\xi} (X, T) - \xi (X, T)\} \{b_{\beta_{1}} (X, T) - \widehat{b}_{\beta_{1, \mathsf{init}}} (X, T)\} ]$ captures only the approximation bias of these fixed estimators, not the sampling variability in the nuisance estimates. Decomposing $b_{\beta_{1}} - \widehat{b}_{\beta_{1, \mathsf{init}}}$ as $(b_{\beta_{1}} - b_{\beta_{1, \mathsf{init}}}) + (b_{\beta_{1, \mathsf{init}}} - \widehat{b}_{\beta_{1, \mathsf{init}}})$ and applying the triangle inequality, we can bound the bias of $\widehat{\psi}^{(1)}(\beta_1)$ under suitable regularity conditions on $\beta_1 \mapsto b_{\beta_1}$ as follows: \begin{align} \Vert \widehat{\xi} - \xi \Vert_{\mathbb{P}, 2} (\Vert \widehat{b}_{\beta_{1, \mathsf{init}}} - b_{\beta_{1, \mathsf{init}}} \Vert_{\mathbb{P}, 2} + \Vert \beta_{1, \mathsf{init}} - \beta_{1} \Vert). \end{align} Since $\beta_{1, \mathsf{init}}$ depends on $\widehat{\xi}$, the bias of $\widehat{\psi}^{(1)}(\beta_1)$ is consequently bounded by \begin{align} \Vert \widehat{\xi} - \xi \Vert_{\mathbb{P}, 2} \Vert \widehat{b}_{\beta_{1, \mathsf{init}}} - b_{\beta_{1, \mathsf{init}}} \Vert_{\mathbb{P}, 2} + \Vert \widehat{\xi} - \xi \Vert_{\mathbb{P}, 2}^{2}. \end{align} If $\widehat{\xi}$ or $\widehat{b}_{\beta_{1, \mathsf{init}}}$ does not converge sufficiently fast to their true values $\xi$ or $b_{\beta_{1, \mathsf{init}}}$, the bias of $\widehat{\psi}^{(1)} (\beta_{1})$ might not be sufficiently small, as indicated by the bounds (ref) and (ref). The higher-order estimator $\widehat{\beta}_{1}^{(2)}$ is motivated by the HOIF theory initially developed in robins2008higher, which is essentially a bias reduction technique. To apply the HOIF theory to estimate $\beta_{1}^{\ast}$, we first fix a set of basis functions of $(X, T)$ of dimension $k$, denoted by $\bar{\phi}_{k}$. Since $T$ is binary, we have $\bar{\phi}_{k} (x, t) = t \bar{{\mathsf{z}}}_{k} (x)$ for some basis functions $\bar{{\mathsf{z}}}_{k}$ over $X$. Given $\bar{\phi}_{k}$, we introduce the following projections of the residuals $\widehat{\xi} - \xi$ and $\widehat{b}_{\beta_{1, \mathsf{init}}} - b_{\beta_{1}}$ onto the space spanned by $\phi_{k}$ (defined in Section (ref)): \begin{align*} & \Pi (\widehat{\xi} - \xi \mid \bar{\phi}_{k}) (x, t) = \bar{\phi}_{k} (x, t)^{\top} \Sigma_{k}^{-1} \mathbb{E} [\{\widehat{\xi} (X, T) - \xi (X, T)\} \bar{\phi}_{k} (X, T)], \\ & \Pi (\widehat{b}_{\beta_{1, \mathsf{init}}} - b_{\beta_{1}} \mid \bar{\phi}_{k}) (x, t) = \bar{\phi}_{k} (x, t)^{\top} \Sigma_{k}^{-1} \mathbb{E} [\{\widehat{b}_{\beta_{1, \mathsf{init}}} (X, T) - b_{\beta_{1}} (X, T)\} \bar{\phi}_{k} (X, T)], \end{align*} where $\Sigma_{k} = \mathbb{E} \{\bar{\phi}_{k} (X, T)^{\otimes 2}\} = \mathbb{E} \{T \bar{{\mathsf{z}}}_{k} (X)^{\otimes 2}\}$. Although it is generally impossible to estimate the bias of an estimator without imposing strong assumptions, the following term can be estimated unbiasedly by a second-order $U$-statistic computable from the data (if $\Sigma_{k}$ is not known, we can plug-in some estimator $\widehat{\Sigma}_{k}$ of $\Sigma_{k}$): \begin{align*} \mathrm{B}_{\psi, k} (\beta_{1}) & \coloneqq \mathbb{E} \{\Pi (\widehat{\xi} - \xi \mid \bar{\phi}_{k}) (X, T) \cdot \Pi (b_{\beta_{1}} - \widehat{b}_{\beta_{1, \mathsf{init}}} \mid \bar{\phi}_{k}) (X, T)\} \\ & = \mathbb{E} [\{\widehat{\xi} (X, T) - \xi (X, T)\} \bar{\phi}_{k} (X, T)^{\top}] \Sigma_{k}^{-1} \mathbb{E} [\{b_{\beta_{1}} (X, T) - \widehat{b}_{\beta_{1, \mathsf{init}}} (X, T)\} \bar{\phi}_{k} (X, T)] \\ & = \mathbb{E} [\{\widehat{\xi} (X, T) - 1\} \bar{{\mathsf{z}}}_{k} (X)^{\top}] \Sigma_{k}^{-1} \mathbb{E} [\{T (\tau - \mathbbm{1} \{Y \leq \beta_{1}\}) - \widehat{b}_{\beta_{1, \mathsf{init}}} (X, T)\} \bar{{\mathsf{z}}}_{k} (X)], \end{align*} as it is a product of two expectations of random variables that can be directly evaluated from the data. The last equality follows from $\bar{\phi}_k(x,t) = t \bar{{\mathsf{z}}}_k(x)$ and $\xi(x,t) = t / \mathbb{P} (T = 1 \mid X = x)$, which give ${\mathbb{E}} \{\xi(X,T) \bar{\phi}_k(X,T)\} = {\mathbb{E}} \{\bar{{\mathsf{z}}}_k(X)\}$, and by the law of iterated expectations: \begin{align*} \mathbb{E} \{b_{\beta_{1}} (X, T) \bar{\phi}_{k} (X, T)\} = \mathbb{E} \{T \cdot \mathbb{E} (\tau - \mathbbm{1} \{Y \leq \beta_{1}\} \mid X, T) \cdot \bar{{\mathsf{z}}}_k(X)\} = \mathbb{E} \{T (\tau - \mathbbm{1} \{Y \leq \beta_{1}\}) \bar{{\mathsf{z}}}_{k} (X)\}. \end{align*} Furthermore, the bias of $\widehat{\psi}^{(1)} (\beta_{1})$ can be decomposed as \begin{align*} \mathbb{E} \{\widehat{\psi}^{(1)} (\beta_{1}) - \psi (\beta_{1})\} = \mathrm{B}_{\psi, k} (\beta_{1}) + \mathbb{E} \{\Pi^{\perp} (\widehat{\xi} - \xi \mid \bar{\phi}_{k}) (X, T) \cdot \Pi^{\perp} (b_{\beta_{1}} - \widehat{b}_{\beta_{1, \mathsf{init}}} \mid \bar{\phi}_{k}) (X, T)\}, \end{align*} where $\Pi^{\perp}$ denotes the orthocomplement of the projection and following robins2008higher, we denote the second term in the decomposition as $\mathrm{TB}_{\psi, k} (\beta_{1})$. $\mathrm{TB}_{\psi, k} (\beta_{1})$ is referred to as the \emph{truncation bias} in robins2008higher, because it represents the remaining bias, after “truncating” the potentially infinite-dimensional residual functions $\widehat{\xi} - \xi$ and $b_{\beta_{1}} - \widehat{b}_{\beta_{1, \mathsf{init}}}$ up to a finite dimension $k$ by the basis $\bar{\phi}_{k}$. Under appropriate assumptions on $\widehat{\xi} - \xi$ and $b_{\beta_{1}} - \widehat{b}_{\beta_{1, \mathsf{init}}}$ and a suitable choice of $\bar{\phi}_{k}$ (to be clarified later), $\mathrm{TB}_{\psi, k} (\beta_{1})$ can be made negligible compared to sampling variability. The key idea of our higher-order estimation approach is to decompose the bias of the first-order estimator $\widehat{\psi}^{(1)} (\beta_{1})$ of $\psi (\beta_{1})$ into two parts $\mathrm{B}_{\psi, k} (\beta_{1}) + \mathrm{TB}_{\psi, k} (\beta_{1})$, so that $\mathrm{B}_{\psi, k} (\beta_{1})$ can be estimated at a sufficiently fast rate by an estimator $\widehat{\mathrm{B}}_{\psi, k}$ and $\mathrm{TB}_{\psi, k} (\beta_{1})$ is negligible. We then solve a bias-reduced estimating equation $\widehat{\psi}^{(1)} - \widehat{\mathrm{B}}_{\psi, k}$ (typically a $U$-statistic) to construct the higher-order estimator $\widehat{\beta}_{1}^{(2)}$. We summarize the construction of $\widehat{\beta}_{1}^{(2)}$ in Figure (ref). \begin{figure} \begin{tikzpicture}[ >=Latex, font=, node distance=1.8cm and 2.0cm, box/.style={ draw, semithick, rectangle, minimum width=2.6cm, minimum height=1.0cm, align=center }, plug/.style={->,semithick}, est/.style={ ->, semithick, decorate, decoration={snake, amplitude=0.9pt, segment length=5pt}, shorten >=-1.1pt }, every node/.style={inner sep=2pt} ] \node[box] (nuis) at (0,1) {nuisance sample}; \node[box, below=0.5cm of nuis] (main) {main sample}; \node[left=0.8cm of nuis] (xi) {$\widehat{\xi}$}; \node[above=0.5cm of nuis] (betainit) {$\beta_{1,\mathrm{init}}$}; \node[right=1.2cm of nuis.north east, anchor=west, yshift=-0.1cm] (bhat) {$\widehat{b}_{\beta_{1,\mathrm{init}}}$}; \node[right=0.8cm of nuis.south east, anchor=west,, yshift=0.1cm] (sigma) {$\widehat{\Sigma}_k$}; \node[ right=2.2cm of main] (grp) {$ \left\{ \begin{array}{c} \widehat{\psi}^{(1)}\\ \widehat{\mathrm{B}}_{\psi,k} \end{array} \right\} $}; \node[right=1.3cm of grp] (phi2) {$\widehat{\psi}^{(2)}_k = \widehat{\psi}^{(1)} - \widehat{\mathrm{B}}_{\psi,k}$}; \node[right=1.3cm of phi2] (beta2) {$\widehat{\beta}^{(2)}_1$}; \draw[est] (nuis.west) -- (xi.east); \draw[est] (nuis.north) -- (betainit.south); \draw[plug] (xi.north) |- (betainit.west); \draw[est] ([, yshift=-0.1cm]nuis.north east) -- (bhat.west); \draw[est] ([, yshift=0.1cm]nuis.south east) -- (sigma.west); \draw[plug] (betainit.east) -| ([xshift=-0.8cm]bhat.north); \draw[plug] (bhat.south) --([yshift=-1.6cm]bhat.south); \draw[plug] (sigma.south) --([yshift=-0.8cm]sigma.south); \draw[plug] ([yshift=-0.1cm]xi.south) -- ([yshift=-2.4cm]xi.south) -- ([yshift=-2.8cm]bhat.south) -- ([yshift=-1.7cm]bhat.south); \draw[plug] (main.east) -- (grp.west); \draw[plug] (grp.east) -- (phi2.west); \node[right=0.4cm of phi2] (arr) {$\Longrightarrow$}; \end{tikzpicture} \caption{A schematic illustration on how to construct the higher-order estimator $\widehat{\beta}_{1}^{(2)}$ of $\beta_{1}^*$ in the example of QTE. $\widehat{\mathrm{B}}_{\psi, k}$ is a second-order $U$-statistic estimator of $\mathrm{B}_{\psi, k}$.} \end{figure} As mentioned in the Introduction, Figure (ref) displays the histograms of different QTE estimators after centered on the true QTE value (i.e., the estimate minus the truth) in a simulation setting where nuisance estimation errors $\|\widehat{\xi} - \xi\|_{\mathbb{P},2}$ and $\|\widehat{b}_{\beta_{1,\mathsf{init}}} - b_{\beta_{1,\mathsf{init}}}\|_{\mathbb{P},2}$ are designed to be large. In the left panel, the LDML estimator exhibits a noticeable bias relative to its sampling variability. In contrast, the right panel shows that the histogram of our higher-order estimators is much more centered around zero, indicating a negligible bias relative to its sampling variability. The variance of the higher-order estimator is only slightly larger than that of the LDML estimator. These results demonstrate that the proposed higher-order estimator can be useful not only in theory but also in practice. \subsection{Higher-order estimators: The general case} In the general case, we begin by noting that the bias of the first-order estimator $\widehat{\psi}^{(1)} (\bm{\beta})$ for the moment equation $\psi(\bm{\beta})$ can also be written as the mean of a product of two nuisance estimation residuals: \begin{align} {\mathbb{E}} \{\widehat{\psi}^{(1)}(\bm{\beta}) - \psi(\bm{\beta})\} = {\mathbb{E}} [\{\widehat{\xi} (X, T) - \xi (X, T)\} \{b_{\bm{\beta}} (X, T) - \widehat{b}_{\bm{\beta}_{\mathsf{init}}} (X, T)\}]. \end{align} To reduce this bias, we project each residual onto a size-$k$ dictionary $\bar{\phi}_{k} \coloneqq \{\phi_{1}, \cdots, \phi_{k}\}$ functions. Specifically, we select $\bar{\phi}_{k} \equiv \bar{{\mathsf{z}}}_{k_{x}} \otimes \bar{{\mathsf{w}}}_{k_{t}}$ as the tensor product of a size-$k_{x}$ dictionary $\bar{{\mathsf{z}}}_{k_{x}} \coloneqq \{\mathsf{z}_{1}, \cdots, \mathsf{z}_{k_{x}}\}$ over ${\mathcal{X}}$ and a size-$k_{t}$ dictionary $\bar{{\mathsf{w}}}_{k_{t}} \coloneqq \{\mathsf{w}_{1}, \cdots, \mathsf{w}_{k_{t}}\}$ over ${\mathcal{T}}$, as in general $T$ is no longer binary. Using the same projection arguments as in Section (ref), the bias decomposes as \begin{align*} {\mathbb{E}} \{\widehat{\psi}^{(1)} (\bm{\beta}) - \psi (\bm{\beta})\} & = {\mathbb{E}} \{\Pi (\widehat{\xi} - \xi \mid \bar{\phi}_{k}) (X, T) \cdot \Pi (b_{\bm{\beta}} - \widehat{b}_{\bm{\beta}_{\mathsf{init}}} \mid \bar{\phi}_{k}) (X, T)\} \\ &\quad + {\mathbb{E}} \{\Pi^{\perp} (\widehat{\xi} - \xi \mid \bar{\phi}_{k}) (X, T) \cdot \Pi^{\perp} (b_{\bm{\beta}} - \widehat{b}_{\bm{\beta}_{\mathsf{init}}} \mid \bar{\phi}_{k}) (X, T)\} \\ & \eqqcolon \mathrm{B}_{\psi, k} (\bm{\beta})+ \mathrm{TB}_{\psi, k} (\bm{\beta}), \end{align*} where $\Pi$ denotes projection onto the span of $\bar{\phi}_k$ and $\Pi^{\perp}$ is its orthocomplement. Moreover, \begin{equation} \mathrm{B}_{\psi, k} (\bm{\beta}) = {\mathbb{E}} [\{\widehat{\xi} (X, T) - \xi (X, T)\} \bar{\phi}_{k} (X, T)^{\top}] \Sigma_{k}^{-1} {\mathbb{E}} [\bar{\phi}_{k} (X, T) \{b_{\bm{\beta}} (X, T) - \widehat{b}_{\bm{\beta}_{\mathsf{init}}} (X, T)\}]. \end{equation} A key insight underlying the improved second-order estimator is that the bias term $\mathrm{B}_{\psi,k}(\bm{\beta})$ admits an unbiased sample analogue, if $\Sigma_k$ is known. First, by the law of iterated expectations, \begin{align*} \mathbb{E}[\bar{\phi}_k (X,T) \{b_{\bm{\beta}}(X,T)-\widehat{b}_{\bm{\beta}_{\mathsf{init}}}(X,T)\}] = \mathbb{E}[\bar{\phi}_k(X,T)\{\Gamma(Y,T,\bm{\beta})-\widehat{b}_{\bm{\beta}_{\mathsf{init}}}(X,T)\}]. \end{align*} Next, recall that $\xi(X,T)=\mathsf{p}_T(T)/\mathsf{p}_{T\mid X}(T\mid X)$. Weighting by $\xi$ transforms the joint distribution of $(X,T)$ into the product of their marginal distributions. Indeed, \begin{align*} \mathbb{E}[\xi(X,T)\bar{\phi}_k(X,T)] = \iint \bar{\phi}_k(x,t)\,\mathsf{p}_T(t)\,\mathsf{p}_X(x) {\mathrm d} t {\mathrm d} x. \end{align*} This identity is crucial: it shows that the expectation of $\xi(X,T)\bar{\phi}_k(X,T)$ can be estimated by averaging $\bar{\phi}_k(X_i,T_j)$ over independent pairs $(X_i,T_j)$ from the main sample. Combining these two observations, we obtain an unbiased estimator of $\mathrm{B}_{\psi,k}(\bm{\beta})$ as the following third-order $U$-statistic: \begin{align*} \widetilde{\mathrm{B}}_{\psi,k}(\bm{\beta}) = \mathbb{U}_{n,3}\Big[\big\{\widehat{\xi}(X_1,T_1)\bar{\phi}_k(X_1,T_1) - \bar{\phi}_k(X_1,T_3)\big\}^{\top} \Sigma_k^{-1} \bar{\phi}_k(X_2,T_2) \big\{\Gamma(Y_2,T_2,\bm{\beta})-\widehat{b}_{\bm{\beta}_{\mathsf{init}}}(X_2,T_2)\big\}\Big]. \end{align*} Recall that the $U$-statistics operator $\mathbb{U}_{n,3}$ is taken over the main sample. We then construct the improved oracle second-order estimator of $\psi (\bm{\beta})$ as \begin{align*} \widetilde{\psi}_{ k}^{(2)} (\bm{\beta}) \coloneqq \widehat{\psi}^{(1)} (\bm{\beta}) - \widetilde{\mathrm{B}}_{\psi,k} (\bm{\beta}). \end{align*} We call $\widetilde{\psi}_{ k}^{(2)}$ an oracle estimator because it requires the knowledge of $\Sigma_{k}$, which is generally unknown. The corresponding oracle second-order estimator $\widetilde{\bm{\beta}}^{(2)}$ for $\bm{\beta}^*$ is defined as the solution to $\widetilde{\psi}_{k}^{(2)} (\bm{\beta}) = 0$. When $\Sigma_{k}$ is replaced by an estimator $\widehat{\Sigma}_{k}$ (see Remark (ref)) computed from the nuisance sample, we obtain the feasible estimators $\widehat{\mathrm{B}}_{\psi, k}$ and $\widehat{\psi}_{k}^{(2)} (\bm{\beta}) \coloneqq \widehat{\psi}^{(1)} (\bm{\beta}) - \widehat{\mathrm{B}}_{\psi, k}(\bm{\beta})$. The feasible estimator $\widehat{\bm{\beta}}^{(2)}$ for $\bm{\beta}^*$ is then defined as the solution to $\widehat{\psi}_k^{(2)}(\bm{\beta}) = 0$. Let $\bar{\psi}_{k} (\bm{\beta}) \coloneqq {\mathbb{E}} [\widetilde{\psi}_{k}^{(2)} (\bm{\beta})]$ denote the truncated parameter of $\psi (\bm{\beta}) $ robins2016technical. The bias of the second-order estimator $\widetilde{\psi}_{k}^{(2)} (\bm{\beta})$ then reduces to \begin{align*} {\mathbb{E}} \{\widetilde{\psi}_{k}^{(2)} (\bm{\beta}) - \psi (\bm{\beta})\} = \bar{\psi}_{k} (\bm{\beta})- \psi (\bm{\beta}) = \mathrm{TB}_{\psi, k} (\bm{\beta}), \end{align*} which can be made much smaller than the original bias. Remark (ref) discusses why this yields an improved convergence rate for $\widetilde{\psi}_{k}^{(2)} (\bm{\beta})$ under \text{H\"{o}lder}-type assumptions in Assumption (ref). Formally analyzing the statistical properties of $\widetilde{\bm{\beta}}^{(2)}$ and $\widehat{\bm{\beta}}^{(2)}$ requires a detailed examination of the impact of the initial estimator $\bm{\beta}_{\mathsf{init}}$ and heavy use of $U$-processes theory, since they are third-order $Z$-estimators. We will make these intuitive arguments precise in Section (ref). \begin{remark} The initial nuisance estimators $\widehat{\xi}$ and $\widehat{b}_{\bm{\beta}_{\mathsf{init}}}$ can be obtained by any flexible method (e.g., kernel smoothing, series estimation, or machine learning) and are treated as fixed functions conditional on the nuisance sample. The dictionary $\bar{\phi}_k$ is introduced solely for the purpose of bias correction: we project the residual functions $\widehat{\xi}-\xi$ and $b_{\bm{\beta}}-\widehat{b}_{\bm{\beta}_{\mathsf{init}}}$ onto the linear span of $\bar{\phi}_k$, and the bias correction term $\widetilde{\mathrm{B}}_{\psi,k}(\bm{\beta})$ estimates the inner product of these projections. The improvement comes from choosing $\bar{\phi}_k$ with sufficient approximation power for these residuals; in particular, the dictionary used for bias correction should be richer than or different from any dictionary used in the initial nuisance estimation. \end{remark} \begin{remark} We now provide intuition for the improved convergence rate of the debiased estimator $\widetilde{\psi}_{k}^{(2)} (\bm{\beta})$. For simplicity, we ignore the impact of $\bm{\beta}_{\mathsf{init}}$ by assuming $\bm{\beta}_{\mathsf{init}} = \bm{\beta}$. Under Assumption (ref) and with an appropriate choice of basis $\bar{\phi}_{k}$ belloni2015some, the truncation bias satisfies $\mathrm{TB}_{\psi, k}(\bm{\beta}) \asymp k^{- 2 s/d}$, where $s$ is the average smoothness defined in Assumption (ref). The variance of $\widetilde{\psi}_{k}^{(2)} (\bm{\beta})$ is of order $\frac{k}{n^{2}} + \frac{1}{n}$. Balancing the truncation bias against the variance by choosing $k \asymp n^{\frac{2 d}{d + 4 s}}$ yields \begin{align*} \left[ {\mathbb{E}} \{\widetilde{\psi}_{k}^{(2)} (\bm{\beta}) - \psi (\bm{\beta})\}^{2} \right]^{1 / 2} \lesssim \left\{ \begin{array}{ll} n^{-1 / 2} & s \geq \frac{d}{4}, \\ n^{- \frac{4 s}{d + 4 s}} & s < \frac{d}{4}. \end{array} \right. \end{align*} Thus, compared to the first-order estimator $\widehat{\psi}^{(1)}$, the debiased estimator $\widetilde{\psi}_k^{(2)}$ relaxes the smoothness requirement for $\sqrt{n}$-consistency. Specifically, $\widetilde{\psi}_k^{(2)}$ achieves the $\sqrt{n}$ rate whenever $s \ge d/4$, whereas $\widehat{\psi}^{(1)}$ fails in the regime $d/4 < s < d/2$ (see the discussion following Proposition (ref)). \end{remark} \begin{remark} Following mcgrath2024nuisance and chen2024method, we focus on the case where $\Sigma_{k}^{-1}$ is known or some estimator $\widehat{\Sigma}_{k}^{-1}$ of $\Sigma_{k}^{-1}$ is sufficiently close to $\Sigma_{k}^{-1}$, which generally requires extra smoothness assumptions on the marginal density $\mathsf{p}_{X}$ of $X$. How the minimax-optimal rate depends on the smoothness of $\mathsf{p}_{X}$ is a long-standing open problem richardson2014causal and is beyond the scope of this paper. For completeness, in Supplementary Material Section (ref), we derive the HOIFs of $\psi (\bm{\beta})$ (or more precisely, of $\bar{\psi}_{k} (\bm{\beta}) $) and construct higher-order estimators of arbitrary order $m$, facilitating future work in this direction. When $k = o (n)$, we briefly discuss in Supplementary Material Section (ref) how to estimate $\bm{\beta}$ using higher-order estimators with $m \asymp \log n$ without any smoothness assumption on $\mathsf{p}_{X}$, where $\Sigma_{k}$ is estimated by its sample average estimator. \end{remark} \subsection{Asymptotic properties of higher-order estimators: The general case} \emph{En route} to deriving statistical properties of our proposed second-order estimators, we further require the following regularity conditions. When stating these conditions, we introduce positive and finite constants $c_{3}, \ldots, c_{8}$ that do not vary with $n$. \end{assumption} \fi \begin{assumption} The support $\mathcal{T}$ of the treatment variable $T$ is a compact subset of ${\mathbb{R}}$. \end{assumption} \begin{assumption} (i) The parameter space ${\mathcal{B}}$ is a compact subset of ${\mathbb{R}}^{p}$ and the true parameter $\bm{\beta}^{*}$ is in the interior of ${\mathcal{B}}$. (ii) $\bm{\beta}^{\ast}$ is the unique solution to the population first-order estimating equation $\psi (\bm{\beta}) = 0$ and for every $\epsilon > 0$, $\inf_{\|\bm{\beta} - \bm{\beta}^{\ast}\| \geq \epsilon} \|\psi (\bm{\beta})\| > 0$. \end{assumption} \begin{assumption} (i) $\psi (\bm{\beta})$ is differentiable with respect to $\bm{\beta}$ at any $\bm{\beta} \in {\mathcal{B}}$; denote its derivative with respect to $\bm{\beta}$ as $\nabla_{\bm{\beta}} \psi (\bm{\beta}) \in {\mathbb{R}}^{p \times p}$. (ii) Each component of $\nabla_{\bm{\beta}} \psi (\bm{\beta})$ is continuous at $\bm{\beta}^{\ast}$. (iii) The eigenvalues of $\nabla_{\bm{\beta}} \psi (\bm{\beta}^{\ast})$ are bounded between constants $c_3$ and $c_4$. \end{assumption} \begin{assumption} (i) $\sup_{\bm{\beta} \in {\mathcal{B}}} \|b_{\bm{\beta}} (\cdot, \cdot)\|_{\infty} \lesssim 1$, $\sup_{\bm{\beta} \in {\mathcal{B}}} \|\widehat{b}_{\bm{\beta}} (\cdot, \cdot)\|_{\infty} \lesssim 1$ and $ \|\widehat{\xi} (\cdot, \cdot) \|_{\infty} \lesssim 1$. (ii) For any $\bm{\beta}_{1}, \bm{\beta}_{2} \in{\mathcal{B}}$, we have $ \| {b}_{\bm{\beta}_{1}}(X,T)-{b}_{\bm{\beta}_{2}} (X,T) \|_{\mathbb{P},2} \lesssim \| \bm{\beta}_{1} - \bm{\beta}_{2}\|.$ (iii) For any $(x, t) \in {\mathcal{X}} \times {\mathcal{T}}$, we have $ \ c_5 \leq \mathsf{p}_{X}(x)\leq c_6 $ and $ \ c_5 \leq \mathsf{p}_{T}(t)\leq c_6$. \end{assumption} \begin{assumption} (i) The function class \ $\left\{ \Gamma (\cdot, \cdot, \bm{\beta}): \bm{\beta} \in {\mathcal{B}} \right\}$ is of VC type with fixed characteristic $(A, v)$. (ii) $\sup_{\bm{\beta} \in {\mathcal{B}}} \| \Gamma (\cdot,\cdot,\bm{\beta}) \|_{\infty} \lesssim 1$. (iii) There exists some constant $\alpha_{0}\in(0,1]$ such that, for any $\bm{\beta}_{1}, \bm{\beta}_{2} \in{\mathcal{B}}$, any small $\delta>0$, ${\mathbb{E}}^{1/2} [\sup_{||\bm{\beta}_{1}-\bm{\beta}_{2}||\lesssim\delta} \left\| \Gamma(Y,T,\bm{\beta}_1)-\Gamma(Y,T,\bm{\beta}_2) \right\|^2 ]\lesssim \delta^{\alpha_{0}}$. \end{assumption} \begin{assumption} (i) For each $j$, the $j^{th}$ basis functions $\mathsf{z}_{j}(\cdot)$ and $\mathsf{w}_{j}(\cdot)$ are supported on sets of probability at most $c_7/k_x$ and $c_7/k_t$, respectively. Let ${\mathcal{S}}_{X,j}$ and ${\mathcal{S}}_{T,j}$ denote the supports of $\mathsf{z}_{j}(\cdot)$ and $\mathsf{w}_{j}(\cdot)$, respectively. (ii) At any point $(x,t)\in{\mathcal{X}} \times {\mathcal{T}}$, at most $c_8$ elements of $\bar{{\mathsf{z}}}_{k_{x}} (x)$ and $\bar{{\mathsf{w}}}_{k_{t}} (t)$ are simultaneously non-zero. (iii) For every $k_{x}$ and $k_{t}$, and for any unit vectors ${\bm{a}}=(a_1,\ldots,a_{k_x})\in {\mathbb{R}}^{k_x}$ and ${\bm{b}}=(b_1,\ldots,b_{k_t})\in {\mathbb{R}}^{k_t}$, we have ${\bm{a}}^{\top}\int_{{\mathcal{S}}_{X,j}} \bar{{\mathsf{z}}}_{k_{x}} (x)\bar{{\mathsf{z}}}_{k_{x}} (x)^{\top} {\mathrm d} x {\bm{a}} \gtrsim a_j^2 $ and ${\bm{b}}^{\top}\int_{{\mathcal{S}}_{T,j}} \bar{{\mathsf{w}}}_{k_{t}} (t) \bar{{\mathsf{w}}}_{k_{t}} (t)^{\top} {\mathrm d} t {\bm{b}} \gtrsim b_j^2$. (iv) There exist two sequences of constants $\zeta_{x}(k_{x})$ and $\zeta_{t}(k_{t})$ such that $\Vert \bar{{\mathsf{w}}}_{k_{t}} (\cdot) \Vert_{\infty} \leq \zeta_{t}(k_{t}) \lesssim \sqrt{k_t}$, $\Vert \bar{{\mathsf{z}}}_{k_{x}} (\cdot) \Vert_{\infty} \leq \zeta_{x}(k_{x}) \lesssim \sqrt{k_x}$, and $\zeta(k)\coloneqq\zeta _{x}(k_{x})\zeta_{t}(k_{t})$. (v) The dimension $k=k_x k_t$ satisfies $k \to \infty$ and $k = o (n^2)$ as $n \to \infty$. \end{assumption} Assumption (ref) requires that the treatment $T$ be bounded. This condition is imposed to simplify the analysis, as most existing technical tools for $U$-statistics and $U$-processes require the $U$-statistic kernel to be bounded. We expect that it could be weakened to light-tailed conditions chakrabortty2025tail, but this is not the main focus of our paper. Assumptions (ref) and (ref) are standard identification and regularity conditions for $Z$-estimators. Assumption (ref) imposes sufficient regularity conditions on the nuisance functions and their estimators. In particular, Assumption (ref)(iii), together with Assumption (ref)(ii), implies that $\xi(\cdot,\cdot)$ is bounded away from zero and infinity, which is commonly required in the existing literature; see ai2021unified and kennedy2017non, among others. Assumption (ref)(i) is a high-level complexity condition on the generalized residual function that defines GTM treatment effect parameters in (ref). In applications, this assumption shall be verified case-by-case. Assumptions (ref)(ii)-(iii) control the envelope functions, a standard requirement in $M/Z$-estimation when the criterion is non-smooth (see, e.g., chen2003estimation,ai2021unified). Notably, Assumption (ref)(iii) is weaker than the Lipschitz continuity condition and accommodates various continuous and discontinuous functions. Fortunately, in various examples, these assumptions are straightforward to verify. In Supplementary Material Section (ref), we show that they hold for both ATE and QTE from Examples (ref) and (ref). Assumptions (ref)(i)--(ii) ensure control over the sup norm of series projection estimators, which are essential for analyzing the $U$-process and establishing variance bounds. These assumptions are satisfied by many commonly used sieve bases, including Cohen–Daubechies–Vial wavelet series, B-splines and local polynomial partition series; see belloni2015some, liu2017semiparametric,cattaneo2020large for detailed discussions. Assumption (ref)(iii) imposes a mild local non-collinearity requirement on the basis functions and is verified by the basis functions mentioned above. This condition is also adopted in cattaneo2020large. Assumption (ref)(iv) is standard in the nonparametric regression literature and is verified for the basis functions mentioned above. Assumption (ref)(v) restricts the basis dimension $k$ from growing too fast (slower than $n^{2}$), a necessary condition to achieve asymptotic normality of the estimators. We now establish the consistency and convergence rate of the oracle second-order estimator $\widetilde{\bm{\beta}}^{(2)}$. \begin{theorem} Suppose that Assumptions (ref) and (ref)--(ref) hold. Let $\widetilde{r}_{n,\bm{\beta}} \to 0$ be a diminishing sequence such that \begin{align*} \bigg( \frac{\sqrt{k}}{n} \vee \frac{1}{\sqrt{n}}\bigg)\cdot \log n + \sup_{\bm{\beta}\in{\mathcal{B}}}\| \Pi^{\perp} [\widehat{\xi} - \xi \mid \bar{\phi}_{k}] \|_{\mathbb{P},2} \cdot \| \Pi^{\perp} [b_{\bm{\beta}} - \widehat{b}_{\bm{\beta}_{\mathsf{init}}} \mid \bar{\phi}_{k}] \|_{\mathbb{P},2} \lesssim \widetilde{r}_{n,\bm{\beta}}, \end{align*} where $\Pi^{\perp}[f \mid \bar{\phi}_{k}]$ denotes the orthocomplement of the projection of $f$ onto the linear span of $\bar{\phi}_{k}$ (see Section (ref)). Then $\widetilde{\bm{\beta}}^{(2)}$ is consistent and satisfies ${\mathbb{E}}[\|\widetilde{\bm{\beta}}^{(2)} - \bm{\beta}^{*}\|] \lesssim \widetilde{r}_{n,\bm{\beta}}$. \end{theorem} The proof of Theorem (ref) is presented in Section (ref) of the Supplementary Material. The convergence rate $\widetilde{r}_{n,\bm{\beta}}$ comprises two components: the first term captures variance, and the second term represents the bias arising from approximating the residuals of the estimated and true nuisance parameters using the basis $\bar{\phi}_{k}$. The next result characterizes the asymptotic distribution of $\widetilde{\bm{\beta}}^{(2)}$. \begin{theorem} Suppose Assumptions (ref) and (ref)--(ref) hold. Further, assume $\widetilde{r}_{n,\bm{\beta}}^{\alpha_{0}}\log n \rightarrow 0$, where $\alpha_{0}$ is the constant from Assumption (ref), and that \begin{equation} \frac{n}{\sqrt{k + n}} \sup_{ \|\bm{\beta} - \bm{\beta}^{\ast}\| \lesssim \widetilde{r}_{n,\bm{\beta}} } \| \Pi^{\perp} [\widehat{\xi} - \xi | \bar{\phi}_{k}] \|_{\mathbb{P},2} \cdot \| \Pi^{\perp} [b_{\bm{\beta}} - \widehat{b}_{\bm{\beta}_{\mathsf{init}}} | \bar{\phi}_{k}] \|_{\mathbb{P},2} \rightarrow 0. \end{equation} Then \begin{align*} \frac{n}{\sqrt{k + n}} (\widetilde{\bm{\beta}}^{(2)} -\bm{\beta}^{*}) \stackrel{d}{\to} N\left(0, \{\nabla_{\bm{\beta}}\psi (\bm{\beta}^{*})\}^{-1} \left\{ V_{1} ( \bm{\beta}^{\ast} ) + \frac{1}{2} V_2 ( \bm{\beta}^{\ast} ) \right\} \{\nabla_{\bm{\beta}}^\top\psi (\bm{\beta}^{*})\}^{-1} \right), \end{align*} where $V_1(\bm{\beta})$ and $V_2(\bm{\beta})$ are defined in (ref)--(ref) in Section (ref) of the Supplementary Material. In particular, the asymptotic variance depends on the relative scaling between $k$ and $n$: $V_2(\bm{\beta}^{\ast})=0$ when $k\ll n$, while $V_1(\bm{\beta}^{\ast})=0$ when $k\gg n$. Consequently, \begin{itemize} • if $k \ll n$, then \begin{align*} \sqrt{n} (\widetilde{\bm{\beta}}^{(2)} -\bm{\beta}^{*}) & \stackrel{d}{\to} N \left( 0, \{\nabla_{\bm{\beta}}\psi (\bm{\beta}^{*})\}^{-1}V_1 (\bm{\beta}^{\ast}) \{\nabla_{\bm{\beta}}^{\top}\psi (\bm{\beta}^{*})\}^{-1} \right); \end{align*} • if $k \gg n$ (recall Assumption (ref) requires $k = o (n^{2})$), then \begin{align*} \frac{n}{\sqrt{k}} (\widetilde{\bm{\beta}}^{(2)} -\bm{\beta}^{*}) & \stackrel{d}{\to} N \left(0, \{\nabla_{\bm{\beta}}\psi (\bm{\beta}^{*})\}^{-1} \frac{1}{2} V_2 (\bm{\beta}^{\ast}) \{\nabla_{\bm{\beta}}^{\top}\psi (\bm{\beta}^{*})\}^{-1} \right); \end{align*} • and if $k / n \rightarrow \tau$, for some $\tau \in (0, \infty)$, \begin{align*} \sqrt{n} (\widetilde{\bm{\beta}}^{(2)} -\bm{\beta}^{*}) \stackrel{d}{\to} N\left(0, (1+\tau) \{\nabla_{\bm{\beta}}\psi (\bm{\beta}^{*})\}^{-1} \left\{ V_{1} ( \bm{\beta}^{\ast} ) + \frac{1}{2} V_2 ( \bm{\beta}^{\ast} ) \right\} \{\nabla_{\bm{\beta}}^{\top}\psi (\bm{\beta}^{*})\}^{-1} \right). \end{align*} \end{itemize} \end{theorem} The complete proof of Theorem (ref) is presented in Supplementary Material Section (ref). We provide a proof sketch as follows. We begin by expanding $\psi (\bm{\beta}^{*})$ around $\widetilde{\bm{\beta}}^{(2)}$: \begin{align*} 0 = \psi (\bm{\beta}^{*}) = \psi (\widetilde{\bm{\beta}}^{(2)}) + \nabla_{\bm{\beta}}\psi (\bm{\beta}^{\dag}) \cdot (\bm{\beta}^{\ast} - \widetilde{\bm{\beta}}^{(2)}), \end{align*} where $\bm{\beta}^{\dag}$ lies between $\bm{\beta}^{*}$ and $\widetilde{\bm{\beta}}^{(2)}$. Next, we decompose $\psi (\widetilde{\bm{\beta}}^{(2)})$ as follows: \begin{align*} \frac{n}{\sqrt{k+n}} \psi (\widetilde{\bm{\beta}}^{(2)}) &= \frac{n}{\sqrt{k+n}} \{\psi (\widetilde{\bm{\beta}}^{(2)}) - \bar{\psi}_{k} (\widetilde{\bm{\beta}}^{(2)})\} - \frac{n}{\sqrt{k+n}} \{\widetilde{\psi}_{k}^{(2)} (\bm{\beta}^{*}) - \bar{\psi}_{k} (\bm{\beta}^{\ast})\}\\ & \quad + \frac{n}{\sqrt{k+n}} \{ (\bar{\psi}_{k} (\widetilde{\bm{\beta}}^{(2)}) - \widetilde{\psi}_{k}^{(2)} (\widetilde{\bm{\beta}}^{(2)})) - (\bar{\psi}_{k} (\bm{\beta}^{*}) - \widetilde{\psi}_{k}^{(2)} (\bm{\beta}^{*}))\} \\ & \eqqcolon A_n + B_n - C_n . \end{align*} Term $A_n$ is bias-related and is $o_{\mathbb{P}}(1)$ under Condition (ref). Term $B_n$ is a centered $U$-statistic, determining the asymptotic distribution of $\widetilde{\bm{\beta}}^{(2)}$. Its asymptotic normality follows from the classical martingale CLT bhattacharya1992class. Since the kernel of this $U$-statistic depends on $k$, the linear component in the Hoeffding decomposition does not always dominate and thus the order of its asymptotic variance depends on the scaling between $n$ and $k$, split into three different regimes as indicated in Theorem (ref). The main technical challenge lies in analyzing the $U$-process term $C_n$, indexed by $\widetilde{\bm{\beta}}^{(2)}$ and defined over a function class changing with the sample size $n$. Existing $U$-process maximal inequalities, such as those of chen2020jackknife, are conservative for our setting. To address this, we establish a refined local maximal inequality in Supplementary Material Section (ref). \begin{remark} Edgeworth expansions offer a different notion of optimality when non-asymptotic or finite sample performance is of concern. cattaneo2025higher study inference for the density-weighted average derivative effect (DWADE) parameter $\beta^{\ast}_0$ (see Remark (ref)) under a randomized experiment, i.e., the propensity score $\mathsf{p}_{T|X}$ is known, using a kernel-based estimator $\widehat{\beta}_0$. Their estimator is a second-order $U$-statistic with Hoeffding decomposition $\widehat{\beta}_0-\beta^{\ast}_0 = \bar{L} + \bar{Q}$, where $\bar{L}$ and $\bar{Q}$ are the linear and quadratic terms, respectively. Classical asymptotics impose bandwidth conditions that make $\bar{Q}$ negligible, yielding an asymptotic linear representation. With smaller bandwidths, $\bar{Q}$ becomes non-negligible, leading to a quadratic approximation where $\mathbb{V}[\widehat{\beta}_0] = \mathbb{V}[\bar{L}] + \mathbb{V}[\bar{Q}]$ and $\mathbb{V}[\widehat{\beta}_0]^{-1/2}(\widehat{\beta}_0-\beta^{\ast}_0) \overset{d}{\to} N(0,1)$. Using Edgeworth expansions, they derive refined distributional approximations that improve finite-sample confidence interval coverage. Our estimator $\widehat{\psi}^{(2)}_k$ of the estimating equation is also a $U$-statistic whose asymptotic distribution depends on quadratic terms when $k \gtrsim n$. However, the goal of cattaneo2025higher is inference refinement (improving finite-sample coverage), while ours is bias correction to achieve rate-optimal estimation under low smoothness conditions. In future work, we can also consider such more refined higher-order distributional approximation if concerned about the finite-sample performance beyond the asymptotic distribution results in Theorem (ref). \end{remark} We next derive the convergence rate of the feasible estimator $\widehat{\bm{\beta}}^{(2)}$. To this end, we need the following assumption on $\widehat{\Sigma}_{k}$. \begin{assumption} We assume that the estimator $\widehat{\Sigma}_{k}$ of $\Sigma_{k}$ satisfies the following condition: \begin{align*} \Vert \widehat{\xi} - \xi \Vert_{\mathbb{P},2} \cdot \left( \Vert \bm{\beta}_{\mathsf{init}} - \bm{\beta}^{\ast} \Vert + \Vert \widehat{b}_{\bm{\beta}_{\mathsf{init}}} - b_{\bm{\beta}_{\mathsf{init}}} \Vert_{\mathbb{P},2} \right) \cdot \Vert \widehat{\Sigma}_{k} - \Sigma_{k} \Vert_{\mathrm{op}} = o_{\mathbb{P}} (n^{- 1 / 2}). \end{align*} \end{assumption} Assumption (ref) ensures that estimating $\Sigma_{k}$ by $\widehat{\Sigma}_{k}$ does not introduce excessive bias. When $\bm{\beta}_{\mathsf{init}}$ is estimated using only $\widehat{\xi}$, we have $\Vert \bm{\beta}_{\mathsf{init}} - \bm{\beta}^{\ast} \Vert \lesssim \Vert \widehat{\xi} - \xi \Vert_{\mathbb{P},2}$. Then the condition in Assumption (ref) reduces to \begin{equation} \left( \Vert \widehat{\xi} - \xi \Vert_{\mathbb{P},2}^{2} + \Vert \widehat{\xi} - \xi \Vert_{\mathbb{P},2} \cdot \Vert \widehat{b}_{ \bm{\beta}_{\mathsf{init}} } - b_{ \bm{\beta}_{\mathsf{init}} } \Vert_{\mathbb{P},2} \right) \cdot \Vert \widehat{\Sigma}_{k} - \Sigma_{k} \Vert_{\mathrm{op}} = o_{\mathbb{P}} (n^{-1 / 2}). \end{equation} Under Assumption (ref), the first two terms satisfy $\Vert \widehat{\xi} - \xi \Vert_{\mathbb{P},2}^{2} \lesssim n^{- \frac{ 2s_{1} }{d + 2 s_{1}}}$ and $\Vert \widehat{\xi} - \xi \Vert_{\mathbb{P},2} \cdot \Vert \widehat{b}_{\bm{\beta}_{\mathsf{init}}} - b_{\bm{\beta}_{\mathsf{init}}} \Vert_{\mathbb{P},2} \lesssim n^{- \frac{s_{1}}{d + 2 s_{1}} - \frac{s_{2}}{d + 2 s_{2}}}$. Consequently, (ref) requires \begin{align*} \Vert \widehat{\Sigma}_{k} - \Sigma_{k} \Vert_{\mathrm{op}} \lesssim n^{- \frac{(d - 2 s_{1} ) \vee 0}{2 (d + 2 s_{1})}} \wedge n^{- \frac{(d^{2} - 4 s_{1} s_{2}) \vee 0}{2 (d + 2 s_{1}) (d + 2 s_{2})}}, \end{align*} which holds if the joint density $\mathsf{p}_{X, T} (\cdot, \cdot)$ is sufficiently smooth. Similar assumptions appear in most recent papers related to HOIF frameworks kennedy2024minimax, bonvini2022fast, mcgrath2024nuisance and can be relaxed using diverging-order $U$-statistic estimators, as detailed in Supplementary Material Section (ref) following liu2017semiparametric or robins2023minimax. Theorem (ref) below establishes the convergence rate of the feasible estimator $\widehat{\bm{\beta}}^{(2)}$ using Theorem (ref) together with Assumption (ref). A proof is provided in Supplementary Material Section (ref). \begin{theorem} Suppose that for some constant $c>0$, $ \Vert {\mathbb{E}} [ \widehat{\psi}^{(2)}_{k}(\bm{\beta}_1) ] - {\mathbb{E}} [ \widehat{\psi}^{(2)}_{k}(\bm{\beta}_2) ] \Vert \geq c \|\bm{\beta}_1 - \bm{\beta}_2\| .$ Under the assumptions of Theorem (ref), we have \begin{align*} {\mathbb{E}}[\Vert \widehat{\bm{\beta}}^{(2)} - \bm{\beta}^{\ast} \Vert] \lesssim \widetilde{r}_{n,\bm{\beta}} + \|\xi - \widehat{\xi}\|_{\mathbb{P},2} \cdot \|b_{\bm{\beta}^*} - \widehat{b}_{\bm{\beta}_{\mathsf{init}} }\|_{\mathbb{P},2} \cdot \|\widehat{\Sigma}_k - \Sigma_k\|_{\mathrm{op}}. \end{align*} If we further impose Assumption (ref), then ${\mathbb{E}}[\Vert \widehat{\bm{\beta}}^{(2)} - \bm{\beta}^{\ast} \Vert] \lesssim \widetilde{r}_{n,\bm{\beta}}$. \end{theorem} The final result in this section characterizes the asymptotic normality of the feasible second-order estimator $\widehat{\bm{\beta}}^{(2)}$, which is a direct consequence of Theorem (ref) and Assumption (ref). \begin{theorem} Under the assumptions of Theorem (ref), together with Assumption (ref), all conclusions in Theorem (ref) continue to hold for $\widehat{\bm{\beta}}^{(2)}$. \end{theorem} \begin{remark} We briefly discuss how the rate of the initial estimator $\bm{\beta}_{\mathsf{init}}$ impacts the results in this section. Condition (ref) in Theorem (ref) specifies the approximation power of $\bar{\phi}_{k}$ for the residuals between estimated and true nuisance parameters. Suppose that $b_{\bm{\beta}}$ is differentiable with respect to $\bm{\beta}$. It is natural to assume that the second-order estimator is at least as close to the truth as the sub-optimal initial estimator, i.e., $\widetilde{r}_{n,\bm{\beta}} \lesssim \|\bm{\beta}_{\mathsf{init}} - \bm{\beta}^{\ast}\|$. Under this assumption, using the triangle inequality and a mean-value expansion in $\bm{\beta}$, we can decompose the left-hand side of condition (ref) as \begin{align} & \ \sup_{\|\bm{\beta} - \bm{\beta}^{\ast}\| \leq \widetilde{r}_{n,\bm{\beta}} } \| \Pi^{\perp} (\widehat{\xi} - \xi \mid \bar{\phi}_{k}) \|_{\mathbb{P},2} \cdot \| \Pi^{\perp} [b_{\bm{\beta}} - \widehat{b}_{\bm{\beta}_{\mathsf{init}}} \mid \bar{\phi}_{k}] \|_{\mathbb{P},2} \\ \leq & \ \| \Pi^{\perp} (\widehat{\xi} - \xi \mid \bar{\phi}_{k}) \|_{\mathbb{P},2} \cdot \Big\{ \sup_{\|\bm{\beta} - \bm{\beta}^{\ast}\| \leq \widetilde{r}_{n,\bm{\beta}} } \| \Pi^{\perp} (b_{\bm{\beta}} - b_{\bm{\beta}_{\mathsf{init}}} \mid \bar{\phi}_{k}) \|_{\mathbb{P},2} + \| \Pi^{\perp} (b_{\bm{\beta}_{\mathsf{init}}} - \widehat{b}_{\bm{\beta}_{\mathsf{init}}} \mid \bar{\phi}_{k}) \|_{\mathbb{P},2} \Big\} \notag \\ \lesssim & \ \| \Pi^{\perp} (\widehat{\xi} - \xi \mid \bar{\phi}_{k}) \|_{\mathbb{P},2} \cdot \Big\{ \sup_{\bm{\beta}^{\dag} \in \mathcal{N}} \| \Pi^{\perp} (\nabla_{\bm{\beta}}b_{\bm{\beta}^{\dag}} \mid \bar{\phi}_{k}) \|_{\mathbb{P},2} \cdot \|\bm{\beta}_{\mathsf{init}} - \bm{\beta}^{\ast}\| + \| \Pi^{\perp} (b_{\bm{\beta}_{\mathsf{init}}} - \widehat{b}_{\bm{\beta}_{\mathsf{init}}} \mid \bar{\phi}_{k}) \|_{\mathbb{P},2}\Big\}. \notag \end{align} Here, $\mathcal{N}\coloneqq \{ \bm{\beta}: \|\bm{\beta}-\bm{\beta}^{*}\| \lesssim \|\bm{\beta}_{\mathsf{init}}-\bm{\beta}^{*}\| \}$, which contains all line segments between $\bm{\beta}_{\mathsf{init}}$ and any $\bm{\beta}$ with $\|\bm{\beta}-\bm{\beta}^{*}\|\le \widetilde r_{n,\bm{\beta}}$. The expression (ref) will be used in Section (ref) to further clarify the rate requirement for $\bm{\beta}_{\mathsf{init}}$ in QTE with binary treatment. \end{remark} \begin{remark} As suggested by a referee, the econometrics literature contains several related proposals for bias correction in two-step semiparametric settings, including jackknife- and bootstrap-based methods cattaneo2019two, cattaneo2018kernel. In Supplementary Material Section (ref), we use the average treatment effect (ATE) as a concrete example to compare our approach with these alternatives. All three methods address the problem that first-stage nuisance estimation can bias the second-stage estimator. However, the sources of bias and the correction mechanisms differ. cattaneo2019two and cattaneo2018kernel aim to estimate the estimating equation $\psi(\bm{\beta})$ itself, with nuisance functions estimated by linear regression and kernel methods, respectively. The resulting estimator is affected by leave-in bias due to the reuse of the same sample in both steps, particularly when the covariate dimension is large or the bandwidth is small. By contrast, we use sample splitting and focus on estimating the approximation bias defined in (ref). Our construction accommodates generic nuisance learners, such as logistic regression, random forests, or neural networks, because we use linear approximation only for the nuisance estimation residuals $\widehat{\xi}-\xi$ and $\widehat{b}_{\bm{\beta}_{\mathsf{init}}}-b_{\bm{\beta}}$, although to obtain sharp theoretical results we need to impose certain smoothness assumptions on $\widehat{\xi} - \xi$ and $\widehat{b}_{\bm{\beta}_{\mathsf{init}}}-b_{\bm{\beta}}$ (see Assumption (ref)). Despite these differences, there are important connections. Specifically, we show that, as cattaneo2019two impose (approximately) linear assumptions on the nuisance functions whereas we do so only for the residuals of the nuisance estimates, both our target parameter and that of cattaneo2019two reduce to a bilinear functional robins2007comment, bruns2026augmented of the form $\mathrm{B}_{\psi,k} (\bm{\beta})$ defined in (ref). We show that the jackknife- and bootstrap-based (under some assumption on the scaling between $k$ and $n$) bias-corrected estimators are, respectively, exactly equivalent and asymptotically equivalent to our second-order $U$-statistic estimator $\widetilde{\mathrm{B}}_{\psi,k}(\bm{\beta})$ of $\mathrm{B}_{\psi,k}(\bm{\beta})$. Moreover, the HOIF framework provides a systematic way to construct higher-order $U$-statistic estimators even when $\Sigma_k$ is unknown (see Supplementary Material Section (ref)), which we believe can also motivate similar higher-order jackknife- or bootstrap-based bias reduction methods koltchinskii2022bootstrap. \end{remark} \begin{remark} As also noted by one referee, both cavaliere2024bootstrap and our paper study settings where a non-negligible asymptotic bias affects the limiting distribution. However, the way how bias is handled differs. cavaliere2024bootstrap consider statistics whose limiting distribution takes the form $T_n = \sqrt{n}(\widehat{\theta}_n - \theta) \stackrel{d}{\to} B_n + Z_1$, where $Z_1$ is centered normal and $B_n$ is a non-vanishing asymptotic bias. They construct a bootstrap analogue $T_n^* := \sqrt{n}(\widehat{\theta}_n^* - \widehat{\theta}_n)$, where $\widehat{\theta}_n^*$ is a bootstrap version of $\widehat{\theta}_n$, such that $T_n^* - \widehat{B}_n \stackrel{d^*}{\to} Z_1$, with $\widehat{B}_n$ a bootstrap-based estimator of $B_n$. When the bootstrap fails to replicate the bias, i.e., $\widehat{B}_n - B_n$ does not converge to zero but converges in distribution to a zero-mean random variable, the standard bootstrap $p$-value $$ \widehat{p}_n := \mathbb{P}^*(T_n^* \le T_n) = \mathbb{P}^*\bigl(T_n^* - \widehat{B}_n \le T_n - \widehat{B}_n\bigr) = \mathbb{P}^*\bigl(T_n^* - \widehat{B}_n \le (T_n - B_n) - (\widehat{B}_n - B_n)\bigr) $$ is no longer asymptotically uniform. Instead, its limiting distribution converges to a distribution $H$ that depends on the joint limiting distribution of $(T_n - \widehat{B}_n, T_n^* - \widehat{B}_n)$ but not on the unknown bias $B_n$ itself. The key insight of cavaliere2024bootstrap is to apply a prepivoting transformation: $H(\widehat{p}_n)$ is asymptotically uniform. Since $H$ is unknown, they estimate it via a double bootstrap to obtain $\widehat{H}_n$, and output the corrected $p$-value $\widetilde{p}_n := \widehat{H}_n(\widehat{p}_n)$. Therefore, their procedure does not directly estimate or remove the bias $B_n$; instead, it bypasses the bias by transforming the bootstrap $p$-value and estimates the distribution function $H$. In contrast, our approach directly estimates the bias itself. As shown in Section (ref), we further decompose the first-order bias $ {\mathbb{E}} [\{\widehat{\xi} (X, T) - \xi (X, T)\} \{b_{\bm{\beta}} (X, T) - \widehat{b}_{\bm{\beta}_{\mathsf{init}}} (X, T)\}]$ into a leading bias component $\mathrm{B}_{\psi,k}(\bm{\beta})$ and a truncation bias $\mathrm{TB}_{\psi,k}(\bm{\beta})$ by projecting the residuals of the estimated nuisance functions onto a sieve basis $\bar{\phi}_k$. The leading bias $\mathrm{B}_{\psi,k}(\bm{\beta})$ is then unbiasedly estimated by a $U$-statistic-based correction term $\widetilde{\mathrm{B}}_{\psi,k}(\bm{\beta})$ when $\Sigma_k$ is known (see Section (ref)). Subtracting this correction term from the first-order estimating equation yields a bias-corrected estimator $\widehat{\bm{\beta}}^{(2)}$ that achieves $\sqrt{n}$-consistency under weaker smoothness conditions than first-order methods (see Theorem (ref) and the discussion in Section (ref)). \end{remark} \subsection{Application: QTE with a binary treatment over \text{H\"{o}lder} smoothness classes} In this section, we revisit Section (ref) and specialize the general results from Theorem (ref) to Theorem (ref) to the QTE setting with binary treatment $T \in \{0,1\}$ (see Example (ref)). To keep the exposition simple, we take the target parameter $\beta_1^*$ to be the $\tau^{th}$ quantile of the potential outcome $Y(1)$, i.e. $\beta_1^* = \inf\{q: \mathbb{P}(Y (1) \leq q) \geq \tau\}$. Under Assumption (ref), $\beta_1^*$ is identified as the unique solution to \begin{align*} \psi(\beta_1) = {\mathbb{E}} (\tau - \mathbbm{1} \{Y(1) \leq \beta_1\}) = {\mathbb{E}} \left\{ \frac{T}{{\mathbb{P}} (T= 1 \mid X )} \cdot (\tau - \mathbbm{1} \{Y \leq \beta_1\}) \right\}. \end{align*} This example fits into our framework with the following correspondences: \begin{align*} & \Gamma(Y(t),t,\beta_1) = \left(\tau - \mathbbm{1} \{Y(t) \leq \beta_1\} \right) \cdot t ,\quad \xi(x,t) = \frac{t}{{\mathbb{P}} (T= 1 \mid X = x )}, \\ & b_{\beta_1} (x,t) = t \cdot {\mathbb{E}} (\tau - \mathbbm{1} \{Y \leq \beta_1 \} \mid X = x, T=t). \end{align*} In the QTE setting, several general regularity conditions for GTMs from Section (ref) (Assumptions (ref)–(ref)) can be simplified. In particular, Assumption (ref) holds automatically because $T$ is binary. Assumption (ref) also holds for QTE under the regularity conditions stated below; see Supplementary Material Section (ref) for verification details. To make this section self-contained, we therefore omit Assumptions (ref) and (ref) and restate the remaining relevant conditions as follows. \begin{assumption} (i) The parameter space ${\mathcal{B}}$ is compact and $\beta_1^*$ lies in its interior; (ii) $\beta^{\ast}_1$ is the unique solution to $\psi (\beta_1)= 0$, and for every $\epsilon > 0$, $\inf_{|\beta_1 - \beta_1^{\ast}| \geq \epsilon} |\psi (\beta_1)| > 0$; and (iii) $\psi (\cdot)$ continuously differentiable on ${\mathcal{B}}$ and its derivative function satisfies $\nabla_{\beta_1} \psi (\beta^{\ast}_1)> 0$. \end{assumption} \begin{assumption} (i) $\sup_{\beta_1 \in {\mathcal{B}}} \|\widehat{b}_{\beta_1} (\cdot, 1)\|_{\infty} \lesssim 1$. (ii) ${\mathbb{P}} ( Y \le \beta_1| X = x, T = 1)$ is Lipschitz continuous in $\beta_1$. (iii) For all $x \in {\mathcal{X}} $, $ \ c_5 \leq \mathsf{p}_{X}(x)\leq c_6$. \end{assumption} \begin{assumption} $\xi (\cdot, 1)$ and $\mathsf{p}_{Y|X,T} (\beta_1|\cdot, T=1)$ are \text{H\"{o}lder} smooth with indices $s_{1}$ and $s_{2}$ respectively, for every $\beta_1 \in {\mathcal{B}}$. Let $s = (s_1 + s_2)/2$ denote the average smoothness. The nuisance estimators $\widehat{\xi}$ and $\widehat{b}_{\beta_1}$ converge to $\xi$ and $b_{\beta_1}$ at the minimax optimal rates ($n^{- \frac{s_{1}}{d + 2 s_{1}}}$ and $n^{- \frac{s_{2}}{d + 2 s_{2}}}$) in $L_{2} (\mathbb{P})$. Moreover, we assume $\widehat{\xi}$ and $\widehat{b}_{\beta_{1,\mathsf{init}}}$ belong to \text{H\"{o}lder} smoothness classes with smoothness indices $s_{1} > 0$ and $s_{2} > 0$, respectively. \end{assumption} \begin{theorem} For the QTE parameter $\beta_1^* = \inf\{q: \mathbb{P}(Y (1) \leq q) \geq \tau\}$, suppose that Assumptions (ref) and (ref)--(ref) hold. Then \begin{itemize} • if $s / d > 1 / 4$ and $k = o (n)$, \begin{align*} \sqrt{n} (\widehat{\beta}^{(2)}_1 - \beta^{\ast}_1) \overset{d}{\rightarrow} N \left( 0, \{\nabla_{\beta_1} \psi (\beta^{\ast}_1)\}^{-1} V_{1} (\beta^{\ast}_1) \{\nabla_{\beta_1}^{\top} \psi (\beta^{\ast}_1)\}^{-1} \right); \end{align*} • if $0 < s / d \leq 1 / 4$ and $k \asymp n^{\frac{2}{1 + 4 s / d}}$, \begin{equation*} {\mathbb{E}}^{1 / 2} [|\widehat{\beta}^{(2)}_1 - \beta^{\ast}_1|^{2}] \lesssim n^{- \frac{4 s / d}{1 + 4 s / d}}. \end{equation*} \end{itemize} \end{theorem} The convergence rates for QTE in Theorem (ref) match those for ATE established in robins2016technical, which shows that these rates are minimax optimal for ATE. The proof of Theorem (ref) is straightforward and proceeds as follows. Under the assumptions of Theorem (ref) with an appropriate choice of basis $\bar{\phi}_{k}$, we have \begin{align*} \| \Pi^{\perp} (\widehat{\xi} - \xi \mid \bar{\phi}_{k}) \|_{\mathbb{P},2} = O(k^{-s_1/d}),\quad \| \Pi^{\perp} (b_{\beta_{1,\mathsf{init}}} - \widehat{b}_{\beta_{1,\mathsf{init}}} \mid \bar{\phi}_{k}) \|_{\mathbb{P},2} = O(k^{-s_2/d}), \end{align*} because the \text{H\"{o}lder}-$s_2$ smoothness of $\mathsf{p}_{Y|X,T}(\beta_1 \mid \cdot, T=1)$ implies that $b_{\beta_1}(\cdot,1)$ is also \text{H\"{o}lder}-$s_2$ smooth for any $\beta_1 \in {\mathcal{B}}$. Furthermore, a direct calculation shows that $\nabla_{\beta_1} b_{\beta_1} (\cdot, 1) = - \mathsf{p}_{Y|X,T}(\beta_1 |\cdot, T = 1 ) $, which is also \text{H\"{o}lder}-$s_{2}$ smooth for any $\beta_1 \in {\mathcal{B}}$. Combining these results with (ref), the bias of $\widehat{\beta}^{(2)}_1$ is bounded above by \begin{align*} k^{- 2 s / d} \cdot | \beta_{1,\mathsf{init}} - \beta^{\ast}_1 | + k^{- 2 s / d} = \left\{ \begin{array}{ll} o (n^{-1 / 2}) & s / d > 1 / 4, k = o (n), \\ O \left( n^{- \frac{4 s / d}{1 + 4 s / d}} \right) & 0 < s / d \leq 1 / 4, k \asymp n^{\frac{2}{1 + 4 s / d}}. \end{array} \right. \end{align*} Therefore, somewhat surprisingly, the higher-order estimator imposes no rate requirement on $\beta_{1,\mathsf{init}}$, as long as $\beta_{1,\mathsf{init}}$ is bounded. This stands in contrast to the first-order LDML estimator, which, as discussed following Proposition (ref), requires $s_1/d > 1/2$. We summarize these comparisons in Table (ref) below. \begin{table}[htbp] \begin{tabular}{l|cc} \toprule & \makecell{\textbf{Smoothness conditions}} & \makecell{\textbf{Rate condition for $\beta_{1,\mathsf{init}}$}} \\ \midrule Our paper & $\displaystyle \frac{s_1/d+s_2/d}{2}>\frac14$ & $\displaystyle \|\beta_{1,\mathsf{init}} - \beta^{*}_1\| = O_{\mathbb{P}}(1)$ \\ kallus2024localized & $\displaystyle {s_{1}}/{d} > \frac{1}{2}$ & $\displaystyle \frac{s_{1}/d}{1 + 2 s_{1}/d} + \frac{s_{2}/d}{1 + 2 s_{2}/d} > \frac{1}{2}$ & $\displaystyle \|\beta_{1,\mathsf{init}} - \beta^{*}_1\| = o_{\mathbb{P}}(n^{-1/4})$ \\ ai2021unified & $\displaystyle {s_1}/{d}>2$ & no $\beta_{1,\mathsf{init}}$ required \\ \bottomrule \end{tabular} \caption{Comparison of smoothness and initial estimator rate conditions for $\sqrt{n}$-consistent estimation of the QTE parameter $\beta_1^* = \inf\{q: \mathbb{P}(Y(1) \leq q) \geq \tau\}$. Here $s_1$ denotes the \text{H\"{o}lder} smoothness index of the propensity score $\xi(\cdot, 1) = 1/\mathbb{P}(T=1 \mid X = \cdot)$, and $s_2$ denotes the Hölder smoothness index of the conditional outcome density $\mathsf{p}_{Y|X,T}(y \mid \cdot, T=t)$ for every $y$.} \end{table} \begin{remark} For ATE in Example (ref), a direct calculation yields $b_{\bm{\beta}}(x,t) = ( ({\mathbb{E}}(Y | X = x, T = 0 ) - \beta_0)\cdot(1-t), ( {\mathbb{E}}(Y |X = x, T = 1 ) - \beta_1)\cdot t )^{\top}$ and $\nabla_{\bm{\beta}}b_{\bm{\beta}} (x,t) = \text{diag}(1-t, t)$. Each component of $\nabla_{\bm{\beta}}b_{\bm{\beta}}$ is infinitely smooth. As a result, condition (ref) imposes no further convergence rate requirements on the initial estimator $\bm{\beta}_{\mathsf{init}}$ of $\bm{\beta}^{\ast}$. \end{remark} \section{Higher-Order Estimators: Toward Infinite-Dimensional Parameters} For continuous or high-dimensional treatments su2019non, fixed-dimensional parameters in (ref) are often inadequate. Motivated by bonvini2022fast and colangelo2026double, we extend our framework to infinite-dimensional GTMs by setting $\omega(\cdot) \equiv \mathsf{p}_T(\cdot)\delta_t(\cdot)$ and letting $p \to \infty$ in (ref). One purpose of this section is to further advocate the fruitful philosophy of turning an infinite-dimensional problem into a finite but diverging-dimensional functional estimation problem kennedy2024minimax. For notational convenience, we denote the infinite-dimensional GTM parameter as $\beta_t^* \equiv \beta^*(t): {\mathcal{T}} \to {\mathbb{R}}$, a function of the treatment value $T=t$. Then $\beta_t^*$ solves \begin{equation*} \int_{u\in{\mathcal{T}}} {\mathbb{E}} \{\Gamma (Y (u) , u, \beta^{\ast}_{t})\} \mathsf{p}_T(u) \delta_{t} (u) {\mathrm d} u \equiv 0 , \ \forall \ t \in {\mathcal{T}}. \end{equation*} Under Assumption (ref), $\beta^{\ast}_{t}$ is identified as the solution to: \begin{equation} \psi_{t} (\beta) \equiv \psi_{t} (\theta_{\beta}) \coloneqq {\mathbb{E}} \left\{ \xi (X, T) \Gamma (Y, T, \beta) \delta_{t} (T) \right\} \equiv 0, \end{equation} where $\xi (x,t) = \mathsf{p}_T(t)/\mathsf{p}_{T|X}(t \mid x)$. Here $\Gamma: {\mathcal{Y}} \times {\mathcal{T}} \times {\mathcal{B}} \rightarrow {\mathbb{R}}$ is a known one-dimensional generalized residual function. To simplify notation, we maintain the conventions of Section (ref): for any $\beta\in{\mathcal{B}}$, we define $b_{\beta}(x,t) = {\mathbb{E}} \{\Gamma (Y, T, \beta) \mid X = x, T=t\}$. It is worth noting that in this section we are interested in estimating the GTM parameter at some fixed treatment value $T = t$. \subsection{Higher-order estimators for infinite-dimensional parameters} Since the Dirac delta function $\delta_{t} (\cdot)$ cannot be evaluated in practice, we approximate it with a kernel function, following bonvini2022fast; see Remark (ref) for details. Let $K (\cdot)$ be a base kernel function and define \begin{align*} K_{t, h} (u) \coloneqq \frac{1}{h} K \left( \frac{u-t}{h} \right) , \quad g(x) \coloneqq \int_{u\in{\mathcal{T}}} K_{t,h} (u) \mathsf{p}_{X,T} ( x, u) {\mathrm d} u. \end{align*} We assume the kernel $K(\cdot)$ satisfies the following condition, which is also imposed by bonvini2022fast. \begin{assumption} $K (\cdot): {\mathcal{T}} \rightarrow {\mathbb{R}} $ is a uniformly bounded kernel function of order $\alpha_1 \wedge \alpha_2$, where $\alpha_1,\alpha_2$ are positive integers specified in Assumption (ref). It satisfies $\int_{{\mathcal{T}}} K (u) {\mathrm d} u = 1$, $\int_{{\mathcal{T}}} u^{i}K (u) {\mathrm d} u = 0$ for $i=1,\ldots, \alpha_1 \wedge \alpha_2-1$, and $\int_{{\mathcal{T}}} u^{\alpha_1 \wedge \alpha_2} K (u) {\mathrm d} u \eqqcolon \kappa_{\alpha_1 \wedge \alpha_2,1} \neq 0$. Moreover, $g(x) \in (\underline{c}_{g}, \bar{c}_g)$ for any $x$ and $h$, where $0 < \underline{c}_g < \bar{c}_g < \infty$ are two universal constants. \end{assumption} Intuitively, $K_{t, h} (\cdot)\to \delta_{t} (\cdot)$ as $h\to 0$, so we replace $\delta_{t}(u)$ in (ref) with $K_{t,h}(u)$. Analogous to (ref), we construct the first-order estimating equation: \begin{align*} \widehat{\psi}_{\beta, t}^{(1)} \coloneqq {\mathbb{U}}_{n, 2} \left[ K_{t,h} (T_1) \widehat{\xi} (X_1, T_1) \{ \Gamma (Y_1, T_1, \beta) - \widehat{b}_{\beta_{\mathsf{init}}} (X_1, T_1) \} + \widehat{b}_{\beta_{\mathsf{init}}} (X_1, T_2)K_{t,h}(T_2) \right]. \end{align*} Following the debiasing argument in Section (ref), we project the localized nuisance errors $\xi(\cdot,t)-\widehat{\xi}(\cdot,t)$ and $b_\beta(\cdot,t)-\widehat{b}_{\beta_{\mathsf{init}}}(\cdot,t)$ onto the dictionary $\bar{{\mathsf{z}}}_{k_x}(x)$ with respect to the weight $g$ (see Supplementary Material Section (ref)). This yields the oracle second-order estimating equation $\widetilde{\psi}_{t,k}^{(2)}(\beta) := \widehat{\psi}_t^{(1)}(\beta) - \widetilde{\mathrm{B}}_{\psi,t,k}(\beta)$, where \begin{align*} &\widetilde{\mathrm{B}}_{\psi, t,k} (\beta)\coloneqq {\mathbb{U}}_{n, 3} \{K_{t,h}(T_{1}) \widehat{\xi} (X_{1}, T_{1}) - K_{t,h} (T_{3} )\} \bar{{\mathsf{z}}}_{k_{x}}^{\top} (X_{1}) \widetilde{\Omega}_{k_x}^{-1}\\ &\qquad\qquad\qquad\qquad\times \bar{{\mathsf{z}}}_{k_{x}} (X_{2}) \{ \Gamma (Y_2, T_2, \beta) - \widehat{b}_{\beta_{\mathsf{init}}} (X_{2}, T_{2})\} K_{t,h} (T_{2}), \end{align*} with $\widetilde{\Omega}_{k_x} \coloneqq {\mathbb{E}} \{K_{t,h}(T) \bar{{\mathsf{z}}}_{k_{x}} (X)\bar{{\mathsf{z}}}_{k_{x}}(X)^{\top}\} = \int_{{\mathcal{X}}} g(x)\bar{{\mathsf{z}}}_{k_{x}} (x)\bar{{\mathsf{z}}}_{k_{x}}(x)^{\top} {\mathrm d} x $ treated as known. Here $\beta_{\mathsf{init}}$ is an initial estimator of $\beta^{\ast}_t$ computed from the nuisance sample. The oracle second-order estimator $\widetilde{\beta}^{(2)}_{t}$ is defined by solving $\widetilde{\psi}_{t,k}^{(2)} (\beta)= 0$. \subsection{Convergence rates of higher-order estimators} To proceed, we impose a set of regularity conditions that closely parallel those for the finite-dimensional case. For brevity, we assume Assumptions (ref) and (ref)–(ref) hold with appropriate notational adjustments. Assumptions (ref) and (ref) are restated below in simplified form. \begin{assumption} For any $t\in{\mathcal{T}}$: (i) the parameter space ${\mathcal{B}}$ is compact and $\beta^{*}_t$ is in the interior of ${\mathcal{B}}$. (ii) $\beta^{\ast}_{t}$ is the unique solution to $\psi_{t} (\beta)= 0$, and $\inf_{|\beta - \beta^{\ast}_{t}| \geq \epsilon} |\psi_{t} (\beta)| > 0$ for every $\epsilon > 0$. (iii) $\psi_{t} (\beta)$ is continuously differentiable in $\beta$ with $\nabla_{\beta} \psi_{t} (\beta^{\ast}_{t})> 0$. \end{assumption} Additionally, we impose the following smoothness conditions on the nuisance functions and their estimators, which facilitate kernel smoothing and are standard in nonparametric estimation. \begin{assumption} The functions $t \mapsto \xi(x, t)$ and $t \mapsto \mathsf{p}_{T}(t)$, $t \mapsto \widehat{\xi}(x, t)$ are $\alpha_{1}$-times continuously differentiable with uniformly bounded derivatives for any $x \in {\mathcal{X}}$. Analogously, $t \mapsto b_{\beta} (x, t)$ and $ t \mapsto \widehat{b}_{\beta} (x, t)$ are $\alpha_{2}$-times continuously differentiable with uniformly bounded derivatives for any $x \in {\mathcal{X}}$ and $\beta\in{\mathcal{B}}$. \end{assumption} Theorem (ref) is the main result of this section, characterizing the convergence rate of $\widetilde{\beta}^{(2)}_t$. To reduce clutter, define $\Delta_{b,t} (X;\beta) \coloneqq b_{\beta} (X,t) - \widehat{b}_{\beta_{\mathsf{init}}} (X,t)$ and $\Delta_{\xi,t} (X) \coloneqq \xi (X,t) - \widehat{\xi} (X,t).$ \begin{theorem} Suppose that Assumptions (ref), (ref)--(ref), and (ref)--(ref) hold. Let $\widetilde{r}_{n,\mathsf{np}} \to 0$ be a diminishing sequence as $n\to \infty$. Suppose that $ h (\log n)^2 \to 0 $, $(\widetilde{r}_{n,\mathsf{np}})^{\alpha_{0}} \log n \to 0$ and \begin{align*} & \left( \frac{\sqrt{k_{x}} }{nh} + \frac{ 1 }{\sqrt{nh}}\right)\log n + h^{\alpha_1 \wedge \alpha_2} + \sup_{ \beta \in{\mathcal{B}} } \left\|\Pi_{g}^{\perp} \{\Delta_{b,t} (\beta) \mid \bar{{\mathsf{z}}}_{k_{x}}\} \right\|_{g,2} \cdot \left\| \Pi_{g}^{\perp} (\Delta_{\xi,t} \mid \bar{{\mathsf{z}}}_{k_{x}}) \right\|_{g,2} \lesssim \widetilde{r}_{n,\mathsf{np}}. \end{align*} Then, $\widetilde{\beta}^{(2)}_{t}$ is a consistent estimator of $\beta^{*}_{t}$ and satisfies \begin{align*} {\mathbb{E}} (|\widetilde{\beta}^{(2)}_{t} - \beta^{*}_{t}|) \lesssim &\ \sup_{|\beta - \beta_t^{\ast}|\leq \widetilde{r}_{n,\mathsf{np}}} \left\|\Pi_{g}^{\perp} \{\Delta_{b,t} (\beta) \mid \bar{{\mathsf{z}}}_{k_{x}}\} \right\|_{g,2} \left\| \Pi_{g}^{\perp} (\Delta_{\xi,t} \mid \bar{{\mathsf{z}}}_{k_{x}}) \right\|_{g,2} + \frac{\sqrt{k_{x}} }{nh} + \frac{ 1 }{\sqrt{nh}} + h^{\alpha_1 \wedge \alpha_2}. \end{align*} Here, $\Pi^{\perp}_{g} (f \mid \bar{{\mathsf{z}}}_{k_x}) (\cdot) \coloneqq f (\cdot) - \int_{{\mathcal{X}}} f(x) \bar{{\mathsf{z}}}_{k_x} (x)^{\top} g(x) {\mathrm d} x \widetilde{\Omega}_{k_x}^{-1} \bar{{\mathsf{z}}}_{k_x} (\cdot)$ denotes the orthocomplement of the $g$-weighted projection onto $\bar{{\mathsf{z}}}_{k_x}$, and $\|f\|_{g,2}^{2} = \int_{{\mathcal{X}}} f(x)^2g(x){\mathrm d} x$. \end{theorem} The proof of Theorem (ref) is given in Supplementary Material Section (ref) and closely follows that of Theorem (ref). Using a similar decomposition, we first establish consistency by showing ${\mathbb{E}} (|\widetilde{\beta}^{(2)}_t - \beta^*_t|) \lesssim \widetilde{r}_{n,\mathsf{np}} \to 0$. Building on this preliminary rate, we refine the analysis to the local region $\{\beta \in {\mathcal{B}}: |\beta - \beta_{t}^{\ast}| \lesssim \widetilde{r}_{n,\mathsf{np}} \}$ to obtain the final convergence rate. The key difference lies in the kernel-weighted $U$-processes introduced by the kernel function. To apply the maximal inequality from Supplementary Material Section (ref), each term in the Hoeffding decomposition must be carefully controlled with bandwidth-dependent envelopes, requiring a refined and delicate argument. These technical calculations distinguish our analysis from that of kennedy2017non, where the target parameter is explicitly defined and does not require such $U$-process techniques. \begin{remark} Following the same argument in Remark (ref), we have \begin{align*} &\sup_{ |\beta - \beta_t^{\ast}|\leq \widetilde{r}_{n,\mathsf{np}} } \|\Pi_{g}^{\perp} \{\Delta_{b,t} (\beta) \mid \bar{{\mathsf{z}}}_{k_{x}}\} \|_{g,2} \| \Pi_{g}^{\perp} (\Delta_{\xi,t} \mid \bar{{\mathsf{z}}}_{k_{x}}) \|_{g,2} \lesssim \| \Pi_{g}^{\perp} \{ \Delta_{b,t} (\beta_{\mathsf{init}}) \mid \bar{{\mathsf{z}}}_{k_{x}}\} \|_{g,2} \| \Pi_{g}^{\perp} (\Delta_{\xi,t} \mid \bar{{\mathsf{z}}}_{k_{x}}) \|_{g,2}\\ &+\ \|\Pi_{g}^{\perp} \{ \nabla_{\beta}b_{\beta^{\dag}} (\cdot,t) \mid \bar{{\mathsf{z}}}_{k_{x}}\} \|_{g,2}\cdot \|\beta_{\mathsf{init}} - \beta^{*}_t\| \cdot \| \Pi_{g}^{\perp} ( \Delta_{\xi,t} \mid \bar{{\mathsf{z}}}_{k_{x}}) \|_{g,2}, \end{align*} where $\beta^\dagger$ lies between $\beta_{\mathsf{init}}$ and $\beta^{\ast}_{t}$. Assume that $\Delta_{\xi,t}(\cdot)$ is $s_{1}$-\text{H\"{o}lder} smooth, while $\Delta_{b,t} (\cdot;\beta_{\mathsf{init}})$ and $\nabla_{\beta}b_{\beta^{\dag}} (\cdot,t)$ are $s_2$-\text{H\"{o}lder} smooth. With an appropriate basis $\bar{{\mathsf{z}}}_{k_{x}}$, we have \begin{align*} &\| \Pi^{\perp} (\Delta_{\xi,t} \mid \bar{{\mathsf{z}}}_{k_{x}}) \|_{g,2} = O(k_{x}^{-s_1/d}),\quad \left\|\Pi_{g}^{\perp} \{\nabla_{\beta}b_{\beta^{\dag}} (\cdot,t) \mid \bar{{\mathsf{z}}}_{k_{x}}\} \right\|_{g,2} = O(k_{x}^{-s_2/d}),\\ &\text{and}\ \| \Pi^{\perp} \{ \Delta_{b,t} (\beta_{\mathsf{init}})\mid \bar{{\mathsf{z}}}_{k_{x}}\} \|_{g,2} = O(k_{x}^{-s_2/d}). \end{align*} Thus, if $|\beta_{\mathsf{init}} - \beta^{*}_t| = O_{\mathbb{P}}(1)$, Theorem (ref) gives \begin{align*} {\mathbb{E}} (| \widetilde{\beta}^{(2)}_{t} - \beta^{*}_{t} |) \lesssim &\ k_{x}^{-\frac{s_1 + s_2}{d}} + \ \frac{\sqrt{k_{x}} }{nh} + \frac{ 1 }{\sqrt{nh}} + h^{\alpha_1 \wedge \alpha_2}. \end{align*} Choosing $k_x \asymp nh$ makes the variance of order $(nh)^{-1/2}$, reducing the rate to $(nh)^{-(s_1+s_2)/d} + (nh)^{-1/2} + h^{\alpha_1 \wedge \alpha_2}$. If the average smoothness satisfies $(s_1+s_2)/(2d) \geq 1/4$ and $h \asymp n^{-1/(2(\alpha_1 \wedge \alpha_2)+1)}$, we obtain the rate $n^{-(\alpha_1 \wedge \alpha_2)/(2(\alpha_1 \wedge \alpha_2)+1)}$. \end{remark} \begin{remark} Our higher-order estimator $\widetilde{\beta}_t^{(2)}$ is related to two existing proposals for average dose–response function (ADRF) estimation. First, colangelo2026double develop a kernel-based DML estimator motivated by an approximate first-order influence function obtained through kernel localization around a target treatment value. Their estimator takes the form \begin{align*} \widehat{\beta}_{\mathrm{CL}} = & \, {\mathbb{U}}_{n, 1}\left\{\frac{K_{t,h}(T_1)\{Y_1-\widehat{\mu}(X_1, t)\}}{ \widehat{\mathsf{p}}_{T \mid X} (T = t \mid X_1 ) }+\widehat{\mu}(X_1, t)\right\}, \end{align*} where $\widehat{\mu}(x,t) \coloneqq {\mathbb{E}} [Y \mid X=x, T=t]$. They establish asymptotic normality of $\widehat{\beta}_{\mathrm{CL}}$ under high-level regularity conditions, including the key product-rate requirement $\sqrt{n h}\|\widehat{\mathsf{p}}_{T \mid X} (T = t \mid \cdot )-\mathsf{p}_{T \mid X} (T = t \mid \cdot )\|_{\mathbb{P},2} \cdot \|\widehat{\mu}(\cdot,t) - \widehat{\mu}(\cdot,t)\|_{\mathbb{P},2}\to 0$. Second, bonvini2022fast propose a second-order estimator (the BK estimator) based on the truncated parameter approach of HOIFs robins2008higher: \begin{align*} \widetilde{\beta}_{\mathrm{BK}} = & \, {\mathbb{U}}_{n, 1}\left\{\frac{K_{t,h}(T_1)\{Y_1-\widehat{\mu}( X_1, t)\}}{ \widehat{\mathsf{p}}_{T \mid X} (T = t \mid X_1 ) }+\widehat{\mu}(t, X_1)\right\} \\ & + \, {\mathbb{U}}_{n, 2}\left\{K_{t,h}(T_{1})\{Y_{1}-\widehat{\mu}(X_{1}, T_{1})\}\bar{{\mathsf{z}}}_{k_{x}}(X_{1})^{\top}\widetilde{\Omega}_{k_x}^{-1}\bar{{\mathsf{z}}}_{k_{x}}(X_{2})\left(\frac{K_{t,h}(T_{2})}{ \widehat{ \mathsf{p}}_{T \mid X} (T_{2} \mid X_{2} )}-1\right)\right\}. \end{align*} In essence, $\widetilde{\beta}_{\mathrm{BK}}$ can be viewed as a higher-order extension of $\widehat{\beta}_{\mathrm{CL}}$ for ADRF. Our oracle second-order estimator is obtained by solving \begin{align*} & 0 = {\mathbb{U}}_{n, 2} \left[ K_{t,h} (T_1) \widehat{\xi} (X_1, T_1) \{ Y_1 - \widehat{\mu} (X_1,T_1) - \beta + \beta_{\mathsf{init}} \} + \{ \widehat{\mu} (X_1, T_2) - \beta_{\mathsf{init}}\}K_{t,h} (T_2) \right]\\ & + {\mathbb{U}}_{n, 3}\big[ \{ K_{t,h}(T_{1})\widehat{\xi} (X_{1}, T_{1}) - K_{t,h} (T_{3} )\}\bar{{\mathsf{z}}}_{k_{x}} (X_{1})^{\top} \widetilde{\Omega}_{k_x}^{-1} \bar{{\mathsf{z}}}_{k_{x}} (X_{2}) \{Y_{2} - \widehat{\mu} (X_{2}, T_{2}) -\beta +\beta_{\mathsf{init}}\}K_{t,h} (T_{2})\big] . \end{align*} Notably, because $b_\beta(X,t) - \widehat{b}_{\beta_{\mathsf{init}}}(X,t) = \mu(X,t) - \beta - \widehat{\mu}(X,t) + \beta_{\mathsf{init}}$ and $\Pi_g^{\perp}[\beta - \beta_{\mathsf{init}} \mid \bar{{\mathsf{z}}}_{k_x}] \equiv 0$ for any $\beta$, the convergence rate in Theorem (ref) simplifies to \begin{align*} \frac{\sqrt{k_{x}} }{nh} + \frac{ 1 }{\sqrt{nh}} + h^{\alpha_1 \wedge \alpha_2} + \left\|\Pi_{g}^{\perp} \left[\mu (\cdot,t) - \widehat{\mu} (\cdot,t) \mid \bar{{\mathsf{z}}}_{k_{x}}\right] \right\|_{g,2} \cdot \left\| \Pi_{g}^{\perp} \left[ \Delta_{\xi,t}\mid \bar{{\mathsf{z}}}_{k_{x}}\right] \right\|_{g,2} \notag, \end{align*} This rate coincides with that of Theorem 1 in bonvini2022fast, confirming that our framework generalizes their result to a broader class of nonparametric GTM parameters, including the QDRF (Example (ref)) as a special case. \end{remark} \begin{remark} When the unknown matrix $\widetilde{\Omega}_{k_x}$ in $\widetilde{\psi}_{t,k}^{(2)}$ is replaced by an estimator $\widehat{\Omega}_{k_x}$ computed from the nuisance sample, we obtain the feasible estimator $\widehat{\psi}_{t,k}^{(2)}$. Correspondingly, we define the feasible estimator $\widehat{\beta}_t^{(2)}$ for $\beta_t^*$ as the solution to $\widehat{\psi}_{t,k}^{(2)}(\beta) = 0$. Following similar arguments as in Theorem (ref), the convergence rate of $\widehat{\beta}_t^{(2)}$ satisfies \begin{equation*} {\mathbb{E}} [| \widehat{\beta}^{(2)}_{t} - \beta^{*}_{t}|] \lesssim\ {\mathbb{E}} [| \widetilde{\beta}^{(2)}_{t} - \beta^{*}_{t}|] + \|\Delta_{\xi,t} \|_{g,2}\cdot \|\Delta_{b,t}(\beta^*_t) \|_{g,2} \cdot \|\widehat{\Omega}_{k_x} - \widetilde{\Omega}_{k_x}\|_{\mathrm{op}}. \end{equation*} This bound mirrors the structure established in Theorem (ref) for the finite-dimensional case, and can serve as a building block for establishing asymptotic normality of the feasible estimator under appropriate conditions on $\widehat{\Omega}_{k_x}$. \end{remark} \section{Numerical Experiments} \subsection{Simulation studies} In this section, we conduct a proof-of-concept simulation study comparing the performance of our proposed higher-order estimator with two competing methods that do not utilize the HOIF framework. We focus on estimating the QTE parameter with binary treatment, $\bm{\beta}^* = (\beta_{0}^*, \beta_{1}^*)^{\top}$, as described in Example (ref). Recall that $\beta_{1}^*$ (resp. $\beta_{0}^*$) denotes the $\tau$-quantile of the treatment (resp. control) group. We choose $\tau = 25\%$ in this simulation. Following xu2022deepmed and liu2017semiparametric, we simulate data from the following two data-generating mechanisms with different numbers of baseline covariates ($d = 1$ for Case 1 and $d = 4$ for Case 2): \begin{itemize} • Case 1: \begin{align*} \left\{ \begin{array}{l} X \sim \mathrm{Uniform} ([-1, 1]), \\ T \mid X \sim \mathrm{Bernoulli} (\mathrm{expit} (\eta (0.5 \cdot X;s))), \\ Y (t) \mid X = (1 + 0.3 \cdot \eta (0.5 \cdot X;s)) \cdot t + 0.2 \cdot \varepsilon, \end{array} \right. \end{align*} • Case 2: \begin{align*} \left\{ \begin{array}{l} X = (X_1,\ldots,X_4)^{\top} \sim \mathrm{Uniform} ([-1, 1]^4), \\ T \mid X \sim \mathrm{Bernoulli} (\mathrm{expit} (\eta (0.5 \cdot X_1;s))), \\ Y (t) \mid X = (1 + 0.3 \cdot \eta (0.5 \cdot X_1;s)) \cdot t + 0.2 \cdot \varepsilon, \end{array} \right. \end{align*} \end{itemize} where $\varepsilon \sim N (0, 1)$, $\mathrm{expit} (x) = 1 / (1 + \exp(-x))$, $\eta(x; s) = \sum_{j \in J, l \in {\mathbb{Z}}} 2 ^{-j (s + 0.25)} w_{j, l}(x)$ with $J = \{0, 3, 6, 9, 10, 16\}$ and $w_{j,l}(\cdot)$ is the D6 father wavelet functions. By construction, $\eta (\cdot; s)$ lies close to the boundary of \text{H\"{o}lder}-smooth functions with smoothness index $s$. We note that this is only an approximation, as exact infinite basis expansions are infeasible in simulations. Cases 1 and 2 share the same true outcome regression and propensity score. The only difference is that Case 2 includes three extra non-contributing covariates. In the analysis, we do not assume such knowledge so the basis $\bar{{\mathsf{z}}}_{k}$ depends on all four covariates by stacking the univariate spline basis blocks together with an intercept term $\bar{{\mathsf{z}}}_{k} (x) = (1,\bar{{\mathsf{z}}}_{k_{1}} (x_{1})^{\top}, \bar{{\mathsf{z}}}_{k_{2}} (x_{2})^{\top}, \bar{{\mathsf{z}}}_{k_{3}} (x_{3})^{\top}, \bar{{\mathsf{z}}}_{k_{4}} (x_{4})^{\top})^{\top}$, which yields a basis vector of dimension $k = \sum_{j = 1}^{4} k_{j}+1$. We vary $s \in \{0.25, 0.4, 0.6\}$ to examine its impact on estimator performance and to validate the theoretical results derived earlier. For each combination of sample size $n$ and smoothness level $s$, we run 1000 replications. In each replication, $2n$ observations are drawn and split equally into a nuisance sample and an estimation sample. We estimate $\bm{\beta}^*$ using three methods: \begin{description} • Our higher-order estimator is defined in Section (ref). For each continuous covariate, we use a degree-1 B-spline expansion with the number of interior knots set to $\lceil n/100\rceil$. As mentioned, we construct the basis $\bar{{\mathsf{z}}}_k(x)$ by stacking the univariate spline basis blocks for all continuous covariates, together with an intercept term. Specifically, in Case 2, $\bar{{\mathsf{z}}}_k(x) = (1, \bar{{\mathsf{z}}}_{k_1}(x_1)^\top, \bar{{\mathsf{z}}}_{k_2}(x_2)^\top, \bar{{\mathsf{z}}}_{k_3}(x_3)^\top, \bar{{\mathsf{z}}}_{k_4}(x_4)^\top)^\top,$ with $k_j=\lceil n/100\rceil+1$ and total dimension $k=\sum_{j=1}^4 k_j + 1$. On the nuisance sample, we estimate the stabilized weight function $\widehat\xi$ via logistic regression, obtain the initial estimator $\bm{\beta}_{\mathrm{init}}$ by stabilized weighting using $\widehat\xi$, and estimate $\widehat b_{\bm{\beta}_{\mathrm{init}}}$ also via logistic regression. • The generalized optimization estimator of ai2021unified, applied to both the nuisance and estimation samples. See Supplementary Material Section (ref) for its explicit form. • The localized debiased machine learning estimator of kallus2024localized, implemented with the same settings as in Section 6.1 of that paper, where all nuisance functions are estimated via logistic regression. \end{description} Both Figure (ref) (for Case 1) and Figure (ref) (for Case 2), displaying simulation results for the QTE $\beta_1^\ast - \beta_0^\ast$, show that when the smoothness parameter $s$ is low (e.g., $s = 0.25$ or more precisely $s$ is close to $0.25$), our higher-order estimator (HOE) substantially outperforms both GOE and LDML, achieving lower MSE and absolute bias, especially when the sample size $n$ is relatively large. In particular, the higher-order estimator is empirically $\sqrt{n}$-consistent even when $s$ is near 0.25 (the lower right panels of Figure (ref) and Figure (ref)). This advantage stems from the difficulty of estimating nuisance functions under low smoothness, which induces large first-stage bias—due to errors in both the nuisance parameter estimates and the initial estimator $\bm{\beta}_{\mathsf{init}}$. Consequently, GOE and LDML fail to achieve $\sqrt{n}$-consistency (see the lower‑right panel). As $s$ increases, all three methods exhibit improvements, with bias, variance, and MSE all decreasing. These empirical findings are consistent with our theoretical results, which indicate that the higher-order estimator enjoys improved statistical properties when the \text{H\"{o}lder} smoothness of the nuisance parameters is low. Figures (ref) and (ref) (for Case 1) and Figures (ref) and (ref) (for Case 2) in Supplementary Material Section (ref) present normal QQ-plots of our higher-order estimators for $\beta_0^*$ and $\beta_1^*$ across all $(n,s)$ combinations. These plots support the theoretical conclusion (Theorem (ref)) that the higher-order estimators are asymptotically normal under appropriate scaling. Finally, we perform sensitivity analyses to examine whether our higher-order estimator is sensitive to the choice of $k$ or to the initial estimator $\bm{\beta}_{\mathsf{init}}$. Recall that, for each continuous covariate, we use a degree-1 B-spline basis with $\lceil n/100\rceil$ interior knots. To examine the sensitivity to the choice of $k$, we multiply $ n/100$ by a factor $\lambda \in \{0.75,1,1.25\}$, so that the number of knots varies over $\lceil \lambda n/100\rceil$, resulting in $k = d\cdot(\lceil \lambda n/100\rceil + 1) + 1$. Figure (ref) in Supplementary Material Section (ref) reports the empirical bias, variance, and mean squared error for HOE QTE estimator across different values of $\lambda$ for both Cases 1 and 2 with $n=2000$ and $s=0.25$, based on 1000 Monte Carlo replications. The figure shows the expected bias--variance tradeoff: as $k$ increases, the empirical bias of the higher-order estimator decreases, whereas the empirical variance increases. In practice, one may choose $k$ such that the point estimates do not differ as much compared to the increase in variance. Here, for both Cases 1 and 2, $\lceil n / 100 \rceil$ is a reasonable choice. A more thorough study of how to choose $k$ is left to a future work that is dedicated to the problem of constructing data-driven basis. To assess sensitivity to the initial estimator $\bm{\beta}_{\mathsf{init}}$, we perturb it deterministically as $ \bm{\beta}_{\mathsf{init}}^{(\delta)} = \bm{\beta}_{\mathsf{init}} + \delta \cdot (-1,1)^\top, $ with $ \delta \in \{-0.5,-0.25,0,0.25,0.5\}$. Here, the perturbation level $\delta$ controls both the magnitude and the direction of the perturbation. Figure (ref) in Supplementary Material Section (ref) reports boxplots of the estimation error of the HOE QTE estimator across different values of $\delta$ for both Cases 1 and 2, in the setting $n=2000$ and $s=0.25$, based on 1000 Monte Carlo replications. The results reported in Figure (ref) in Supplementary Material Section (ref) suggest that HOE is fairly robust to moderate perturbations in the initial estimator. However, when the perturbation becomes too large, the estimator becomes less stable, as reflected by the increased dispersion and a slightly greater number of outliers. \begin{figure}[htbp] \caption{Results for the simulation study Case 1 of estimating the QTE, $\beta_1^*-\beta_0^*$, where $\beta_{1}^*$ (resp. $\beta_{0}^*$) denotes the $25\%$-quantile of the treatment (resp. control) group, using different methods, based on 1000 Monte Carlo replications. } \end{figure} \begin{figure}[htbp] \caption{Results for the simulation study Case 2 of estimating QTE, $\beta_1^*-\beta_0^*$, where $\beta_{1}^*$ (resp. $\beta_{0}^*$) denotes the $25\%$-quantile of the treatment (resp. control) group, using different methods, based on 1000 Monte Carlo replications. } \end{figure} \subsection{Real data analysis} In this section, we apply our higher-order estimators to estimate the QTE of 401(k) plan eligibility on household wealth, using the Survey of Income and Program Participation (SIPP) data studied in benjamin2003does,abadie2003semiparametric, chernozhukov2004effects,kallus2024localized. The sample contains 9{,}915 observations. The treatment $T$ indicates eligibility for a 401(k) plan, and the outcome $Y$ is net financial assets; our goal is to assess whether 401(k) eligibility increases household saving. The dataset consists of a set of pre-treatment covariates $X$, including age, income, family size, education, marital status, two-earner status, DB pension status, IRA participation, and homeownership (see chernozhukov2004effects, Section 3, for details). Recently, kallus2024localized analyzed this dataset using LDML estimators. We implement our higher-order debiased estimators on the same dataset to provide a proof-of-concept for their application in real-world settings. As in our simulation studies, we split the sample into two folds. With the nuisance sample, we compute the initial estimator $\widehat{\bm{\beta}}_{\mathsf{init}}$ by stabilized weighting, and estimate both the stabilized weights $\widehat\xi$ and the generalized outcome regression model $\widehat{b}_{\bm{\beta}_{\mathsf{init}}}$ using four learners: random forest, boosting, LASSO, and a one-hidden-layer neural network. Random forests are implemented using the R package \texttt{randomForest}, boosting using \texttt{gbm}, LASSO using \texttt{hdm}, and the one-hidden-layer neural network using \texttt{nnet}. For LASSO, we include polynomial terms of degree 6 for age, 8 for income, 4 for education, and 2 for family size, together with all binary covariates, and all pairwise interactions among these terms, resulting in a total of 275 predictors. In the debiasing step, we use degree-2 B-spline expansions for the continuous covariates income, age, family size, and education, with interior knot numbers $(\lceil n/200\rceil,\lceil n/200\rceil,4,2)=(25,25,4,2)$, respectively. We then construct $\bar{{\mathsf{z}}}_k(x)$ by stacking the resulting univariate spline basis blocks with the binary covariates and an intercept term. This choice accounts for the heterogeneous supports of the continuous covariates: income and age vary over wide ranges, while family size and education take values on smaller grids. It also ensures numerical stability of $\widehat{\Sigma}_k$ by avoiding very small eigenvalues. For comparison, we also report LDML estimates under the same settings as Section 6.2 of kallus2024localized (with $K=5$), using a fixed two-fold split without re-randomizing folds to maintain transparency. \begin{table}[ht] \begin{tabular}{clcccc} \toprule $\tau$ & & Forest & Neural Net & Boosting & LASSO \\ \midrule \multirow{2}{*}{0.25} & LDML & 1.06 & 1.02 & 1.04 & 1.06 \\ & HOE & 1.01 & 1.09 & 1.09 & 1.11 \\ \midrule \multirow{2}{*}{0.50} & LDML & 4.95 & 4.46 & 4.45 & 4.75 \\ & HOE & 4.74 & 4.37 & 4.40 & 4.35 \\ \midrule \multirow{2}{*}{0.75}& LDML & 13.60 & 11.90 & 12.16 & 12.65 \\ & HOE & 12.68 & 11.77 & 12.15 & 11.90 \\ \bottomrule \end{tabular} \caption{The QTE of 401(k) eligibility in thousand dollars estimated by LDML and higher-order estimates (HOE) using different regression methods (LASSO, neural network, boosting, and random forests). Here $\tau \in \{25\%,50\%,75\%\}$ denotes the quantile level of QTE.} \end{table} Table (ref) reports the point estimates of 25%, 50%, and 75% QTEs. Overall, QTE estimates are generally stable among different nuisance learners, suggesting that the bias of the LDML estimator may be relatively small in this application. The higher-order estimates are close to their LDML counterparts, further strengthening this conclusion. Consequently, the difference between the higher-order estimator and the LDML estimator provides empirical evidence on the magnitude of the potential bias of the LDML estimator liu2024assumption. We examine sensitivity to both the choice of $k$ and the initial estimator $\bm{\beta}_{\mathsf{init}}$. To assess sensitivity to $k$, we increase the number of interior knots for age by multiplying $n/200$ by a factor $\lambda \in \{0.75,1,1.25\}$ while keeping other knot numbers fixed. The higher-order estimators perform reasonably well in this range, as shown in Table (ref) in Supplementary Material Section (ref). Stability is ensured by the numerical inversion of $\widehat{\Sigma}_k$ when $k$ is relatively small compared to $n$, whereas the performance degrades if $\widehat{\Sigma}_k$ approaches singularity. Therefore, a heuristic strategy to explore in the future is to increase $k$ (via knots or degree) until numerical instability occurs. Developing such data-driven rules for selecting $k$ and constructing $\bar{{\mathsf{z}}}_k$ remains an important direction. Sensitivity to $\bm{\beta}_{\mathsf{init}}$ is examined similarly to the simulation studies. We find that the higher-order estimator is generally robust to small perturbations of $\widehat{\bm{\beta}}_{\mathsf{init}}$, as shown in Table (ref) in Supplementary Material Section (ref). \section{Concluding Remarks} In this paper, we have constructed higher-order estimators for implicitly defined parameters frequently encountered in econometrics, including average treatment effects and quantile treatment effects as special cases. By applying the HOIFs framework to implicitly defined parameters, we have contributed to the existing econometrics and causal inference literature. Our proposed estimators achieve improved convergence rates compared to first-order methods, which requires only that the initial estimator be consistent rather than imposing stringent rate conditions. For binary-treatment QTE, our estimator attains $\sqrt{n}$-consistency under the Hölder smoothness condition $s/d > 1/4$, substantially relaxing the requirements of existing approaches and matching the minimal Hölder smoothness condition for the ATE to be $\sqrt{n}$-estimable. We have also extended the framework to infinite-dimensional GTMs, providing a unified approach for estimating nonparametric causal functions such as dose–response curves. Several directions remain for future research. First, while we focused on second-order estimators assuming that $\Sigma_k$ can be estimated sufficiently accurately, relaxing this assumption would require further advances in the theory of higher-order $U$-processes. Second, extending our framework to other complexity-reducing assumptions, such as sparsity classes athey2018approximate, remains an open question. Third, with slight modifications, our framework is also expected to be able to address endogeneity problems (including nonparametric and many weak IV models) breunig2024adaptive. Finally, it would be interesting to compare our higher-order debiased estimator with other related bias-correction or bias-aware methods armstrong2020bias, zheng2025perturbed, bonhomme2026higher. \putbib[Master.bib]