Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
114,944 characters · 20 sections · 30 citation commands
HIGHER-ORDER NEYMAN ORTHOGONALITY IN MOMENT-CONDITION MODELS
\def\spacingset#1 \spacingset{1}
{
}
{\bf Keywords:} Neyman-orthogonality, higher-order bias correction, moment conditions, GMM, U-statistics.
\onehalfspacing
\allowdisplaybreaks \setcounter{equation}{0}
Estimation and inference in the presence of nuisance parameters is a central problem in econometrics and statistics. The generalized method of moments Hansen1982 provides a unifying framework for a vast range of estimators, from ordinary least squares and instrumental-variable estimation to maximum likelihood and causal-inference procedures. In all these settings, the parameter of interest must be estimated in the presence of additional unknown quantities. The accuracy with which these nuisance parameters are estimated can have a substantial effect on inference for the parameter of interest.
In such a setting, it is helpful to work with estimating equations that are orthogonal to the nuisance parameters in the sense of Neyman1959. (First-order) orthogonality means that the first derivative of the expected estimating equation with respect to the nuisance parameter vanishes at the true parameter values.\footnote{Neyman-orthogonal estimating equations play a central role in semiparametric estimation BickelKlaassenRitovWellner1993,Newey1994. They are also at the heart of the recent literature on debiased inference in high-dimensional models BelloniChernozhukovHansen2014,JavanmardMontanari2014,vandeGeerBuhlmannRitovDezeure2014,ZhangZhang2014 and on double/debiased machine learning ChernozhukovChetverikovDemirerDufloHansenNeweyRobins2018.} First-order Neyman orthogonality, when combined with sample splitting, permits construction of $\sqrt{n}$-consistent estimators as long as the nuisance parameter is estimated at a rate faster than $n^{-\nicefrac{1}{4}}$, where $n$ is the sample size ChernozhukovChetverikovDemirerDufloHansenNeweyRobins2018.
The faster-than-$n^{-\nicefrac{1}{4}}$ requirement is binding in many important problems. It fails when the nuisance parameter is high-dimensional relative to the sample size, when its convergence rate is slow, or both. One example is panel data with individual fixed effects, where the incidental parameter problem persists under first-order orthogonality NeymanScott1948,LiLindsayWaterman2003,HahnNewey2004,JochmansWeidner2019,KlineSaggioSoelvsten2020. Another example is high-dimensional regression (e.g., MikushevaSolvsten2025), where approximate sparsity conditions are needed for post-double-Lasso inference to be valid WuthrichZhu2023,SurCandes2019. Two further examples are instrumental-variable estimation, where machine-learning-based variable selection introduces biases that first-order corrections can fail to offset AngristFrandsen2022, and nonparametric models with moderate-dimensional nuisance functions, where minimax rates for estimating those functions fall short of $n^{-\nicefrac{1}{4}}$ RobinsLiTchetgenTchetgenvanderVaart2008,vanderVaart2014.
A natural response is to go beyond first-order orthogonality. MackeySyrgkanisZadik2018 formalized the concept of $q$-th order Neyman orthogonality. The condition requires that all derivatives of the expected estimating equation with respect to the nuisance parameter, up to order $q$, vanish at the true parameter values. They showed that $q$-th order orthogonality, again combined with sample splitting, permits valid inference when the nuisance parameter converges at rate $n^{-\nicefrac{1}{2(q+1)}}$; for example, the rate requirement drops from $n^{-\nicefrac{1}{4}}$ to $n^{-\nicefrac{1}{6}}$ in going from first to second order. Constructing moment functions that are orthogonal beyond first order has nevertheless proved difficult. This paper provides general constructions of $q$-th order Neyman-orthogonal moment functions for parameters that satisfy parametric unconditional moment conditions, with a separate target functional. We derive a closed-form, non-recursive expression for the $q$-th order orthogonal moment function in a general parametric setting where both the target moment function and the one that identifies the nuisance parameter can be nonlinear.
The orthogonal moment function depends on multiple independent copies of the data, and it features a generalized inverse of a Jacobian matrix as an additional nuisance parameter whose dimension does not grow with the orthogonalization order $q$. We first provide explicit moment functions in the case where the moment function for the (original) nuisance parameter is affine; this covers linear regression, instrumental-variable estimation with linear moment conditions, and grouped data models with unit-specific intercepts or slopes. We then extend our construction to the case where both moment functions are nonlinear, in which case additional correction terms are needed. Our explicit expression for orthogonal moment functions in the nonlinear case is based on rooted trees. A noteworthy feature of our approach is that the dimension of the additional nuisance parameter can be dramatically reduced: by replacing the moment function with a transformed version constructed from subdeterminants of the Jacobian matrix, the nuisance can be reduced to a single scalar, even under overidentification.
Our paper is related to several recent contributions. MackeySyrgkanisZadik2018 provide an explicit second-order orthogonal moment for the partially linear regression model, but only when first-stage errors are non-Gaussian; they prove an impossibility result under Gaussianity. Their constructive results do not extend to general models. Our constructions circumvent this impossibility because we rely on independent copies of the data and on ordinary finite-dimensional derivatives, rather than on functional derivatives. An approach to constructing $q$-th order orthogonal estimating equations for conditional-likelihood models was provided in BonhommeJochmansWeidner2025. In recent independent work, ChetverikovSorensenTsyvinski2026 propose a recursive formula for Z-estimation problems and apply it to a “triple-Lasso” estimator with second-order orthogonality. Like our construction, theirs uses independent copies of the data. They also provide numerical evidence that higher-order orthogonality can substantially improve coverage of confidence intervals. Their construction differs from ours in several respects. In particular, they introduce new nuisance parameters at each recursion step, so the nuisance dimension grows with the order of orthogonality. They also restrict attention to exactly-identified moment conditions. In a different direction, RobinsLiTchetgenTchetgenvanderVaart2008, vanderVaart2014, and RobinsLiMukherjeeTchetgenvanderVaart2017 develop higher-order influence functions and U-statistic-based estimators for nonparametric models. Those results rely on nonparametric von Mises expansions that hold only approximately. They also do not provide explicit constructions for parameters defined by finite-dimensional moment conditions.
The rest of the paper is organized as follows. Section (ref) describes the setup. Section (ref) illustrates the main ideas in a heterogeneous-coefficient model for grouped data. Section (ref) presents the general construction for the affine case and the nonlinear case, in turn, and also gives a nuisance-dimension reduction technique. Section (ref) reports implementation and simulation evidence for the heterogeneous coefficients example. Section (ref) establishes the asymptotic distribution of the resulting estimator. Finally, the Appendix contains proofs and additional formulas.
We observe a sample $\{W_i : i = 1, \ldots, n\}$ from an unknown distribution $P_0$. Consider a model defined by two moment conditions,
where $\theta_0 \in \mathbb{R}^{d_\theta}$ is the parameter of interest and $\eta_0 \in \mathbb{R}^{d_\eta}$ is a finite-dimensional nuisance parameter. The functions $g(W_i, \theta, \eta) \in \mathbb{R}^{d_g}$ and $m(W_i, \theta, \eta) \in \mathbb{R}$ are known. All expectations are taken with respect to $P_0$.
The first moment condition identifies $\eta_0$ for given $\theta_0$. We assume $d_g \geq d_\eta$, which allows the nuisance parameter to be either exactly identified or overidentified. The second moment condition is used to identify $\theta_0$. We take $m$ to be scalar throughout, purely for notational convenience. In practice one will often have $m(W_i, \theta, \eta) \in \mathbb{R}^{d_m}$ with $d_m > 1$, in which case our orthogonalization construction is applied to each component of $m$ separately. Identification of $\theta_0$ in a given application requires $d_m \geq d_\theta$, and the target-moment system in $\theta$ may itself be exactly identified or overidentified. These identification considerations do not affect the orthogonalization construction that we present.
Our goal is to construct a Neyman-orthogonal moment function $\psi$ that replaces the target moment function $m$. The orthogonal moment function has expected derivatives in $\eta$ that vanish up to order $q$ at the true parameter values, which makes it less sensitive to estimation error in $\eta$. To achieve orthogonality to a given order, our construction relies, in addition to $m$, on a moment condition for the nuisance parameter, $\mathbb{E}[g(W_i, \theta_0, \eta_0)] = 0$.\footnote{In our construction, the parameter $\theta$ plays a passive role. It is carried through as an argument of $m$ and $g$. Only after the orthogonal moment function has been constructed does $\theta$ enter the analysis as the parameter to be estimated, through the orthogonal moment equations. Thus, one could drop $\theta$ from most of the formulas in Sections (ref) through (ref) without changing their content, but we retain it throughout the main text to avoid notational confusion.}
Lastly, in our setup, $\eta_0$ is a vector of finite dimension (which may grow with the sample size). The possible extension to infinite-dimensional nuisance parameters, as in the literature on higher-order influence functions RobinsLiTchetgenTchetgenvanderVaart2008,vanderVaart2014, is left for future work.
Let $\partial_\eta$ be the gradient with respect to $\eta$. Our construction requires a left inverse of the $d_g \times d_\eta$ Jacobian matrix
Let $\Lambda_0$ be a $d_\eta \times d_g$ matrix satisfying $\Lambda_0 J_0 = I_{d_\eta}$. When $d_g = d_\eta$ and $J_0$ is invertible, $\Lambda_0 = J_0^{-1}$. When $d_g > d_\eta$ and $J_0$ has full column rank, we can set $\Lambda_0 = (J_0'\, \Omega^{-1}\, J_0)^{-1} J_0'\, \Omega^{-1}$ for a positive-definite weight matrix $\Omega$.
The matrix $\Lambda_0$ enters as an additional nuisance parameter in our construction. We parametrize it as $\Lambda(\lambda)$ for a finite-dimensional parameter $\lambda \in \mathbb{R}^{d_\lambda}$, with true value $\lambda_0$ satisfying $\Lambda(\lambda_0) = \Lambda_0$. The simplest choice is $\lambda = \mathrm{vec}(\Lambda)$, giving $d_\lambda = d_\eta \cdot d_g$.
While the dimension $d_\lambda$ does not enter the orthogonalization construction itself, $d_{\lambda}$ does matter in applications, because the accuracy with which $\widehat\Lambda$ can be estimated depends on $d_\lambda$, with a smaller $d_\lambda$ potentially leading to a more accurate preliminary estimate. For this reason, Section (ref) develops techniques that reduce $d_\lambda$, to a single scalar if desired, by replacing the moment function $g$ with a suitably transformed version.
We first briefly illustrate the setup with two cross-sectional examples, and we will study a grouped data model in more detail in the next section.
Let $\psi(W_1, \ldots, W_L; \theta, \eta, \lambda) \in \mathbb{R}$ be a moment function that depends on $L \geq 1$ independent copies of $W_i$, the parameter of interest $\theta$, the nuisance parameter $\eta$, and the Jacobian-inverse parameter $\lambda$. Define the population moment
When $q = 1$, Definition (ref) reduces to the standard Neyman orthogonality condition. When $q = 0$, no orthogonality is required and $\psi(W_i; \theta, \eta) = m(W_i, \theta, \eta)$ suffices (with $L=1$). The condition (ref) imposes orthogonality in $\eta$ and $\lambda$ jointly. Beyond the pure $\eta$-derivatives and pure $\lambda$-derivatives of $\Psi$, it also requires all mixed derivatives $\partial_\eta^\alpha \partial_\lambda^\beta \Psi$ with $|\alpha| \geq 1$, $|\beta| \geq 1$, and $|\alpha| + |\beta| \leq q$ to vanish at the true values. The pure $\eta$-conditions control bias from estimating the original nuisance parameter, the pure $\lambda$-conditions control bias from estimating the Jacobian inverse, and the mixed conditions control bias arising from the interaction of the two estimation errors in the Taylor expansion of $\Psi$. All three kinds of condition are necessary for $q$-th order orthogonality.
\paragraph{Estimation.} Given an i.i.d.\ sample $\{W_i : i = 1, \ldots, n\}$, the natural sample analogue of the population moment $\Psi(\theta, \eta, \lambda)$ is the U-statistic \[ \widehat\Psi_n(\theta, \eta, \lambda) := \frac{1}{n(n-1)\cdots(n-L+1)} \sum_{\substack{i_1, \ldots, i_L = 1 \\ \text{all distinct}}}^{n} \psi(W_{i_1}, \ldots, W_{i_L}; \theta, \eta, \lambda). \] An estimator $\widehat\theta$ of $\theta_0$ is obtained by solving $\widehat\Psi_n(\widehat\theta, \widehat\eta, \widehat\lambda) = 0$, where $\widehat\eta$ and $\widehat\lambda$ are preliminary estimators of the nuisance parameters. To ensure that the estimation error in $(\widehat\eta, \widehat\lambda)$ does not contaminate the U-statistic, the sample is split into two independent parts. One part is used to construct $\widehat\eta$ and $\widehat\lambda$, and the other is used to evaluate $\widehat\Psi_n$. Efficiency can be improved by cross-fitting, that is by swapping the roles of the two sample parts and averaging, as in ChernozhukovChetverikovDemirerDufloHansenNeweyRobins2018.
The key consequence of $q$-th order orthogonality is that, under appropriate regularity conditions, the estimator satisfies $\widehat\theta - \theta_0 = O_p(\|\widehat\eta - \eta_0\|^{q+1} + \|\widehat\lambda - \lambda_0\|^{q+1} + n^{-\nicefrac{1}{2}})$, since all lower-order terms in the Taylor expansion of $\Psi$ around $(\theta_0, \eta_0, \lambda_0)$ are of smaller order because of the orthogonality conditions. Valid inference is therefore possible even when the nuisance parameters converge at rates much slower than $n^{-\nicefrac{1}{4}}$, provided $q$ is chosen large enough. We provide more details on the asymptotic theory in Section (ref).
Before presenting the general construction of higher-order orthogonal moment functions, we illustrate the main ideas in a concrete example.
Consider a grouped data model with heterogeneous coefficients,
for a scalar outcome $Y_{it}$ and a $d_\eta$-dimensional covariate vector $X_{it}$, under the uncorrelatedness condition
Such models have been studied extensively under the mean independence condition $\mathbb{E}[U_{it} \mid X_{i1}, \ldots, X_{iT}] = 0$ (e.g., Chamberlain1992). Condition (ref) is weaker, since it does not restrict the conditional mean of $Y_{it}$ given $X_{it}$ to be linear. The coefficients $\eta_{i0}$ are the best linear predictors, \[ \eta_{i0} = \bigg(\mathbb{E}\bigg[\sum_{t=1}^T X_{it} X_{it}'\bigg]\bigg)^{-1} \mathbb{E}\bigg[\sum_{t=1}^T X_{it} Y_{it}\bigg]. \] We focus on the target parameter \[ \theta_0 = \frac{1}{N} \sum_{i=1}^N m(\eta_{i0}), \] where $m$ is a smooth scalar function. The mean and second moment of the components of $\eta_{i0}$ are natural choices, for example to quantify the level and dispersion of the coefficients across units. As an application, KlineRoseWalters2022 study employment discrimination by sending randomized fictitious job applications to firms and measuring callbacks $Y_{it}$. There, $\eta_{i0}$ captures how each firm $i$ responds to applicant characteristics $X_{it}$ such as race, gender, and age, and moments of $\eta_{i0}$ quantify the extent and variation of discrimination across firms.
This model fits the framework of Section (ref). For each unit $i$, write the observation at $t$ as $W_{it} = (Y_{it}, X_{it})$. The nuisance-identifying moment function is \[ g(W_{it}, \theta, \eta_i) = X_{it}(Y_{it} - X_{it}' \eta_i), \] and condition (ref) is equivalent to $\mathbb{E}[g(W_{it}, \theta_0, \eta_{i0})] = 0$ for every $t$. The Jacobian-inverse nuisance parameter is $\Lambda_{i0} = -(\mathbb{E}[X_{it} X_{it}'])^{-1}$, and the target-moment function for a given unit $i$ is $m_i(W_{it}, \theta, \eta_i) = m(\eta_i) - \theta$.
The observations $W_{it}$ (i.e., job applications) are independent across $t$ for each $i$ (i.e., within each firm). The moment conditions for $\eta_{i0}$ are affine in $\eta_i$, so this example falls within the affine case of Section (ref) below. The ordinary least-squares estimators \[ \widehat\eta_i^{\rm OLS} = \bigg(\sum_{t=1}^T X_{it} X_{it}'\bigg)^{-1} \sum_{t=1}^T X_{it} Y_{it},\quad i=1,...,N, \] are noisy when $T$ is small, and the plug-in estimator $\widehat\theta^{\rm plug\text{-}in} = (1/N) \sum_{i=1}^N m(\widehat\eta_i^{\rm OLS})$ is biased. Our approach constructs alternative estimators with reduced bias.
To apply the framework from Section (ref) to this setting, let us fix a single unit $i$ and consider a sample of observations for $t = 1, \ldots, T$. Each observation $W_{it}$ plays the role of $W_i$ in Section (ref), the pair $(\eta_{i0}, \Lambda_{i0})$ plays the role of $(\eta_0, \Lambda_0)$, and the orthogonal moment function $\psi$ becomes a function of several observations for the same unit. The U-statistic constructions that follow are therefore applied across $t$ within a unit. The unit dimension $i = 1, \ldots, N$ plays no role at the orthogonalization stage. It enters only through the final step, in which the unit-level orthogonal moment functions are averaged over $i$ to form the estimator of $\theta_0$.\footnote{An equivalent approach is to stack $(\eta_{10}, \ldots, \eta_{N0})$ into a single nuisance vector of dimension $N d_\eta$ and treat the full data set of $n = NT$ observations as one sample to which the Section (ref) construction is applied directly. The two routes produce the same orthogonal moment functions. We adopt the unit-by-unit view because it keeps the effective nuisance dimension bounded as $N$ grows and because it is conceptually simpler.}
Since the problem stratifies across units, we focus on a single unit $i$ and describe the construction at the unit level. It is useful to introduce the reparameterization $A_{i0} = -\Lambda_{i0}^{-1}$ and $b_{i0} = -\Lambda_{i0}^{-1} \eta_{i0}$, with $a_{i0} = \mathrm{vech}(A_{i0})$ denoting the vector that contains the elements of $A_{i0}$.
The construction relies on three observations.
The first observation is that $A_{i0}$ and $b_{i0}$ admit simple unbiased estimators:
The second observation is that any polynomial $P(a_{i0}, b_{i0})$ admits an unbiased estimator. Each element of $(a_{i0}, b_{i0})$ has an unbiased estimator by (ref). Since observations are independent across $t$, one can combine estimators from different $t$ observations to produce unbiased estimators of powers and products. For example, when $d_\eta = 1$, an unbiased estimator of $a_{i0}^2 = (\mathbb{E}[X_{it}^2])^2$ is $X_{i1}^2 X_{i2}^2$, since $\mathbb{E}[X_{i1}^2 X_{i2}^2] = \mathbb{E}[X_{i1}^2]\, \mathbb{E}[X_{i2}^2] = a_{i0}^2$.
The third observation is that $m(\eta_{i0})$ can be approximated by a polynomial in $(a_{i0}, b_{i0})$. Writing $m(\eta_{i0}) = \varphi(a_{i0}, b_{i0})$ with $\varphi(a, b) = m(A^{-1} b)$, the function $\varphi$ has a singularity at $\det A = 0$ and so cannot be globally represented as a polynomial. However, given preliminary estimators $\widehat\Lambda_i$ and $\widehat\eta_i$ that are independent of the remaining data, one can expand locally: setting $\widehat a_i = -\mathrm{vec}(\widehat\Lambda_i^{-1})$ and $\widehat b_i = -\widehat\Lambda_i^{-1} \widehat\eta_i$,
where the remainder $R_{q,i}$ depends on powers of order $q+1$ in $a_{i0} - \widehat a_i$ and $b_{i0} - \widehat b_i$. Each term $(a_{i0} - \widehat a_i)^\alpha (b_{i0} - \widehat b_i)^\beta$ admits an unbiased estimator built from independent observations. This yields an approximately unbiased estimator of $m(\eta_{i0})$ based on a $q$-th order orthogonal moment function.
In practice, the preliminary estimators $\widehat\Lambda_i$ and $\widehat\eta_i$ are obtained by holding out a set of $t$ observations. Efficiency can be improved by cross-fitting. A feature of this construction is that the nuisance parameters $\widehat\Lambda_i$ and $\widehat\eta_i$ do not change with the order of orthogonalization, $q$.
Consider first the special case without covariates,
where $U_{it}$ are i.i.d.\ but not necessarily normally distributed. This is the NeymanScott1948 model without the normality assumption. Here $d_\eta = 1$ and $\Lambda_{i0} = -1$ is known, so $\eta_{i0}$ is the only nuisance parameter.
A $q$-th order Taylor expansion of $m(\eta_{i0})$ around a preliminary estimator $\widehat\eta_i$ gives \[ m(\eta_{i0}) = \sum_{k=0}^q \frac{\partial_\eta^k m(\widehat\eta_i)}{k!}\, (\eta_{i0} - \widehat\eta_i)^k + R_{q,i}. \] Expanding $(\eta_{i0} - \widehat\eta_i)^k = \sum_{j=0}^k \binom{k}{j} \eta_{i0}^j\, (-\widehat\eta_i)^{k-j}$ and interchanging the order of summation, \[ m(\eta_{i0}) = \sum_{j=0}^q \bigg(\sum_{k=j}^q \frac{\partial_\eta^k m(\widehat\eta_i)}{k!}\, \binom{k}{j}\, (-\widehat\eta_i)^{k-j}\bigg)\, \eta_{i0}^j + R_{q,i}. \] Since, by independence and $\mathbb{E}[U_{it}] = 0$,
an approximately unbiased estimator of $m(\eta_{i0})$ is
A lower-variance estimator is obtained by averaging (ref) over all subsets of $q$ observations, yielding a U-statistic.
\paragraph{Orthogonality.} Define the moment function
where $W_i = (Y_{i1}, \ldots, Y_{iT})$. Using ((ref)), reversing the order of summation, and applying the binomial theorem, we have
At $(\theta_{i0}, \eta_{i0})$ with $\theta_{i0} = m(\eta_{i0})$, the right-hand side equals $m(\eta_{i0}) - \theta_{i0} = 0$. Write $f(\eta_i) = \sum_{k=0}^q \frac{\partial_\eta^k m(\eta_i)}{k!}\, (\eta_{i0} - \eta_i)^k$. A direct computation using the product rule shows that the derivative telescopes:
which vanishes at $\eta_i = \eta_{i0}$. By Leibniz' rule, the $j$-th derivative of $f$ at $\eta_i = \eta_{i0}$ also vanishes for all $1 \leq j \leq q$, since every term contains a factor $(\eta_{i0} - \eta_i)^{q+1-j+\ell}$ for some $\ell \geq 0$. The moment function $\psi$ is, therefore, $q$-th order Neyman-orthogonal in the sense of Definition (ref).
We now consider the model (ref) with a scalar covariate ($d_\eta = 1$) and construct a second-order orthogonal moment function ($q = 2$). Unlike the previous example, the Jacobian inverse $\Lambda_{i0} =- (\mathbb{E}[X_{it}^2])^{-1}$ is now unknown and enters as an additional nuisance parameter.
A second-order Taylor expansion of $m(\eta_{i0})$ around $\widehat\eta_i$ gives
To express $\eta_{i0} - \widehat\eta_i$ in terms of quantities with unbiased estimators, write $a_{i0} =- \Lambda_{i0}^{-1} = \mathbb{E}[X_{it}^2]$ and $b_{i0} = -\Lambda_{i0}^{-1} \eta_{i0} = \mathbb{E}[X_{it} Y_{it}]$. Expanding $\eta_{i0} = b_{i0}/a_{i0}$ around $(\widehat a_i, \widehat b_i) = (-\widehat\Lambda_i^{-1}, -\widehat\Lambda_i^{-1}\widehat\eta_i)$ to second order, \[ \eta_{i0} - \widehat\eta_i \approx -\widehat\Lambda_i(b_{i0} - \widehat\eta_i\, a_{i0}) - \widehat\Lambda_i(1 + \widehat\Lambda_i\, a_{i0})(b_{i0} - \widehat\eta_i\, a_{i0}). \] Substituting into ((ref)) and replacing $a_{i0}$ and $b_{i0}$ by their unbiased estimators ($X_{it}^2$ and $X_{it} Y_{it}$, respectively), we obtain the approximately unbiased estimator of $m(\eta_{i0})$
The resulting moment function is
One can verify that $\mathbb{E}[\psi(W_i; \theta_{i0}, \eta_i, \Lambda_i)]$ and all its first and second derivatives in $(\eta_i, \Lambda_i)$ at $(\eta_{i0}, \Lambda_{i0})$ vanish. The moment function is therefore second-order Neyman-orthogonal.
Both constructions in this section are special cases of the general affine formula in Theorem (ref) below. The next section presents the general construction for arbitrary $d_\eta$ and arbitrary order $q$, first for affine $g$ and then for the fully nonlinear case.
We now present the general construction of higher-order orthogonal moment functions. The target moment function $m(W_i, \theta, \eta)$ is allowed to be nonlinear in $\eta$ throughout. We first treat the case where the nuisance-identifying moment function $g(W_i, \theta, \eta)$ is affine in $\eta$, and then extend to the case where $g$ is also nonlinear.
Assume throughout this subsection that the moment function $g(W_i, \theta, \eta)$ is affine in $\eta$, so that $\partial_\eta g(W_i, \theta, \eta)$ does not depend on $\eta$, whereas the target moment function $m(W_i, \theta, \eta)$ can be nonlinear in $\eta$. This covers many important cases. Linear regression and linear instrumental variables (Example (ref)) have $g$ and $m$ both affine in $\eta$, while the average effect of a generated regressor (Example (ref)) has $g$ affine and $m$ nonlinear in $\eta$.
We begin by introducing notation for multilinear derivatives that is used throughout this section. For $r \geq 1$, write $\partial_\eta^r m(W_i, \theta, \eta)$ for the $r$-th derivative tensor of $m$ with respect to $\eta$. For vectors $v^{(1)}, \ldots, v^{(r)} \in \mathbb{R}^{d_\eta}$, denote the multilinear action by \[ \langle \partial_\eta^r m(W_i, \theta, \eta),\, v^{(1)} \otimes \cdots \otimes v^{(r)} \rangle = \sum_{j_1, \ldots, j_r} \frac{\partial^r m(W_i, \theta, \eta)}{\partial \eta_{j_1} \cdots \partial \eta_{j_r}} v^{(1)}_{j_1} \cdots v^{(r)}_{j_r}. \] For $r = 1$ this gives the directional derivative $\partial_\eta m(W_i, \theta, \eta)' v^{(1)}$. For $r = 2$ it gives the quadratic form $v^{(1)\prime} \partial_{\eta\eta'}^2 m(W_i, \theta, \eta)\, v^{(2)}$.
For an integer $q \geq 0$, the orthogonal moment function is built from factors $Z_k$, each using $k$ independent copies of $W_i$. The outer sum aggregates over all ways to distribute at most $q$ copies among $r$ such factors, with the $r$-th derivative of $m$ providing the contraction at the root. Concretely, define
where the tensor product $\bigotimes_{s=1}^r Z_{k_s}$ denotes $Z_{k_1} \otimes \cdots \otimes Z_{k_r}$ (i.e., $\langle \partial_\eta^r m, \bigotimes_s Z_{k_s} \rangle$ is the $r$-linear contraction of $\partial_\eta^r m$ with the vectors $Z_{k_1}, \ldots, Z_{k_r}$), the argument $W = (W_1, W^{(1)}, \allowbreak \ldots , W^{(q)})$ consists of $q+1$ independent copies of $W_i$, and
Each block $W^{(s)} = (W_{(s,1)}, \ldots, W_{(s,k_s)})$ uses $k_s$ independent copies, the copy $W_1$ enters only through $\partial_\eta^r m(W_1, \theta, \eta)$, and the $r = 0$ term uses the empty-product convention.
Each $Z_{k_s}$ pre-multiplies $k_s - 1$ Jacobian factors $\Lambda \partial_\eta g$ onto a terminal factor $\Lambda g$. The binomial weight $(-1)^{\sum_{s=1}^r k_s} \binom{q}{\sum_{s=1}^r k_s}$ ensures that all monomials of degree $1$ through $q$ in the nuisance estimation error cancel in expectation.
\paragraph{Explicit formulas for small $q$.} To make formula (ref) concrete, we display $\psi^{(q)}$ for $q = 0, 1, 2$. Write $m := m(W_1, \theta, \eta)$, $m_\eta := \partial_\eta m(W_1, \theta, \eta)$, and $m_{\eta\eta} := \partial_{\eta\eta'}^2 m(W_1, \theta, \eta)$.
\paragraph{Case $q = 0$.} \[ \psi^{(0)}(W_1; \theta, \eta) = m(W_1, \theta, \eta). \]
\paragraph{Case $q = 1$.} \[ \psi^{(1)}(W_1, W_2; \theta, \eta, \lambda) = m - m_\eta' \,\Lambda(\lambda)\, g(W_2, \theta, \eta). \]
\paragraph{Case $q = 2$.}
The $q = 2$ moment function uses $L = 3$ independent copies. Explicit formulas for $q = 3$ and $q = 4$, using $L = 4$ and $L = 5$ copies respectively, are given in Appendix (ref).
The proof is in Appendix (ref).
When the moment function $g(W_i, \theta, \eta)$ is nonlinear in $\eta$, additional correction terms are needed beyond the affine orthogonal moment function. The target $m$ was already allowed to be nonlinear in the previous subsection, but we now also allow for nonlinearity of $g$. We denote the affine orthogonal moment function from (ref) by $\psi^{(q)}_{\rm aff}$ and write the corrected orthogonal moment function as
The reason $\psi^{(q)}_{\rm aff}$ fails to be orthogonal for nonlinear $g$ is that the proof of Theorem (ref) relies on $\Lambda_0\, \mathbb{E}[g(W_i, \theta_0, \eta_0 + \delta)]$ being linear in $\delta$. When $g$ is nonlinear, this expression contains quadratic and higher-order terms in $\delta$ that the affine construction does not cancel. The correction terms involve the higher-order derivatives $\partial_\eta^p g(W_i, \theta, \eta)$ for $p \geq 2$. For $q = 1$ no correction is needed: $\psi^{(1)} = \psi^{(1)}_{\rm aff}$ is first-order Neyman-orthogonal even when $g$ is nonlinear.
Before turning to the construction, we set up notation used throughout this subsection. For $p \geq 1$, the $p$-th derivative tensor $\partial_\eta^p g(W_i, \theta, \eta) \in \mathbb{R}^{d_g} \otimes (\mathbb{R}^{d_\eta})^{\otimes p}$ has one $d_g$-output slot and $p$ input slots in $\mathbb{R}^{d_\eta}$, with entries $[\partial_\eta^p g]_{k, j_1, \ldots, j_p} = \partial^p g_k(W_i, \theta, \eta) / \partial \eta_{j_1} \cdots \partial \eta_{j_p}$. For $p = 1$ this is the $d_g \times d_\eta$ Jacobian. To express contractions of such tensors with vectors in $\mathbb{R}^{d_\eta}$, we extend the angle-bracket notation of Section (ref) as follows: for a tensor $A$ with $p$ input slots in $\mathbb{R}^{d_\eta}$ and output in a vector space $V$, and for $v_1, \ldots, v_p \in \mathbb{R}^{d_\eta}$, \[ A[v_1, \ldots, v_p] \;\in\; V \] denotes the contraction of $A$ along its inputs, leaving the $V$-output unchanged. For the scalar function $m$, the bracket and the angle-bracket of Section (ref) agree, $\partial_\eta^r m\,[v_1, \ldots, v_r] = \langle \partial_\eta^r m, v_1 \otimes \cdots \otimes v_r \rangle \in \mathbb{R}$. For the vector-valued $g$, $\partial_\eta^p g\,[v_1, \ldots, v_p] \in \mathbb{R}^{d_g}$, and pre-multiplying by $\Lambda(\lambda)$ contracts the $d_g$-output slot to give a vector in $\mathbb{R}^{d_\eta}$.
The orthogonal moment function for the general nonlinear case can be derived from the affine formula via a change of nuisance parameter $\phi = f(\eta)$, and the goal of this subsection is to show that derivation for $q = 2$. Let $f : \mathbb{R}^{d_\eta} \rightarrow \mathbb{R}^{d_\eta}$ be a parameter transformation, define $ \overline g(W_i, \theta, \phi) := g\!\left(W_i, \theta, f^{-1}(\phi)\right). $
We want to apply the affine moment function construction to $\overline g$, but it is important to note that $\overline g$ being affine in $\phi$ is not a necessary condition for that construction to be applicable. Inspecting the proof of Theorem (ref) in the appendix, one sees that all that is actually needed at $q = 2$ is the condition
where $\partial_\phi^2 \overline g(W_i, \theta, \phi) \in \mathbb{R}^{d_g} \otimes (\mathbb{R}^{d_\eta})^{\otimes 2}$ follows the same convention as $\partial_\eta^p g$ above, with $\eta$ replaced by $\phi$ and $g$ by $\overline g$. Pre-multiplication by $\Lambda_0$ contracts the $d_g$-output slot, and (ref) requires the resulting $d_\eta \times d_\eta \times d_\eta$ tensor to vanish. Define $T_2 := \Lambda_0\, \mathbb{E}\!\left[\partial_\eta^2 g(W_i, \theta_0, \eta_0)\right] \;\in\; \mathbb{R}^{d_\eta \times d_\eta \times d_\eta}$. Then a simple explicit choice for $f$ that guarantees (ref) is given by
Because (ref) holds for this choice of $f$,\footnote{From (ref) we have $\phi_0 = \eta_0$ and $\partial_\eta f(\eta_0) = I_{d_\eta}$, so by the inverse function theorem $f^{-1}$ exists locally with $\partial_\phi f^{-1}(\phi_0) = I_{d_\eta}$ and $\partial_\phi^2 f^{-1}(\phi_0) = -T_2$. The chain rule applied to $\overline g(W_i, \theta_0, \phi) = g(W_i, \theta_0, f^{-1}(\phi))$ at $\phi_0$ then gives $\partial_\phi^2 \overline g(W_i, \theta_0, \phi_0) = \partial_\eta^2 g(W_i, \theta_0, \eta_0) + \partial_\eta g(W_i, \theta_0, \eta_0)\cdot(-T_2)$, and pre-multiplying by $\Lambda_0$ and taking expectations yields $T_2 - T_2 = 0$.} the affine $q = 2$ moment function from (ref), with $g$ replaced by $\overline g$ and $\eta$ replaced by $\phi$, is second-order Neyman-orthogonal at $(\theta_0, \phi_0, \lambda_0)$. Defining $\overline m(W_i, \theta, \phi) := m(W_i, \theta, f^{-1}(\phi))$ analogously to $\overline g$, this orthogonal moment function reads, evaluated at $\phi = f(\eta)$,
By construction $\overline\psi^{(2)}$ is second-order Neyman-orthogonal at $(\theta_0, \eta_0, \lambda_0)$, but its expression in $\eta$-coordinates still involves $\overline m, \overline g$ and their $\phi$-derivatives, which is not yet a closed-form expression in $\eta$. We obtain such an expression by chain-rule expansion of each $\overline m, \overline g$ derivative at $\phi = f(\eta_0) = \eta_0$, using $\partial_\phi f^{-1}(\phi_0) = I_{d_\eta}$ and $\partial_\phi^2 f^{-1}(\phi_0) = -T_2$. The first three terms of $\overline\psi^{(2)}$ reduce to the corresponding terms of $\psi^{(2)}_{\rm aff}$, and the second-order chain rule for $\overline m$ produces the identity \[ \partial_{\phi\phi'}^2 \overline m(W_1, \theta_0, \phi_0) \;=\; \partial_{\eta\eta'}^2 m(W_1, \theta_0, \eta_0) \;-\; \partial_\eta m(W_1, \theta_0, \eta_0)\cdot T_2, \] which splits the fourth term of $\overline\psi^{(2)}$ into the fourth term of $\psi^{(2)}_{\rm aff}$ plus a $T_2$-correction. The result is that the value of $\overline\psi^{(2)}$ at the truth $(\theta_0, \eta_0, \lambda_0)$, and its $\eta$- and $\lambda$-derivatives at the truth up to total order two, coincide with those of the moment function\footnote{The two functions $\overline\psi^{(2)}$ and $\psi^{(2)}$ are not equal as functions of $(\eta, \lambda)$ --- the chain-rule identities used to derive (ref) hold only at $\eta = \eta_0$, since $\partial_\phi f^{-1}(f(\eta))$ depends on $\eta$ for $\eta \neq \eta_0$ --- but they have matching Taylor expansions at the truth up to order two, which is what orthogonality requires.}
However, the expression for $\psi^{(2)}(W, \theta, \eta, \lambda)$ in the last display is not yet a usable moment function, because $T_2$ is unknown. To make it feasible we replace $T_2 = \Lambda_0\, \mathbb{E}[\partial_\eta^2 g(W_4, \theta_0, \eta_0)]$ by $\Lambda(\lambda)\, \partial_\eta^2 g(W_4, \theta, \eta)$, where $W_4$ is an independent copy of the data. This delivers a moment function that uses $L = 4$ independent copies of $W_i$ and that by construction is second-order Neyman-orthogonal:
where we swapped the labeling of $W_1,\ldots,W_4$ in the final term. The replacement of $T_2$ by $\Lambda(\lambda)\, \partial_\eta^2 g(W_i, \theta, \eta)$ introduces additional dependence on $(\eta, \lambda)$ that was not present in the $\phi$-coordinate derivation, but this does not affect second-order Neyman-orthogonality.
We now introduce notation that turns out to be useful for expressing $\psi^{(q)}$ at general $q$. We motivate it by rewriting $\psi^{(2)}$ from (ref). Figure (ref) shows how the five terms of $\psi^{(2)}$ can each be represented by a rooted tree built from three kinds of nodes. The root (shaded blue) carries either $m$ or one of its $\eta$-derivatives $\partial_\eta^r m$, where $r$ is the number of children of the root. The leaves (shaded gray) each carry $\Lambda\, g$. The non-root, non-leaf nodes (unshaded) each carry $\Lambda\, \partial_\eta^p g$, where $p$ is the number of children of the node. Each node is evaluated on its own independent copy of $W_i$. The tree, which we denote by $\tau$, then uniquely determines the corresponding term in $\psi^{(2)}$, which we denote by $\kappa_\tau(W, \theta, \eta, \lambda)$.
It turns out that all terms in the orthogonal moment function $\psi^{(q)}$ can be indexed by rooted trees in the same way for any order $q \in \{1, 2, 3, \ldots\}$,
where $\mathcal{T}_q$ is a finite set of rooted trees, $c_{q,\tau} \in \mathbb{R}$ is a combinatorial coefficient, and the kernel $\kappa_\tau(W, \theta, \eta, \lambda)$ is determined by $\tau$ via the rules illustrated in Figure (ref).
A rooted tree $\tau$ is a finite tree (a connected acyclic graph) with one node singled out as the root. Let $d(\tau)$ denote the number of non-root nodes of $\tau$ with at most one child. The middle column of Figure (ref) reports $d(\tau)$ for each of the five trees: the trivial tree has $d = 0$, the root with one leaf child has $d = 1$, and the remaining three trees all have $d = 2$. The set $\mathcal{T}_q$ in (ref) is then
For $q=0$, $\mathcal{T}_0$ only contains the single-node tree. For $q=1$, $\mathcal{T}_1$ contains the single-node tree and the tree with one root and one leaf child. For $q=2$, the five elements of $\mathcal{T}_2$ are shown in Figure (ref). For $q = 3$ we have $|\mathcal{T}_3| = 13$, shown in Figure (ref). The cardinalities grow rapidly, with $|\mathcal{T}_4| = 40$ and $|\mathcal{T}_5| = 130$.
Each rooted tree $\tau \in \mathcal{T}_q$ translates into a kernel $\kappa_\tau(W, \theta, \eta, \lambda)$ by populating the nodes with the algebraic objects already used in $\psi^{(2)}$. The root carries $m$ or one of its $\eta$-derivatives, with the order of differentiation equal to the number of children of the root. Each leaf carries $\Lambda\, g$. Each non-root, non-leaf node carries $\Lambda\, \partial_\eta^p g$, where $p$ is the number of children of the node. Each node is evaluated on its own independent copy of $W_i$, and the resulting expression is the contraction along the parent--child edges of the tree, exactly as in Figure (ref). Applying this construction to the seven trees in the top row of Figure (ref) gives the affine moment function $\psi^{(3)}_{\rm aff}$, where for compactness we suppress the arguments $(\theta, \eta)$ in $m$ and $g$ and the argument $\lambda$ in $\Lambda$:
while the six trees in the bottom row of Figure (ref) correspond to the additional correction terms that are needed when $g$ is nonlinear in $\eta$,
The combinatorial pre-factors in (ref) and (ref) are exactly the values needed for third-order Neyman orthogonality. The affine pre-factors can be obtained from (ref) by collecting the ordered chain-length tuples $(k_1, \ldots, k_r)$ that correspond to the same unordered tree. The correction term pre-factors are the values of the closed-form $c_{q,\tau}$ that we introduce next. Note that the substantive content of the tree formalism is our closed form expression for $c_{q,\tau}$, while the index set $\mathcal{T}_q$ and the construction rule for $\kappa_\tau$ are purely notational bookkeeping tools.
To state the closed form for $c_{q,\tau}$, we first need to define two further integers associated with each rooted tree $\tau$. We have already introduced $d(\tau)$ as the number of non-root nodes of $\tau$ with at most one child. Let $|\tau|$ denote the total number of non-root nodes of $\tau$ (equivalently, the number of edges of $\tau$). Let $\mathrm{Aut}(\tau)$ denote the automorphism group of $\tau$, that is, the set of bijections of the nodes of $\tau$ that fix the root and keep the tree structure unchanged. Finally, let $|\mathrm{Aut}(\tau)|$ be its order (i.e.\ the cardinality of $\mathrm{Aut}(\tau)$). Then,
Figure (ref) shows the values of $|\tau|$, $d(\tau)$, and $|\mathrm{Aut}(\tau)|$ for all rooted trees $\tau$ with $d(\tau) \leq 3$. The values of $|\tau|$ and $d(\tau)$ are read off by counting the relevant nodes. The values of $|\mathrm{Aut}(\tau)|$ in this figure are also straightforward: for many of the trees, only the trivial bijection preserves the tree structure, so $|\mathrm{Aut}(\tau)| = 1$. For each of the remaining trees, exactly one node has $c$ isomorphic leaf children, which can be permuted in $c!$ ways, giving $|\mathrm{Aut}(\tau)| = c!$. For larger values of $d(\tau)$ the automorphism group can of course be more complicated.
As a check, applying (ref) at $q = 3$ to the values in Figure (ref) reproduces the 13 coefficients in (ref) and (ref). For example, the root with three leaf children gives $c_{3,\tau} = -\binom{3}{3}/3! = -1/6$ (matching the coefficient of $\partial_\eta^3 m[\Lambda g, \Lambda g, \Lambda g]$ in (ref)), etc.\footnote{ It is also instructive to verify the consistency between (ref) and the affine moment function in (ref). The trees that contribute to $\psi^{(q)}_{\rm aff}$ are precisely those in which only the root has more than one child, equivalently those with $|\tau| = d(\tau)$. For such a tree, let $r \geq 0$ be the number of children of the root, and let $k_s \geq 1$ be the chain length of each child for $s \in \{1, \ldots, r\}$. We then have $|\tau| = d(\tau) = k_1 + \cdots + k_r$. The automorphism group permutes branches of equal length, so if the multiset $\{k_1, \ldots, k_r\}$ has $n_\ell$ elements equal to $\ell$, then $|\mathrm{Aut}(\tau)| = \prod_\ell n_\ell!$, and (ref) simplifies to
The same unordered tree $\tau$ is the image of $r! / \prod_\ell n_\ell!$ ordered branch tuples $(k_1, \ldots, k_r)$, each contributing the term $\frac{(-1)^{|\tau|}}{r!}\, \binom{q}{|\tau|}\, \kappa_\tau$ to the affine moment function $\psi^{(q)}_{\rm aff}$ in (ref). Summing across these tuples cancels the $r!$ and reproduces exactly $c_{q,\tau}\, \kappa_\tau$ with $c_{q,\tau}$ as in (ref), confirming that $\psi^{(q)}_{\rm aff}$ in (ref) matches the affine part of (ref). }
Tree-indexed expansions of this form, including the appearance of the symmetry factor $1/|\mathrm{Aut}(\tau)|$, are a well-known structure in the literature on multivariate Taylor expansions, see Appendix (ref) for a discussion of the Butcher $B$-series. That appendix also gives a simple recursion for computing $|\mathrm{Aut}(\tau)|$.\footnote{ Using this recursion we have numerically verified that (ref) delivers orthogonal moment functions for $q$ up to $10$ in the scalar case $d_\eta = d_g = 1$, and for $q$ up to $8$ in the multivariate case. Enumerating all $|\mathcal{T}_{10}| = 110,135$ trees and calculating all corresponding $c_{q,\tau}$ is possible within seconds.}
With $c_{q,\tau}$ in hand, we can state our main result for the nonlinear case.
The proof is in Appendix (ref).
The expression (ref) applies to both the affine and correction parts of $\psi^{(q)}$. Trees with $|\tau| = d(\tau)$ are exactly those in which every non-root node is either a leaf or a one-child node; for these, (ref) reduces to ((ref)), the binomial weights of the affine score (ref). Trees with $|\tau| > d(\tau)$ --- those with at least one non-root node with at least two children --- are the correction trees, which contribute terms involving the higher-order derivatives $\partial_\eta^p g$ for $p \geq 2$. When $g$ is affine in $\eta$, all such kernels vanish, (ref) reduces to $\psi^{(q)} = \psi^{(q)}_{\rm aff}$, and Theorem (ref) recovers Theorem (ref).
The total number of independent copies of $W_i$ used by $\psi^{(q)}$ equals one plus the maximum number of nodes in any tree $\tau \in \mathcal{T}_q$, which is bounded by $L \leq 2q$ for $q \geq 2$ (compared with $L = 1 + q$ in the affine case). The bound is attained by the binary tree with $d(\tau) = q$ in which every non-root, non-leaf node has exactly two children.
To reduce the dimension of the nuisance parameter $\lambda$, the key observation is that the moment function $g(W_i, \theta, \eta)$ in (ref) can be replaced by a different function $\widetilde g(W, \theta, \eta)$ that may depend on multiple independent copies of $W_i$.\footnote{As before, here we write $W$ without subscript when the argument requires several copies.} The replacement is valid as long as $\mathbb{E}[\widetilde g(W, \theta_0, \eta_0)] = 0$ and the corresponding Jacobian $\widetilde J_0 = \mathbb{E}[\partial_\eta \widetilde g(W, \theta_0, \eta_0)]$ has a left inverse $\Lambda_0$ with $\Lambda_0 \widetilde J_0 = I_{d_\eta}$. The orthogonality construction in Sections (ref)--(ref) applies with $g$ replaced by $\widetilde g$ throughout. The advantage is that $\widetilde g$ can be designed so that the nuisance parameter $\lambda$ has much lower dimension.
There is, however, a trade-off: reducing $d_\lambda$ comes at the cost of increasing the number of independent copies $L$ needed by the orthogonal moment function. The original moment function $g(W_i, \theta, \eta)$ depends on a single observation, so the $q$-th order affine orthogonal moment function uses $L = 1 + q$ copies. The transformed function $\widetilde g(W, \theta, \eta)$ may itself require multiple copies (e.g., $d_\eta$ copies for the determinant construction below), and each appearance of $\widetilde g$ in the orthogonal moment function then uses that many copies, increasing $L$ accordingly. In practice, the choice of $\widetilde g$ balances the difficulty of estimating a high-dimensional nuisance $\lambda$ (large $d_\lambda$, small $L$) against the variance cost of a higher-order U-statistic (small $d_\lambda$, large $L$).
In Sections (ref) and (ref), the orthogonal moment functions were stated in terms of a moment function $g$ that depends on a single observation $W_j$. When $\widetilde g$ requires a block of several independent copies, the construction is applied in exactly the same way, with one bookkeeping change: each occurrence of $g(W_j, \theta, \eta)$ or $\partial_\eta g(W_j, \theta, \eta)$ in the formulas of those sections is replaced by $\widetilde g$ or $\partial_\eta \widetilde g$ evaluated on a fresh block of independent copies of the size required by $\widetilde g$. We give examples of constructions for $\widetilde g$ below.
\paragraph{Determinant construction (exactly identified case).} Suppose $d_g = d_\eta$ and $J_0$ is invertible. Define $\widetilde g(W_1, \ldots, W_{d_\eta}; \theta, \eta)$ componentwise by
where $r = 1,\ldots,d_\eta$, and the $d_\eta \times d_\eta$ matrix inside the determinant has $g(W_r, \theta, \eta)$ in its $r$-th column and $\partial_{\eta_j} g(W_j, \theta, \eta)$ in its $j$-th column for $j \neq r$. Then $\mathbb{E}[\widetilde g_r] = 0$ follows from $\mathbb{E}[g(W_i, \theta_0, \eta_0)] = 0$ and the multilinearity of the determinant. The Jacobian satisfies $\mathbb{E}[\partial_\eta \widetilde g(W, \theta_0, \eta_0)] = \det(J_0) \cdot I_{d_\eta}$, so $\Lambda(\lambda) = \lambda \cdot I_{d_\eta}$ with \[ \lambda_0 = \frac{1}{\det(J_0)}. \] This gives $d_\lambda = 1$. The trade-off is that each evaluation of $\widetilde g$ uses $d_\eta$ independent copies of $W_i$, so the total number of copies for the $q$-th order moment function becomes $L = 1 + q \cdot d_\eta$.
\paragraph{Determinant construction (overidentified case).} When $d_g > d_\eta$, the Jacobian $J_0$ is no longer square, so the construction above does not directly apply. The idea is to first project the $d_g$-dimensional vectors $g$ and $\partial_{\eta_j} g$ down to $d_\eta$ dimensions using independent copies of $\partial_\eta g$, and then take the determinant. Concretely, define $\widetilde g(W_1, \ldots, W_{2d_\eta}; \theta, \eta)$ componentwise by
where $r = 1, \ldots, d_\eta$. Each column of the $d_\eta \times d_\eta$ matrix is formed by pre-multiplying a $d_g$-vector ($\partial_{\eta_j} g$ or $g$) by the $d_\eta \times d_g$ matrix $(\partial_\eta g)'$ from an independent copy, yielding a $d_\eta$-vector. The $r$-th column uses $g(W_{2r}, \theta, \eta)$ (without a $\partial_{\eta_r}$ derivative), while all other columns use $\partial_{\eta_j} g(W_{2j}, \theta, \eta)$. Then $\mathbb{E}[\widetilde g_r] = 0$ by the same multilinearity argument as before, the Jacobian satisfies $\mathbb{E}[\partial_\eta \widetilde g(W, \theta_0, \eta_0)] = \det(J_0' J_0) \cdot I_{d_\eta}$, and \[ \Lambda(\lambda) = \lambda \cdot I_{d_\eta}, \qquad \lambda_0 = \frac{1}{\det(J_0' J_0)}, \] again giving $d_\lambda = 1$. Each evaluation of $\widetilde g$ now uses $2 d_\eta$ independent copies of $W_i$.
In this section we present an implementation of our orthogonal moment functions in the setup of Section (ref), together with Monte Carlo simulation results. The simulation design is a stylized version of the setup in KlineRoseWalters2022 (though we do not attempt to calibrate the design to their empirical results). There are $N$ firms, and $T$ job applications per firm. Applications feature the race $X_{it,1}\in\{0,1\}$ of applicants and their gender $X_{it,2}\in\{0,1\}$. Binary callback indicators are generated as
where $\varepsilon_{it}$ are i.i.d. standard logistic independent of $X_{it,1},X_{it,2}$, and observations are independent across $i$ and $t$.
Job applications are randomly generated within firms but -- for example due to implementation constraints -- assignment is heterogeneous across firms. Specifically, we assume that firms belong to one of two types $Z_i\in\{1,2\}$, where $Z_i$ are i.i.d. Bernoulli with probability 1/2. When $Z_i=1$, $(X_{it,1}=1,X_{it,2}=1)$ and $(X_{it,1}=0,X_{it,2}=0)$ are each drawn with probability 3/8, and $(X_{it,1}=0,X_{it,2}=1)$ and $(X_{it,1}=1,X_{it,2}=0)$ are each drawn with probability 1/8. When $Z_i=2$, $(X_{it,1}=1,X_{it,2}=1)$ and $(X_{it,1}=0,X_{it,2}=0)$ are each drawn with probability 1/8, and $(X_{it,1}=0,X_{it,2}=1)$ and $(X_{it,1}=1,X_{it,2}=0)$ are each drawn with probability 3/8. For $j\in\{0,1,2\}$, $\beta_{ij}=Z_i+V_{ij}$, where the $V_{ij}$ are independent standard normal, independent of $Z_i$.
Let $X_{it}=(1,X_{it,1},X_{it,2})'$, and let $$\eta_{i0}=\left(\mathbb{E}\left[\sum_{t=1}^TX_{it}X_{it}'\right] \right)^{-1} \mathbb{E}\left[\sum_{t=1}^TX_{it}Y_{it}\right].$$ We estimate the parameters $$\theta_{10}=\frac{1}{N}\sum_{i=1}^N\eta_{i10},\quad \theta_{20}=\frac{1}{N}\sum_{i=1}^N\eta_{i10}^2,$$ which together are informative about the level and dispersion of race-based employment discrimination (here, the coefficient $\eta_{i10}$ corresponds to race). The true values in the DGP are $\theta_{10}=-0.0844$, which corresponds to a $8.44$ percentage point lower callback probability for blacks compared to whites, and $\theta_{20}=0.0177$, which corresponds to a standard deviation of $0.1030$.
We report results based on $1000$ simulations. We compute the OLS estimators as $$\widehat\eta_i=\left(\sum_{t=1}^TX_{it}X_{it}'\right)^{-1} \sum_{t=1}^TX_{it}Y_{it},$$ and report the resulting estimates of $\theta_1$ and $\theta_2$ (denoted as OLS). Note that, given the design, $\sum_{t=1}^TX_{it}X_{it}'$ is singular with positive probability. We report averages of $\widehat\eta_{i1}$ and $\widehat\eta_{i1}^2$ on the subset of units for which singularity does not occur, and compare those to the averages of $\eta_{i10}$ and $\eta_{i10}^2$ for the same subset of units. We proceed similarly for the other estimators described below.
Next, we compute orthogonal estimators to order 2 (denoted as ORTH). Let, for all $s_1,s_2\in \{1,...,T\}$, $$\widehat\eta_{i,-(s_1,s_2)}=\left(\sum_{t=1}^T\boldsymbol{1}\{t\neq s_1,s_2\}X_{it}X_{it}'\right)^{-1} \sum_{t=1}^T\boldsymbol{1}\{t\neq s_1,s_2\}X_{it}Y_{it},$$ and
For $\theta_1$, the orthogonal estimator is the second element (i.e., the one corresponding to $\eta_{i1}$) of
where $I_3$ is the $3\times 3$ identity matrix, averaged across units. For $\theta_2$, the orthogonal estimator is, for $D$ the $3\times 3$ matrix with a one in the position $(2,2)$ and zeroes everywhere else,
also averaged across units.
In addition to estimators based on the plug-in estimates ((ref)) of $\Lambda_i$, we also report estimators based on empirical Bayes regularization. Let $p_{ik\ell}=\Pr(X_{it1}=k,X_{it2}=\ell)$ for all $(k,\ell)\in\{0,1\}^2$. To regularize $\Lambda_i$, we endow $(p_{i00},p_{i01},p_{i10},p_{i11})$ with a Dirichlet prior with parameters $(\alpha_{00},\alpha_{01},\alpha_{10},\alpha_{11})$. Let $\alpha=\alpha_{00}+\alpha_{01}+\alpha_{10}+\alpha_{11}$, $\pi_{k\ell}=\alpha_{k\ell}/\alpha$, and $$\Pi=\left(
\right).$$ We estimate $-\Lambda_i^{-1}=\mathbb{E}[X_{it}X_{it}']$ as the posterior mean $$-\widehat \Lambda_i^{-1}=\frac{T}{T+\alpha}\left(\frac{1}{T}\sum_{t=1}^TX_{it}X_{it}'\right)+\frac{\alpha}{T+\alpha}\Pi.$$
Finally, to select $(\alpha_{00},\alpha_{01},\alpha_{10},\alpha_{11})$, we maximize the marginal likelihood across firms with respect to $\alpha$, while fixing $\widehat{\pi}_{k\ell}=\frac{1}{NT}\sum_i\sum_t \boldsymbol{1}\{X_{it1}=k,X_{it2}=\ell\}$ to the pooled mean estimator. We estimate $\alpha$ as $$\widehat{\alpha}=\underset{\alpha>0}{\mbox{argmax}}\, N[\ln\Gamma (\alpha)-\ln\Gamma (\alpha+T)]+\sum_{i=1}^N\sum_{(k,\ell)\in\{0,1\}^2}\left[\ln\Gamma (\widehat\pi_{k\ell}\alpha+n_{ik\ell})-\ln\Gamma (\widehat\pi_{k\ell}\alpha)\right],$$ where $n_{ik\ell}=\sum_{t=1}^T\boldsymbol{1}\{X_{it1}=k,X_{it2}=\ell\}$. We apply the regularization scheme to both the OLS and ORTH estimates. In the latter case, we apply the above formulas to subsets of observations, leaving out two $t$ indices and averaging across pairs of indices ex post.
Figures (ref) and (ref) report the Monte Carlo performance of the OLS (in blue) and ORTH (in red) estimators for $\theta_1$ and $\theta_2$, respectively, as the number of applications per firm $T$ varies between $T=20$ and $T=100$. Each figure plots the mean estimate across simulations (solid line), the true value (dashed horizontal line), and 90% simulation bands defined by the 5th and 95th quantiles across simulations (shaded regions). The plots on the left-hand side show the estimators based on regularized $\Lambda_i$'s, while the plots on the right-hand side correspond to un-regularized $\Lambda_i$'s.
Starting with $\theta_{10}$ -- the average of $\eta_{i10}$ -- the right graph in Figure (ref) shows that both the OLS and the second-order orthogonalized estimator are close to unbiased irrespective of $T$. The left graph, in turn, shows that while OLS is slightly biased as a result of the regularization of $\Lambda_{i0}$, the orthogonal estimator remains virtually unbiased. In addition, the confidence bands show that the orthogonalization is associated to some increase in sampling variability.
Shifting attention to $\theta_{20}$ -- the average of $\eta_{i10}^2$ -- the right graph in Figure (ref) shows that both estimators are biased when one does not rely on a regularized estimator of $\Lambda_{i0}$. Moreover, the bias is higher for the orthogonal estimator, and it only becomes lower than the bias of OLS when $T\geq 60$. Lastly, the variance increases associated with the orthogonalization is larger than in the case of $\theta_{10}$, especially for lower values of $T$. Hence, while orthogonalization reduces bias for large enough $T$, an estimator based on un-regularized estimates of $\Lambda_{i0}$ may perform poorly for lower values of $T$.
The left graph in Figure (ref) shows that the situation improves greatly when relying on regularized estimates of $\Lambda_{i0}$. The bias of OLS is reduced relative to the un-regularized case, without a noticeable increase in variability. Moreover, the bias of the orthogonal estimator decreases substantially relative to the un-regularized case, and in fact the orthogonal estimator has small bias, lower than the one of OLS for all values of $T$. This suggests that the combination of orthogonal moments -- to reap the benefits of higher-order Neyman orthogonality -- and regularization of the nuisance parameters -- to ensure lower estimation error -- can be particularly appealing in applications.
In this section we establish the asymptotic distribution of the estimator based on a $q$-th order Neyman-orthogonal moment function. The assumptions and result cover grouped data applications in which each unit contributes a within-unit U-statistic over within-group observations, including the heterogeneous-coefficients application of Sections (ref) and (ref) as a special case. We state the result for a $d_\theta$-dimensional target parameter $\theta_0$, with an orthogonal moment function $\psi$ taking values in $\mathbb{R}^{d_\psi}$ for some $d_\psi \geq d_\theta$. Under the assumptions below, the estimator is $\sqrt{N}$-consistent and asymptotically normal, where $N$ denotes the number of units (groups), even when the preliminary nuisance estimates converge at a rate slower than $\sqrt{N}$.
The analysis uses the unit-by-unit view introduced in Section (ref). For each unit $i\in\{1,\ldots,N\}$, the $T$ observations $W_{i1}, \ldots, W_{iT}$ form the sample on which the orthogonal moment function $\psi(W_i; \theta, \nu_i)$ is constructed, and preliminary estimation of the unit-specific nuisance $\nu_{i0} = (\eta_{i0}, \lambda_{i0})$ is carried out on a held-out subset of observations for the same unit. The cross-sectional dimension $i = 1, \ldots, N$ enters only at the final step, through the sample average that defines the estimator of $\theta_0$. As a consequence, the $\sqrt{N}$-rate and the asymptotic normality in Theorem (ref) below are driven by a cross-sectional central limit theorem. The within-unit structure of $\psi$ affects only the asymptotic variance.
The analysis closely parallels Section 6 of BonhommeJochmansWeidner2025. The main adaptation is that the nuisance parameter $\eta_i$ in that paper is replaced by the pair $\nu_i = (\eta_i, \lambda_i)$, and that the moment function is our orthogonal U-statistic rather than a conditional-likelihood score. Because the limiting distribution is driven by the cross-sectional average, with the within-unit structure only entering through the variance of the moment function, the proof of BonhommeJochmansWeidner2025 applies after a straightforward notational translation. We record this translation in Appendix (ref) rather than reproducing the full proof.
We observe an independent cross-section of $N$ units. Each unit contributes a bundle $W_i = (W_{i1}, \ldots, W_{iT})$ of $T$ observations. Both $N$ and $T$ grow as $N \to \infty$; we suppress the dependence of $T$ on $N$ in the notation. For each unit $i$, let $\nu_{i0} = (\eta_{i0}, \lambda_{i0})$ denote the true value of the combined nuisance parameter, where $\lambda_{i0}$ is the finite-dimensional parameter encoding the unit-level Jacobian inverse $\Lambda_{i0}$ as in Section (ref). The dimension $d_{\nu,i}$ is bounded uniformly in $i$ and $N$.
The target parameter $\theta_0 \in \mathbb{R}^{d_\theta}$ is identified by
where $\psi(W_i; \theta, \nu_i) \in \mathbb{R}^{d_\psi}$ is a $q$-th order Neyman-orthogonal moment function with respect to $\nu_i$ in the sense of Definition (ref). We allow $d_\psi \geq d_\theta$ to permit overidentification. In applications such as the heterogeneous coefficients example of Section (ref), $\theta_0 = (1/N) \sum_i m(\eta_{i0})$ is a sample-dependent parameter; the analysis below covers this case, with the understanding that $\theta_0$ may depend on $N$.
Let $\widehat\nu_i$ be a preliminary estimator of $\nu_{i0}$. We assume throughout that $\widehat\nu_i$ is constructed from a held-out subset of observations, independent of the observations used to evaluate $\psi$. Concretely, for each unit $i$ we partition the within-unit index $\{1, \ldots, T\}$ into disjoint subsets $\mathcal{S}_1$ and $\mathcal{S}_2$, estimate $\widehat\nu_i$ from $\{W_{it} : t \in \mathcal{S}_1\}$, and evaluate $\psi$ using observations in $\mathcal{S}_2$. The sample is split within each unit, not across units. Efficiency can be improved by cross-fitting, by swapping the roles of $\mathcal{S}_1$ and $\mathcal{S}_2$ and averaging the resulting moment functions. The analysis below applies to either a single split or the cross-fitted estimator.
The estimator $\widehat\theta$ is defined by the GMM problem
where $\Theta \subseteq \mathbb{R}^{d_\theta}$ is the parameter space, $\Omega$ is a symmetric positive-definite $d_\psi \times d_\psi$ weight matrix, and $\|x\|_\Omega^2 = x'\Omega x$. When $d_\psi = d_\theta$ (just-identification), $\widehat\theta$ solves the moment equation $(1/N) \sum_i \psi(W_i; \widehat\theta, \widehat\nu_i) = 0$.
Part (ii) mirrors the bounded-nuisance-dimension condition of BonhommeJochmansWeidner2025. In the grouped data example of Section (ref), $d_{\nu,i} = d_\eta + d_\eta(d_\eta +1)/2$, which is independent of $N$ and $T$. Part (iv) is the key rate condition on the preliminary estimator. In the grouped data example, $\widehat\nu_i$ is an OLS-based estimator computed on $|\mathcal{S}_1|$ observations, so $\mathbb{E}[\|\widehat\nu_i - \nu_{i0}\|^{2(q+1)}] = O(T^{-(q+1)})$ if $|\mathcal{S}_1|/T$ tends to a non-zero constant, and the rate condition is satisfied whenever $N = o(T^{q+1})$. Higher-order orthogonality thus permits inference in grouped settings where the group sizes are small relative to the number of groups. Part (v) is a standard rank condition.
Part (ii) is achieved by sample splitting within each unit, as described above. Part (iii) requires independence across units. It is straightforward to modify the variance formula in Theorem (ref) below to accommodate particular forms of dependence across units, by substituting an appropriate expression for $V$.
The proof is in Appendix (ref).
The asymptotic variance depends on the orthogonality order $q$ through the form of the asymptotic variance, although we have left this dependence implicit in the notation. As in BonhommeJochmansWeidner2025, the use of cross-fitting mitigates the variability introduced by a single split.