Debiased Machine Learning: Identification, Estimation, and Shape Constraints

Debiased Machine Learning: Identification, Estimation, and Shape Constraints

\pdfbookmark[1]{Title}{title}

Introduction

Debiased machine learning (DML) has emerged as a powerful framework for statistical inference in empirical research, particularly in applications with big data AhrensChernozhukovHansenKozburSchafferWiemann2026DML. In these settings, the parameter of interest $\theta_0$ is identified by a moment restriction that depends on a first step nuisance parameter $\gamma_0$. Thus, inference on $\theta_0$ must account for estimation error in $\gamma_0$, a challenge exacerbated when $\gamma_0$ is high dimensional. While machine learning is well suited for estimating such high dimensional objects, it tends to induce regularization and overfitting biases that may lead to biased and $\sqrt{n}$-inconsistent estimation of $\theta_0$ ChernozhukovChetverikovDemirerDufloHansenNeweyRobins2018Double,ChernozhukovEscancianoIchimuraNewey2022LocalRobust. DML alleviates these biases by orthogonalizing the moment with respect to $\gamma_0$, coupled with cross-fitting (a form of sample splitting). At the core of the general DML framework in ChernozhukovEscancianoIchimuraNewey2022LocalRobust,ChernozhukovNeweySingh2022Automatic is an additional nuisance parameter $\alpha_0$, known as the Riesz representer, which is in general necessary to achieve orthogonality through the so-called first step influence function. As argued in the literature, since $\alpha_0$ may lack a closed-form expression or involve inverting unknowns, it is important to estimate $\alpha_0$ in an automatic way without the knowledge of its analytic form.

Despite several existing strategies for the automatic estimation of $\alpha_0$ ChernozhukovNeweyQuintasSyrgkanis2024RieszReg, the identification of $\alpha_0$, to the best of our knowledge, remains undeveloped at a general level, particularly for first steps defined by models with endogeneity. As our first main contribution, we establish conditions under which $\alpha_0$ is identified by the moment conditions underlying the automatic estimation procedures in ChernozhukovEscancianoIchimuraNewey2022LocalRobust,ChernozhukovNeweySingh2022Automatic. The conditions we impose are mild in the sense that they consist of standard smoothness requirements, a linearity condition with respect to $\alpha_0$ that reflects a common feature of the first step influence function, and a {\it coercivity} condition that is not only sufficient but also necessary in a certain sense.

While the identifying equation may in principle be exploited to directly estimate $\alpha_0$, it is important to obtain a tractable extremal reformulation, which is then amenable to penalized optimization as common in machine learning. Our second contribution shows that $\alpha_0$ is identified precisely when it uniquely optimizes a quadratic functional $\mathfrak C_0$. Crucially, the optimization only depends on the knowledge of the functional forms of the original moment restriction and the first step influence function but not the analytic expression of $\alpha_0$, thereby enabling automatic estimation. Our result also reveals that the quadratic characterization recurring in the literature (see, e.g., ChernozhukovEscancianoIchimuraNewey2022LocalRobust,ChernozhukovNeweySingh2022Automatic,ChernozhukovNeweySinghSyrgkanis2024Adversarial,ChernozhukovNeweyQuintasSyrgkanis2024RieszReg and Singh2024KRRR) arises from a more general level.

As our third contribution, we develop an automatic estimation procedure of $\alpha_0$ based on the empirical analog of $\mathfrak C_0$, known as Riesz regression, which generalizes ChernozhukovEscancianoIchimuraNewey2022LocalRobust,ChernozhukovNeweySingh2022Automatic,ChernozhukovNeweyQuintasSyrgkanis2024RieszReg to settings where $\gamma_0$ may be defined by models with endogeneity, e.g., nonparametric instrumental variable (NPIV) models. The generalization appears nontrivial and poses nontrivial challenges because the parameter space of $\alpha_0$ may be unknown. We derive the convergence rates of the resulting estimator under conditions that allow for penalization and nonlinear sieve approximations (e.g., neural networks). At the technical level, a key ingredient in characterizing the rates is the modulus of continuity of an empirical process, which connects to the critical radius in statistical learning Wainwright2019HighStats,ChernozhukovNeweySinghSyrgkanis2024Adversarial.

Our final contribution aims to integrate shape constraints into DML. While machine learning is powerful for processing big data, it is no panacea in delivering statistical guarantees and equally suffers from the curse of dimensionality dictated by the optimal convergence rates Stone1982RateGlobal,HallHorowitz2005NPIV,ChenReiss2011Rate. Therefore, in high dimensional settings, additional structures must be exploited to ensure that the convergence rates are fast enough. Following Stone1985Additive,Stone1994Spline, some recent studies in deep learning exploit “low dimensional” structures such as hierarchical interaction Schmidt2020ReLU,KohlerLanger2021RateDNN. These structures, however, may not always be easy to motivate. Shape constraints, on the other hand, are deeply grounded in economics because extensive economic theories are formulated as shape constraints Matzkin1994Handbook,ChetverikovSantosAzeem2018Shape. Thus, they may be viewed as alternative structures that help improve estimation precision. In developing our general theory, we operate under the restriction $\gamma_0\in\Gamma$ for some possibly nonlinear parameter space $\Gamma$ which in turn restricts $\alpha_0$. Our identification, characterization, and estimation results for $\alpha_0$ thus allow for general nonlinear shape constraints, though we stress that these results appear novel even without shape constraints.

We complement our contributions on Riesz regression and shape constraints by simulation studies and empirical applications. We focus on implementation via deep learning, not only to maintain a coherent treatment but also because neural networks provide a relatively flexible and tractable way to incorporate shape constraints. Our simulation designs include endogenous as well as exogenous first steps. In both cases, incorporating shape constraints overall improves the performance of point estimates and confidence intervals, with larger gains as more constraints are imposed. Our first application is a difference-in-differences (DiD) analysis of Medicaid expansions and mortality, where monotonicity is imposed based on institutional knowledge well documented in the literature. The second application studies wage growth and working hours, imposing convexity as an implication of economic theory.

This paper contributes to the extensive and rapidly evolving literature on debiased machine learning. Building on the semiparametric literature Newey1994AsymptoticVar,RobinsRotnitzkyZhao1994Regression,IchimuraNewey2022IF, ChernozhukovEscancianoIchimuraNewey2022LocalRobust,ChernozhukovNeweySingh2022Automatic developed a general framework for automatic DML while leaving the identification of $\alpha_0$ for future work. We thus complement their work by laying the identification foundation. The automatic estimation of $\alpha_0$ in ChernozhukovEscancianoIchimuraNewey2022LocalRobust,ChernozhukovNeweySingh2022Automatic is based on a LASSO minimum distance criterion where $\gamma_0$ is defined by a “generalized linear regression” without endogeneity and subject to linear constraints. ChernozhukovNeweyQuintasSyrgkanis2024RieszReg obtained the quadratic characterization of $\alpha_0$ for the same class of first steps and developed the Riesz regression that facilitates the use of nonlinear machine learners---see also ChernozhukovNeweySingh2022GlobalLocalDML, ChernozhukovNeweySinghSyrgkanis2024Adversarial, and Singh2024KRRR for constructions in related settings. Our quadratic characterization is more general in that $\theta_0$ may be defined by nonlinear moment conditions, $\gamma_0$ may be defined by models with endogeneity, and the parameter space of $\gamma_0$ may be nonlinear.

There is relatively less work on automatic DML for endogenous first steps, even under linear constraints. While the setup in ChernozhukovEscancianoIchimuraNewey2022LocalRobust does accommodate endogenous $\gamma_0$, they primarily focus on exogenous $\gamma_0$ after showing how first step influence functions can be used to construct orthogonal moment conditions. We build on their work by taking the orthogonality as given. Bakhitov2026PGMM developed a LASSO minimum distance estimator of $\alpha_0$ for functionals of the NPIV regression. Existing simulation evidence, however, suggests that allowing for nonlinear machine learners can materially affect inference; see, e.g., ChernozhukovChetverikovDemirerDufloHansenNeweyRobins2018Double and AhrensChernozhukovHansenKozburSchafferWiemann2026DML. BennettKallusMaoNeweySyrgkanisUehara2025StrongID studied linear functionals of first steps that may be partially identified by linear conditional moment restrictions under a “strong identification” condition. When specialized to this setup, our general theory implies that the strong identification condition is stronger than necessary for the identification of $\alpha_0$---see also Appendix (ref) for more comparisons. Bruns2025TSML developed a two-stage machine learning estimator of the NPIV regression by adaptively learning the instrument basis. Our identification, characterization, and estimation results appear novel and more general while also allowing for nonlinear constraints on $\gamma_0$.

Our work is also related to the literature on shape constraints which have long played prominent roles in economics Marshall1890Principles. Statistical studies on shape constraints emerged in the 1950s Hildreth1954Concave,Grenander1956II,BrunkEwingUtz1957Minimize and have since developed into a broad and influential area of research---see Matzkin1994Handbook, ChetverikovSantosAzeem2018Shape, and SamworthSen2018Shape for reviews. Our use of shape restrictions is driven by their potential in improving estimation BlundellHorowitzParey2012Measuring,BlundellHorowitzParey2017Slutsky,ChetverikovWlihelm2017MonNPIV,HorowitzLee2017Shape. Shape constraints may also be viewed as another form of regularization, as emphasized in ChetverikovWlihelm2017MonNPIV who study NPIV models under monotonicity. This perspective is particularly helpful in deep learning which often depends on various forms of hyperparameters and regularization and yet a general theory on data-driven tuning is currently unavailable. We develop a general framework of identification and estimation that brings shape constraints into DML and our simulation results show that shape constraints may improve estimation and inference even with crude tuning.

Finally, we introduce some notation and concepts. For $a,b\in\mathbf R$, set $a\vee b\equiv\max\{a,b\}$ and let $a\lesssim b$ mean $a\le cb$ for some constant $c>0$. For a set $A$ in a vector space $\mathbf H$ with seminorm $\|\cdot\|_{\mathbf H}$,\footnote{A seminorm $\|\cdot\|_{\mathbf H}\colon \mathbf H\to\mathbf R$ is such that $\|a+b\|_{\mathbf H}\le \|a\|_{\mathbf H}+\|b\|_{\mathbf H}$ and $\|ta\|_{\mathbf H}=|t|\|a\|_{\mathbf H}$ for all $t\in\mathbf R$ and $a,b\in\mathbf H$. If in addition $a=0$ whenever $\|a\|_{\mathbf H}=0$, then $\|\cdot\|_{\mathbf H}$ is a norm.} its interior $A^\circ$ is the largest open set contained in $A$, its convex conical hull $\mathrm{con}(A)$ is the smallest convex cone in $\mathbf H$ containing $A$, its linear span $\mathrm{lin}(A)$ is the smallest subspace containing $A$, and its closed linear span $\overline{\mathrm{lin}}(A)$ is the closure of $\mathrm{lin}(A)$. The distance from $a\in\mathbf H$ to $A$ is $d_{\mathbf H}(a,A)\equiv \inf_{a'\in A}\|a-a'\|_{\mathbf H}$. A map $\phi\colon \Gamma\subset\mathbf H\to\mathbf R$ is Gateaux differentiable at $\gamma_0\in\Gamma$ tangentially to a set $H_0\subset\mathbf H$ if there is a linear map $\nabla_\gamma\phi(\gamma_0)\colon\mathrm{lin}(H_0)\to\mathbf R$ such that $\lim_{t\downarrow 0}\{\phi(\gamma_0+th)-\phi(\gamma_0)\}/t=\nabla_\gamma\phi(\gamma_0)[h]$ whenever $h\in H_0$ and $\gamma_0+th\in\Gamma$ for all small $t\ge 0$. For $X\in\mathcal X$ with law $P$, let $L^2(X)\equiv\{f\colon\mathcal X\to\mathbf R\colon \|f\|_{P,2}<\infty\}$ with $\|f\|_{P,2}\equiv\{E[|f(X)|^2]\}^{1/2}$ and $L^\infty(X)\equiv\{f\colon\mathcal X\to\mathbf R\colon \|f\|_{P,\infty}<\infty\}$ with $\|f\|_{P,\infty}\equiv\inf\{a\ge 0\colon P(|f(X)|>a)=0\}$. For $\{\phi(\cdot,\gamma)\colon\gamma\in\Gamma\}\subset L^2(X)$, we say that $\phi(X,\cdot)$ is mean-square continuous at $\gamma_0$ if the map $\gamma\mapsto\phi(\cdot,\gamma)\in L^2(X)$ is continuous at $\gamma_0$.

The remainder of the paper is organized as follows. Section (ref) introduces the setup and related examples. Section (ref) develops our general framework of automatic DML under shape constraints. Section (ref) discusses implementation details and presents simulation studies, while Section (ref) showcases empirical applications. Section (ref) concludes. All proofs and supporting results are relegated to Appendices (ref), (ref), and (ref).

The Setup and Examples

As in the lirature Newey1994AsymptoticVar,ChernozhukovEscancianoIchimuraNewey2022LocalRobust, we work with a generalized method of moments (GMM) model that accommodates numerous causal and structural parameters. Let $X\in\mathcal X$ be a vector of observables with $\mathcal X$ a sample space and $\theta_0\in\Theta\subset\mathbf R^{d_\theta}$ be a parameter of interest such that

admits a unique solution at $\theta=\theta_0$, where $g\colon \mathcal X\times \Theta\times\Gamma^\dag\to\mathbf R^{d_g}$ is a known map with $d_g\ge d_\theta$ and $\gamma_0\in\Gamma^\dag\subset \mathbf H$ is a nuisance parameter in a vector space $\mathbf H$. The generality of our setup also stems from the flexibility in specifying $\gamma_0$ and its parameter space $\mathbf H$. First, we do not take a stand on the nature of $\gamma_0$ as it may be parametric or nonparametric, causal or structural, or a vector of nuisance parameters. Second, $\mathbf H$ may be equipped with a strong or weak topology depending on the application. There are settings where it is more appropriate to work with a (weak) seminorm than a (strong) norm, and vice versa---see Remark (ref) for details.

Our general theory extends the existing literature ChernozhukovEscancianoIchimuraNewey2022LocalRobust,ChernozhukovNeweySingh2022Automatic,ChernozhukovNeweyQuintasSyrgkanis2024RieszReg by accommodating two important additional features of the setup in (ref). First, the first step $\gamma_0$ entering the nonlinear moment conditions (ref) may be defined by a model with endogeneity; e.g., $\gamma_0$ may be an NPIV regression.Such endogeneity poses substantial challenges for identifying and estimating the debiasing nuisance $\alpha_0$, because the relevant parameter space for $\alpha_0$ is then unknown. Second, $\gamma_0$ may be subject to nonlinear shape constraints, which is encoded by restricting $\gamma_0\in\Gamma\subset\Gamma^\dag$. For example, if $\gamma_0$ is monotone, then $\Gamma$ may be taken to be the set of monotone functions in $\Gamma^\dag$. These shape constraints, in turn, induce corresponding restrictions on $\alpha_0$, as our theory below makes clear. Although a few studies allow for endogenous first steps in special settings (see the introduction), they focus on estimation. By contrast, we study not only the estimation of $\alpha_0$ but also its identification and characterization, and we do so in a more general framework in the sense that $\theta_0$ is defined implicitly by a general GMM model (ref) and $\gamma_0$ is a generic first step subject to possibly nonlinear shape constraints. To our knowledge, the incorporation of general nonlinear shape constraints appears to be new to the DML literature. We also note that our identification result below for $\alpha_0$ appears novel even for exogenous first steps.

We follow the blueprint for automatic DML laid out in ChernozhukovEscancianoIchimuraNewey2022LocalRobust, which is based on the notion of a first step influence function. Intuitively, this function captures the effect of first step estimation errors on the moment condition in (ref). Since it plays a fundamental role in our analysis, we next provide a brief review.

The First Step Influence Function

First step estimation errors may be viewed as finite sample counterparts of local perturbations in $\gamma_0$ at the population level. In order to describe these perturbations, we assume that $X$ is distributed according to $P\in\mathbf P$ where $\mathbf P$ is a collection of distributions that possibly generate $X$ while identifying $\theta_0$ via the model (ref). Let $P_t=(1-t)P+tQ$ for a distribution $Q$ such that $P_t\in\mathbf P$ and the associated first step, denoted $\gamma (P_t)$, belongs to $\Gamma^\dag$ for each $t\in[0,1]$. Then the first step influence function is a function $\phi\colon\mathcal X\times\Theta\times\Gamma^\dag\times\mathbf A\to\mathbf R^{d_g}$ with $\mathbf A$ a suitable space that satisfies, for some $\alpha_0\in\mathbf A$ and all $Q$ such that each $P_t\in\mathbf P$ and $\gamma (P_t)\in\Gamma^\dag$, the restrictions (i) all entries of $\phi(\cdot,\theta_0,\gamma_0,\alpha_0)$ are in $L^2(X)$ and (ii)

where $|_{t=0+}$ indicates the right derivative at $0$. Note that (ref) implies $E[\phi(X,\theta_0,\gamma_0,\alpha_0)]=0$. The presence of an additional nuisance $\alpha_0$ is a common feature in DML.

The parameter $\alpha_0$ plays a key role in DML because it helps orthogonalize the influences of the nuisance parameters through the first step influence function, without compromising the identification of $\theta_0$. Specifically, the modified moment condition

continues to uniquely identify $\theta_0$ and is insensitive to small changes in both $\gamma_0$ and $\alpha_0$. In particular, ChernozhukovEscancianoIchimuraNewey2022LocalRobust establish under regularity conditions the following local orthogonality condition with respect to $\gamma_0$:

with $\delta$ ranging over $\mathbf H$ such that $\gamma_0+t\delta\in\Gamma^\dag$. A stronger orthogonality condition holds for $\alpha_0$ but it is (ref) that matters for our discussions below.

We do not aim to derive first step influence functions but instead take them as given by leveraging existing work; see, e.g., ChernozhukovEscancianoIchimuraNewey2022LocalRobust, IchimuraNewey2022IF, and references therein. We stress that, in line with ChernozhukovEscancianoIchimuraNewey2022LocalRobust, we view first step influence functions as means of correcting first step biases rather than as a route to deriving efficient influence functions. Indeed, the asymptotic theory in Section (ref) merely relies on the orthogonality---see the discussions of Assumption (ref). A common structure of first step influence functions is that they are linear in $\alpha_0$, which we illustrate through examples below.

Examples

To fix ideas, we introduce examples that play predominant roles in empirical work. Importantly, the first step influence functions in these settings are all linear in $\alpha$. While we focus on these examples for conciseness, we note that such linearity also appears in related settings such as moment conditions with density as first steps Newey1994AsymptoticVar, dynamic discrete choice models ChernozhukovEscancianoIchimuraNewey2022LocalRobust, inference on support functions in a class of partially identified models Semenova2023SetDML, and models with generated regressors EscancianoPerez2025DML.

The first example concerns a parameter $\theta_0$ that depends on a first step $\gamma_0$ defined by exogenous orthogonal conditions.

Our second example is a special case of Example (ref), which we single out due to its importance in the study of welfare analysis HausmanNewey2017Welfare.

Our third example generalizes Example (ref) to settings where the first step is defined by endogenous orthogonality conditions.

The General Framework

We develop our general theory in three steps. In Section (ref), we establish the identification of $\alpha_0$ when the nuisance parameter $\gamma_0$ is subject to possibly nonlinear shape constraints in addition to those already incorporated in $\Gamma^\dag$, and then characterize $\alpha_0$ as the unique solution to a quadratic optimization problem. In Section (ref), we study automatic estimation of $\alpha_0$ based on the quadratic characterization. In Section (ref), we obtain the asymptotic distribution of the debiased GMM estimator $\hat\theta_n$.

Identification and Characterization

The parameter $\alpha_0$ originates from deriving the first step influence function in (ref). As mentioned in the introduction, $\alpha_0$ may not admit a closed-form expression, especially under shape constraints on $\gamma_0$. When it does have a closed-form expression, it may involve inverting unknown objects (e.g., densities) that lead to estimators with poor finite sample performance. These concerns motivate automatic estimators, which tend to perform better, as shown in, e.g., ChernozhukovNeweyQuintasSyrgkanis2022RieszNet and AhrensChernozhukovHansenKozburSchafferWiemann2026DML. Thus, it is important to estimate $\alpha_0$ in an automatic way without requiring knowledge of its analytical form. A general strategy for the automatic estimation of $\alpha_0$ is proposed by ChernozhukovEscancianoIchimuraNewey2022LocalRobust who developed a LASSO minimum distance estimator of $\alpha_0$ based on sample analogs of the moment conditions in (ref). This, however, raises an important question: Does the system (ref) admit a unique solution?

To address this question, we must choose the set of feasible perturbations that reflect the shape restriction(s) incorporated in $\Gamma$. Our theory in Section (ref) suggests requiring (ref) to hold for all $\delta\in \Gamma-\gamma_0\equiv\{\gamma-\gamma_0\colon \gamma\in\Gamma\}$---see discussions of Assumption (ref)(ii) in Section (ref). To this end, we assume throughout that $\Gamma^\dag$ is linear or at least convex so that $\gamma_0+t\delta\in\Gamma^\dag$ for all $t\in[0,1]$ and $\delta\in \Gamma-\gamma_0$, though no such requirements are imposed on $\Gamma$. Since $\alpha_0$ is defined for each moment equation, below we focus on a single moment equation by pretending that $g$ is real-valued (i.e., $d_g=1$). To formalize our discussions, we introduce the first assumption.

Assumptions (ref)(i)(iii) simply impose minimal smoothness on the original moment function and the first step influence function so that (ref) is well-defined, though we note that Fr\'{e}chet differentiability is required when establishing the asymptotic distribution of the debiased GMM estimator $\hat\theta_n$. Assumptions (ref)(ii)(iv) demand that the resulting derivatives possess certain composition structures. In particular, $\Pi_0$ serves as the link when $\gamma_0$ and $\alpha_0$ are functions of different variables. This is motivated by a necessary condition for $\sqrt n$-estimability of $\theta_0$---see SeveriniTripathi2012Efficiency and also IchimuraNewey2022IF. In settings such as Examples (ref) and (ref) where $\gamma_0$ and $\alpha_0$ are functions of the same variable, one may simply set $\Pi_0$ to be the identity map. The linearity of $\Upsilon_0$ with respect to $\alpha$ in Assumption (ref)(iv) is inherited from the same property of first step influence functions in Section (ref).

Our first main result below formally establishes under Assumptions (ref) that (ref) admits a unique solution $\alpha_0$ under an additional condition on $\Upsilon_0$. Following ChernozhukovNeweySingh2022Automatic, we refer to this unique solution as the Riesz representer.

Theorem (ref) holds regardless of whether the parameter space $\Gamma$ is linear or not. Note that $\alpha_0$ is subject to the shape constraint on $\gamma_0$ via $\alpha_0\in\mathbf A_0$ (the closure of the range $\Pi_0(\mathbf H_0)$). To appreciate the coercivity condition, suppose that $\mathbf A_0$ is finite dimensional so that we may identify $\Upsilon_0$ with some matrix $\Phi_0$. Then (ref) implies but is not implied by the invertibility of $\Phi_0$. In a certain sense, coercivity is not only sufficient but also necessary---see Remark (ref). Theorem (ref) is a consequence of the Lax-Milgram theorem LaxMilgram1954Parabolic, which is a foundational tool in partial differential equations. On a technical note, while the Lax-Milgram theorem may be viewed as a generalization of the Riesz representation theorem, the latter is not directly applicable by defining an inner product through $\Upsilon_0$, unless $\Upsilon_0$ is symmetric and has a definite sign (i.e., positive or negative).\footnote{Recall that $\Upsilon_0\colon \mathbf A_0\times\mathbf A_0\to\mathbf R$ is symmetric if $\Upsilon_0(a,b)=\Upsilon_0(b,a)$ for all $a,b\in\mathbf A_0$, positive if $\Upsilon_0(\alpha,\alpha)\ge 0$ for all $\alpha\in\mathbf A_0$, and negative if $\Upsilon_0(\alpha,\alpha)\le 0$ for all $\alpha\in\mathbf A_0$.} These additional properties, however, allow us to obtain an extremal characterization of $\alpha_0$, which we turn to next.

Theorem (ref) suggests that we may estimate $\alpha_0$ directly based on the equation (ref). However, it is desirable to characterize $\alpha_0$ in a way that facilitates regularization which is critical in machine learning. Our next theorem relates $\alpha_0$ to an optimization problem, thereby enabling regularization via penalization.

Theorem (ref) states that the Riesz representer is identified by (ref) as $\delta$ ranges over $\Gamma-\gamma_0$ precisely when it is the unique solution over $\mathbf A_0$ that optimizes the quadratic functional $\mathfrak C_0$. Hence, automatic estimation of $\alpha_0$ based on optimizing $\mathfrak C_0$ over $\mathbf A_0$ is possible because this problem only depends on the moment function (through $\Upsilon_0$, $\Psi_0$ and $\Pi_0$). ChernozhukovNeweyQuintasSyrgkanis2024RieszReg obtained essentially the same characterization when $\theta_0$ is of the form $E[m(X,\gamma_0)]$ for $\gamma_0$ defined by a generalized regression without endogeneity under linear constraints. Finally, symmetry and positiveness/negativeness of $\Upsilon_0$ play different roles in Theorem (ref)---see Remark (ref).\footnote{We thank Whitney K.\ Newey for suggesting that we clarify the role played by symmetry.}

Examples Revisited

In this section, we revisit the examples in Section (ref) by verifying the assumptions in Theorem (ref). We omit Example (ref) as it is a special case of Example (ref).

Riesz Regression

In this section, we develop a general procedure for estimating the Riesz representer $\alpha_0$ based on Theorem (ref). As in the literature, the estimation is automatic in the sense that only the knowledge of the functional forms of the moment function and the first step influence function are required but not the analytic expression of $\alpha_0$ itself. We call the procedure a Riesz regression following ChernozhukovNeweyQuintasSyrgkanis2024RieszReg, which is a constrained quadratic optimization problem.

As a first step, we approximate the problem by a “parametrized” version. Despite the notation, the linear span $\mathbf H_0\equiv\mathrm{lin}(\Gamma-\gamma_0)$ does not involve $\gamma_0$ because $\mathrm{lin}(\Gamma-\gamma_0)=\mathrm{lin}(\Gamma-\Gamma)$, which further simplifies to $\mathrm{con}(\Gamma)-\mathrm{con}(\Gamma)$ if $0\in\Gamma$, to $\Gamma-\Gamma$ if $\Gamma$ is a convex cone, and to $\Gamma$ if $\Gamma$ is linear; see Lemma (ref). For $\Delta\Gamma_n\subset \mathrm{lin}(\Gamma-\Gamma)$ a sieve space of $\mathrm{lin}(\Gamma-\Gamma)$, we then approximate $\mathbf A_0$ by $\mathbf A_{0,n}\equiv\Pi_0(\Delta\Gamma_n)$. In Section (ref), we shall provide details of concrete constructions.

In order to obtain fruitful results, we next formalize common structures of the quadratic functional $\mathfrak C_0$ by imposing the following assumption. As in Section (ref), we keep pretending $d_g=1$ in what follows.

Assumption (ref) essentially reformulates conditions in Theorem (ref) when the maps $\Upsilon_0$ and $\Psi_0$ take the forms of expectations, with two caveats. First, Assumption (ref)(iii) requires $\Upsilon_0$ to be negative, but this is without loss of generality and one may instead consider positive $\Upsilon_0$. Second, the conditions on $\Upsilon_0$ and $\Psi_0$ are strengthened to hold over $\mathbf A$ instead of $\mathbf A_0$. This is not necessary when $\Pi_0$ is known but is needed when $\Pi_0$ and hence $\mathbf A_0$ are unknown and therefore must be estimated.

To estimate $\alpha_0$, we first consider the case when $\Pi_0$ is known (e.g., the identity map). Given Assumptions (ref) and a sample $\{X_i\}_{i=1}^n$ of $X$, we define:

where $\hat\nu_n\in\mathbf B_n\subset\mathbf B$, and let $\hat\alpha_n\in\mathbf A_{0,n}$ satisfy

where $\delta_n$ is the order of optimization error, and $J_n\colon\mathbf A_{0,n}\to\mathbf R_+$ is a penalty function with the amount of penalization dictated by the tuning parameter $\lambda_n\ge 0$. One may set $\lambda_n=0$, resulting in an estimator without penalization. We also note that, in certain settings such as NPIV models, the linear term $\Psi_0$ in $\mathfrak C_0$ may be automatically estimated without needing to estimate $\nu_0$---see Remark (ref).

Establishing statistical guarantees of $\hat\alpha_n$ requires restrictions on the sieve, the penalization, and the sample. For ease of notation, let $\mathbb G_nf\equiv \sum_{i=1}^{n}\{f(X_i)-E[f(X)]\}/\sqrt n$ for a generic function $f$ and $f(X,\alpha,\nu)\equiv l(X,\alpha,\alpha,\nu)/2+h(X,\alpha,\nu)$ so that $\mathbb G_nf(\alpha,\nu)=\sum_{i=1}^{n}\{f(X_i,\alpha,\nu)-E[f(X,\alpha,\nu)]\}/\sqrt n$. For each $\delta>0$ let $U_{0,n}(\delta)$ be the set

where $\varrho_n\ge 1$ and $\alpha_{0,n}\in\mathbf A_{0,n}$ (to be specified). Thus, $U_{0,n}(\delta)$ may be loosely viewed as a (partially punctured) neighborhood of $(\alpha_{0,n},\nu_0)$. We now impose:

Assumption (ref)(i) implies that $\Delta\Gamma_n$ and hence $\mathbf A_{0,n}$ are bounded by possibly diverging sequences. This is consistent with common choices of sieves (including neural networks). Assumption (ref)(ii) formalizes $\alpha_{0,n}$ as an approximation of $\alpha_0$, which may be taken as the maximizer of $\mathfrak C_0$ over $\mathbf A_{0,n}$ in certain settings; see Lemma (ref). Specific approximation errors are known for a variety of sieve spaces; see, e.g., Chen2007Handbook for classical sieves and Yarotsky2017Error, Yang2025ReLU, NaglerLanger2026Optimal, and references therein for deep neural networks.

Assumption (ref)(iii) regulates the penalization function $J_n$ and the parameter $\lambda_n$. Assumption (ref)(iv) is known as the modulus of continuity condition VaartWellner1996Book, which, intuitively speaking, restricts the “local size” of the sieve $\mathbf A_{0,n}$ at $\alpha_{0,n}$. In turn, the modulus $\omega_n$ determines the order $\delta_n$ via Assumption (ref)(v). We note that $\delta_n$ may be viewed as the critical radius, a concept that plays important roles in modern machine learning---see Remark (ref). Finally, Assumption (ref) introduces the sample and a consistent estimator $\hat\nu_n$ of $\nu_0$.

Theorem (ref) implies that one may employ a sieve $\mathbf A_{0,n}$ with a diverging bound while still having $\|\hat\alpha_n-\alpha_0\|_{\mathbf A}= O_p(\delta_n)$, as long as $\hat\nu_n$ converges to $\nu_0$ fast enough so that $\|\hat\nu_n-\nu_0\|_{\mathbf B}=O_p(\delta_n/\varrho_n)$. The sieve space $\mathbf A_{0,n}$ determines the convergence rate through its critical radius as well as its approximation capability. The convergence rate of $\hat\alpha_n$ is also impacted by the tuning parameter $\lambda_n$. Developing a data driven choice of $\lambda_n$ is an important issue that is beyond the scope of this paper.

Next, we turn to the case when $\Pi_0$ is unknown. Let $\mathbf D_n$ be a set of linear maps from $\mathbf H_0$ to $\mathbf A$ and define $\|\Pi\|_{op,n}\equiv \sup_{\eta\in\Delta\Gamma_n\colon\|\eta\|_{\mathbf H}\le 1}\|\Pi(\eta)\|_{\mathbf A}$ for any linear $\Pi\colon \mathbf H_0\to\mathbf A$. Accordingly, we redefine for each $\delta>0$ the set $U_{0,n}(\delta)$ as:

where $\eta_{0,n}\in\Delta\Gamma_n$ is such that $\alpha_{0,n}=\Pi_0(\eta_{0,n})$. We now introduce our final assumption in this section as follows.

Assumption (ref)(i) is a mild simplifying condition so that $\|\cdot\|_{op,n}$ may be defined over the unit ball in $\Delta\Gamma_n$. Assumption (ref)(ii) introduces a $\varrho_n$-consistent estimator $\hat\Pi_n$ of $\Pi_0$. Assumption (ref)(iii) is an analog of Assumption (ref)(iv) that accommodates unknown $\Pi_0$. Assumption (ref)(iv) is needed because our estimator $\hat\alpha_n$ below may lie outside the parameter space $\mathbf A_0$ of $\alpha_0$ (due to the estimator error in $\hat\Pi_n$). The condition $\hat\Pi_n(\Delta\Gamma_n)\subset\mathbf A_0$ with probability approaching one means that we may analyze $\hat\alpha_n$ as if it belongs to $\mathbf A_0$. Alternatively, the condition $\mathfrak C_0(\alpha)-\mathfrak C_0(\alpha_0) \lesssim\|\alpha-\alpha_0\|_{\mathbf A}^\tau$ whenever $d_{\mathbf A}(\alpha,\mathbf A_0)\le\epsilon$ regulates the curvature of $\mathfrak C_0$ over a small enlargement of $\mathbf A_0$ so that the estimator error of $\hat\Pi_n$ is under control. In Example (ref), this curvature condition automatically holds with $\tau=2$---see Appendix (ref) for details.

Given the estimator $\hat\Pi_n$, we estimate the sieve $\mathbf A_{0,n}$ by $\hat{\mathbf A}_n\equiv\{\hat\Pi_n(\eta)\colon \eta\in\Delta\Gamma_n\}$. In turn, we define the estimator $\hat\alpha_n\equiv \hat\Pi_n(\hat\eta_n)$ where $\hat\eta_n\in\Delta\Gamma_n$ is such that

where $\bar J_n\colon\Delta\Gamma_n\to\mathbf R_+$ relates to the previous penalty $J_n$ via $J_n(\Pi(\eta))=\bar J_n(\eta)$ for any $\Pi\in\mathbf D_n$ and $\eta\in\Delta\Gamma_n$, i.e., $J_n(\Pi(\eta))$ depends on $\Pi(\eta)$ only through $\eta$ (as is common in machine learning). As previously, the linear term in $\mathfrak C_0$ may be estimated automatically in some settings without requiring an estimator of $\nu_0$---see Remark (ref).

The next theorem obtains the convergence rate of $\hat\alpha_n$ when $\Pi_0$ is unknown.

Theorem (ref) shows that the convergence rate of $\hat\alpha_n$ is in addition controlled by the convergence rate of $\hat\Pi_n$. Note that Theorem (ref) includes Theorem (ref) as a special case by setting $\mathbf D_n=\{\Pi_0\}$. In Example (ref), $\Pi_0$ is the conditional expectation operator and so $\hat\Pi_n$ may be constructed by nonparametric regression methods. While the convergence rates of various classical nonparametric regression estimators are well understood BlundellChenKristensen2007Engel,DarollesFanFlorensRenault2011NPIV,Horowitz2011NPIV, those of modern machine learners are less established. Recent advances on convergence rates of nonparametric regression estimators obtained by deep learning include Schmidt2020ReLU, FarrellLiangMisra2021DNN, and KohlerLanger2021RateDNN. Finally, in certain settings such as NPIV models, the $\hat\nu_n$ term does not appear in the order because $\Upsilon_0$ does not involve $\nu_0$ and $\Psi_0$ may be estimated automatically as explained in Remark (ref).

Asymptotic Normality

While the asymptotic normality of the debiased GMM estimator $\hat\theta_n$ is a relatively standard result, we present it for completeness. To this end, we impose assumptions similar to but more primitive than those in ChernozhukovEscancianoIchimuraNewey2022LocalRobust as we separate analytic conditions from probabilistic ones. Following the literature ChernozhukovChetverikovDemirerDufloHansenNeweyRobins2018Double, we combine the orthogonal moment in (ref) with cross-fitting---the latter in particular helps further alleviate overfitting biases inherited from first step learners.

Formally, let $\bigcup_{k=1}^K I_k$ be a $K$-fold partition of $\{1,\ldots,n\}$ for some $K>1$ and denote by $n_k$ the sample size of the $k$-th fold. Suppose that $\hat\gamma_{-k}$ and $\hat\alpha_{-k}$ are estimators of $\gamma_0$ and $\alpha_0$ respectively that are based on observations {\it not} from the $k$-th fold. Further assume that $\tilde\theta_{-k}$ is a preliminary estimator of $\theta_0$ based on observations not in the $k$-th fold. For notational simplicity, define $\psi(X,\theta,\gamma,\alpha) \equiv g(X,\theta,\gamma) +\phi(X,\theta,\gamma,\alpha)$, $\psi_0(X)\equiv \psi(X,\theta_0,\gamma_0,\alpha_0)$, and, for a sample $\{X_i\}_{i=1}^n$,

Since each moment function is associated with a (potentially) different Riesz representer, $\alpha_0$ and $\hat\alpha_{-k}$ in this section are therefore understood as $d_g\times 1$ vectors of functions. For a $d_g\times d_g$ weighting matrix $\hat\Omega_n$, the debiased GMM estimator is then

Having introduced the notation, we impose the assumptions in this section.

Assumption (ref) places standard restrictions on the GMM model. Assumption (ref)(iv) in particular is a sufficient condition to deal with the triangular array nature of each fold $I_k$. Assumption (ref) imposes conditions on the augmented moment function. Assumption (ref)(i) is a minimal requirement. Assumption (ref)(ii) entails the orthogonality of $E[\psi(X,\theta_0,\cdot,\alpha_0)]$ at $\gamma_0$, i.e., $\nabla_\gamma E[\psi(X,\theta_0,\gamma_0,\alpha_0)][\gamma-\gamma_0]=0$ for all $\gamma\in V_0$. Given $E[\psi_0(X)]=0$ and the orthogonality, Assumption (ref)(ii) holds if $E[\psi(X,\theta_0,\cdot,\alpha_0)]$ is twice continuously Fr\'{e}chet differentiable around $\gamma_0$. Assumption (ref)(iii) is a global orthogonality property of the first step influence function ChernozhukovEscancianoIchimuraNewey2022LocalRobust. The continuity in Assumptions (ref)(iv)(v)(vi) follows if $g$ and $\phi$ are mean-square continuous at the truth. Assumption (ref)(vii) is mild since $E[\Delta_0(X,\theta,\gamma,\alpha)]$ may be viewed as a discrete version of the second derivative of $E[\phi(X,\cdot,\cdot,\cdot)]$ at $(\theta_0,\gamma_0,\alpha_0)$, with the discrepancy of order $o(\|\theta-\theta_0\|^2+\|\gamma-\gamma_0\|_{\mathbf H}^2+\|\alpha-\alpha_0\|_{\mathbf A}^2)$.

Assumption (ref) formalizes restrictions on the sample and relevant estimators. Assumption (ref)(i) is standard, while Assumption (ref)(ii) is a high level condition imposed for the sake of transparency of the theory---see Lemma (ref) for lower level conditions. Assumption (ref)(iii) demands the usual faster-than-$n^{-1/4}$ rate of convergence on the first step estimators $\hat\gamma_{-k}$ and $\hat\alpha_{-k}$ as well as the preliminary estimator $\tilde\theta_{-k}$. Note that when the augmented moment function is doubly robust (so that $E[\psi(X,\theta_0,\gamma,\alpha_0)]=E[\psi(X,\theta_0,\gamma_0,\alpha)]=0$) and $\|E[\Delta_0(X,\theta,\gamma,\alpha)]\|=O(\|\gamma-\gamma_0\|_{\mathbf H}\|\alpha-\alpha_0\|_{\mathbf A})$, it suffices to have $\|\hat\gamma_{-k}-\gamma_0\|_{\mathbf H} \|\hat\alpha_{-k}-\alpha_0\|_{\mathbf A}=o_p(n^{-1/2})$. Assumption (ref)(iv) is also minimal by requiring the sample size in each fold to be large. Assumption (ref)(v) is the common consistency requirement on the weighting matrix.

Assumption (ref) imposes additional conditions for variance estimation, which requires sample analogs of $\Sigma_0\equiv \mathrm{Var}[\psi_0(X)]$ and $G_0\equiv E[\nabla_\theta g(X,\theta_0,\gamma_0)]$ defined as:

for $\hat\psi_{-k}(X_i)\equiv \psi(X_i,\tilde\theta_{-k},\hat\gamma_{-k},\hat\alpha_{-k})$. In particular, Assumptions (ref)(i)(ii) and (iii) ensure the consistency of $\hat\Sigma_n$ and $\hat G_n$ respectively.

Proposition (ref) formally establishes the asymptotic normality of the debiased GMM estimator and provides a consistent estimator for its asymptotic variance.

Implementation and Simulations

In this section, we outline some implementation details and conduct simulation studies. In particular, we discuss how monotonicity and convexity may be enforced in neural networks, while noting that much work remains to develop architectures for other important constraints in economics, e.g., the Slutsky restriction, supermodularity, quasi-concavity, and joint constraints. We pay attention to finite sample gains from imposing shape constraints on first steps in both exogenous and endogenous settings.

Implementation

Algorithm (ref) is a general recipe for computing the debiased-GMM estimator $\hat\theta_n$ with optimal weighting and its variance estimator $\hat V_{\theta,n}$. It starts with splitting the sample into $K$-folds. As in ChernozhukovChetverikovDemirerDufloHansenNeweyRobins2018Double,ChernozhukovEscancianoIchimuraNewey2022LocalRobust,ChernozhukovNeweySingh2022Automatic, the general theory dictates that $K$ be small so that the sample size in each fold remains relatively large. In practice, one may choose $K=5$ or $K=10$ as recommended by ChernozhukovEscancianoIchimuraNewey2022LocalRobust,ChernozhukovNeweySingh2022Automatic. Given the split sample, we then construct preliminary estimators for the parameters $\gamma_0$, $\alpha_0$ and possibly also for $\theta_0$ which may enter both the weighting matrix $\Omega_0=\Sigma_0^{-1}$ and the first step influence function $\phi$. In just-identified models (i.e., $d_\theta=d_g$), there is no need to obtain the preliminary estimator $\tilde\theta_{-k}$ if $\phi$ does not involve $\theta_0$.

A key step in estimating $\gamma_0$ and $\alpha_0$ is incorporating the relevant shape constraint. A generic strategy is to post-process an unconstrained estimator BonakdarpourChatterjeeBarberLafferty2018Reshape,ChenChernozhukovFernandezKostyshakLuo2021Shape. While this approach works well in low dimensional settings, it becomes impractical as dimensionality increases. A strand of the neural network literature instead develops architectures that enforce shape constraints by construction, thereby yielding global guarantees and improved computational tractability. Since these architectures do not appear well known in economics, below we review those employed in our numerical exercises, assuming that the reader has a basic familiarity with neural networks---see, e.g., GoodfellowBengioCourville2016DL for an introduction.

The literature primarily focuses on monotonicity and convexity. Figure (ref)(a) presents the monotonic neural network of SartorSinigagliaSusto2025MNN in which (i) the weights are constrained to be nonnegative; (ii) the activation functions sequentially alternate between $\sigma\colon\mathbf R\to\mathbf R$ and its point reflection $\sigma^\dag(\cdot)\equiv-\sigma(-\cdot)$; (iii) $\sigma$ is weakly increasing and admits a finite limit at either $\infty$ or $-\infty$ (or both); and (iv) there are at least three hidden layers. Crucially, compared to prior monotonic networks, this architecture permits the use of unbounded activation functions (e.g., ReLU) while also retaining the universal approximation property. In many applications, however, one is interested in partially monotonic functions. Figure (ref) illustrates a two-branch neural network that is monotonic in a subvector $(z_1^{\mathrm{c}},\ldots,z_{p_{\mathrm{c}}}^{\mathrm{c}})$ of the input $z=(z_1,\ldots,z_{d_z})$.

Convexity requires a more careful treatment. Although a number of partially convex architectures have been developed, including AmosXuKolter2017ICNN, KimKim2022Convex, and SchallerBemporadBoyd2025Convex, they do not appear to perform well for inference, as documented in ChenChoiFang2026SCDL. We therefore adopt the two-branch convex architecture proposed by ChenChoiFang2026SCDL, where the convex branch is a pointwise max-affine subnet with multiple outputs, and the head is similar to the convex branch except that the inputs from the convex branch are constrained to be nonnegative so that the overall architecture is partially convex in $(z_1^{\mathrm{c}},\ldots,z_{p_{\mathrm{c}}}^{\mathrm{c}})$---see Figures (ref)(b) and (ref).

\tikzset{layer sep/.store in=\LayerSep,layer sep=0.70, neuron/.style={circle,draw,minimum size=5.9mm,inner sep=0pt}, networkblock/.style={draw,minimum width=3.2cm,minimum height=1.1cm, align=center,font=,inner sep=3pt}, networkheadblock/.style={draw,minimum width=4.1cm,minimum height=3.15cm, align=center,font=,inner sep=3pt}, nodeunit/.style={minimum size=5.9mm,inner sep=0pt, text height=1.25ex,text depth=1.25ex,font=}, labin/.style={font=,inner sep=0pt}, paneltitle/.style={font=,inner sep=0pt}, conn/.style={-{Stealth[length=1.3mm,width=0.9mm]},ultra thin}, network dots/.pic={ \foreach \dy in {-0.08,0,0.08}{\fill (0,\dy) circle[radius=0.45pt];} }} \makeatletter \pgfplotsset{ colormap/BuGn-9, auto title/.style={ title={(\@alph{\pgfplots@group@current@plot}) #1}, }, NetworkGroup/.style={ group style={horizontal sep=0pt}, width=0.46\textwidth, height=0.40\textwidth, xmin=-2.0,xmax=2.0, ymin=-1.45,ymax=2.85, hide axis, clip=false, scale only axis, title style={paneltitle,at={(axis description cs:0.5,-0.11)},anchor=north}, }, } \makeatother

Simulation Studies

In this section, we conduct Monte Carlo simulations to illustrate our estimation procedure and demonstrate the potential gains of imposing shape constraints. Below, we let $\Phi$ be the standard normal cdf and denote by $\mathrm{Uni}[0,1]$, $\mathrm{Ber}(\varrho)$, and $N(0,1)$ the uniform distribution on $[0,1]$, the Bernoulli distribution with success probability $\varrho$, and the standard normal distribution respectively.

Average Treatment Effects

Let $Z\in\mathbf R^{d_z}$ be a vector of i.i.d.\ entries drawn from $\mathrm{Uni}[0,1]$, $D\sim \mathrm{Ber}(\varrho_0(Z))$ conditional on $Z$, $U\sim N(0,1)$ drawn independently of $D$ and $Z$, and

Let $Z_{\mathrm{c}}\in\mathbf R^{d_{\mathrm{c}}}$ and $Z_{\mathrm{f}}\in\mathbf R^{d_{\mathrm{f}}}$ denote the constrained and free inputs respectively, and, for $Z=(Z_{\mathrm{c}}^\intercal,Z_{\mathrm{f}}^\intercal)^\intercal$, set $\varrho_0(Z)=1/[1+\exp\{-Z^\intercal\mathsf a_\varrho\}]$ and

where $q(Z_{\mathrm{f}})$ collects all second-degree polynomial terms of $Z_{\mathrm{f}}$ (squares and pairwise interactions), and each of the vectors $\mathsf a_\varrho$, $\mathsf a_{\mu,\mathrm{c}}$, $\mathsf b_{\mu,\mathrm{c}}$, $\mathsf a_{\mu,\mathrm{f}}$, $\mathsf b_{\mu,\mathrm{f}}$, $\mathsf a_{\tau,\mathrm{c}}$, $\mathsf b_{\tau,\mathrm{c}}$, and $\mathsf a_{\tau,\mathrm{f}}$ consists of i.i.d.\ entries drawn from uniform distributions whose supports are given in Table (ref). By design, both $\mu_0(Z)$ and $\tau_0(Z)$ are monotone and concave in $Z_{\mathrm{c}}$. Consequently, $(D,Z)\mapsto\gamma_0(D,Z)=\mu_0(Z)+D\tau_0(Z)$ satisfies the same shape restrictions. In the simulations, we impose monotonicity and concavity separately in the first step estimation to investigate their respective impacts on the performance of the resulting debiased estimators. To impose concavity, we apply the convex architecture to estimate $-\gamma_0$ and then flip the sign of the resulting estimator.

We also consider deviations from the above baseline design. First, we vary the signal-to-noise ratio by drawing $U \sim N(0,\sigma_U^2)$ with $\sigma_U^2 \in \{0.5,1,5\}$. The case $\sigma_U^2=1$ serves as the baseline, while $\sigma_U^2=0.5$ and $\sigma_U^2=5$ correspond to higher and lower signal-to-noise ratios, respectively. Second, we vary the “strength” of monotonicity/concavity by considering three supports for the entries of $\mathsf a_{\mu,\mathrm{c}}$, $\mathsf b_{\mu,\mathrm{c}}$, $\mathsf a_{\tau,\mathrm{c}}$, and $\mathsf b_{\tau,\mathrm{c}}$: $[0.1,0.2]$ (weak), $[0.2,0.3]$ (baseline), and $[0.3,0.4]$ (strong). For brevity, we report only the results for the baseline specification with $\sigma_U^2=1$. The results under $\sigma_U^2\in\{0.5,5\}$ are qualitatively similar, suggesting that our findings are robust to alternative signal-to-noise ratios. These additional results are available upon request.

For each design, we draw $5000$ independent i.i.d.\ samples of size $1000$, each with $d_z=40$ covariates that comprise $d_{\mathrm{c}}=20$ constrained ones and $d_{\mathrm{f}}=20$ free ones. We then estimate $\gamma_0$ and $\alpha_0$ using partially monotonic/convex neural networks, varying the number of imposed constrained inputs over $p_{\mathrm{c}}=0,1,\ldots,20$. In the partially monotonic architecture, both the monotone and free branches consist of two hidden layers with width $32$ for each and $32$ outputs, while the head network has two hidden layers with the same width for each. In the partially convex architecture, the free branch consists of three hidden layers with width $128$ for each and $8$ outputs, the convex branch consists of $32$ affine nodes with $4$ outputs, while the head has $128$ affine nodes.

Training of the networks is carried out using the AdamW optimizer LoshchilovHutter2018Decoupled, a state-of-the-art variant of stochastic gradient descent. We further incorporate dropout, early stopping, and $\ell^2$-penalization for the sake of regularization. For each sample, we construct $K=5$ folds such that treated and untreated individuals are approximately evenly distributed across the folds. Given existing theoretical and simulation studies (see, e.g., ChernozhukovNeweyQuintasSyrgkanis2024RieszReg and references therein), we benchmark our results against DML via Riesz regression but without shape constraints, which corresponds to $p_{\mathrm{c}}=0$ in the plots below.

\pgfplotstableread{ num_mon bias std mean se coverage length 0 0.009039872 0.072466711 1.381145572 0.071396392 0.952 0.279873857 1 0.000266501 0.070234476 1.372372201 0.070270935 0.947 0.275462063 5 0.002036933 0.069527558 1.374142633 0.069300152 0.9462 0.271656596 9 0.022116964 0.069382357 1.394222664 0.069177566 0.9374 0.27117606 13 0.012232518 0.06857148 1.384338218 0.068721717 0.9472 0.269389129 17 0.016108279 0.068317802 1.388213979 0.068128666 0.9414 0.267064371 20 0.011119099 0.067298628 1.383224799 0.067299259 0.9474 0.263813094 }\MonLM

\pgfplotstableread{ num_mon bias std mean se coverage length 0 0.007204591 0.073929027 2.013739087 0.073076758 0.9482 0.286460891 1 -0.004257891 0.072093505 2.002276606 0.071664658 0.9474 0.280925458 5 -0.002056687 0.070354959 2.004477809 0.069958324 0.9458 0.27423663 9 0.01161494 0.069722074 2.018149436 0.069426529 0.9428 0.272151993 13 0.006463038 0.069340462 2.012997534 0.068891871 0.9464 0.270056133 17 0.010298526 0.068708109 2.016833023 0.068281159 0.9438 0.267662145 20 0.00921005 0.067774701 2.015744546 0.067630856 0.9472 0.265112955 }\MonMM

\pgfplotstableread{ num_mon bias std mean se coverage length 0 0.004767242 0.07633053 2.607255209 0.075294589 0.9508 0.29515479 1 -0.00881513 0.072886081 2.593672837 0.073390909 0.9476 0.287692365 5 -0.006978 0.071773499 2.595509967 0.071036219 0.9408 0.27846198 9 0.003671135 0.07078449 2.606159101 0.070159375 0.9474 0.27502475 13 0.000102286 0.069905226 2.602590253 0.069331072 0.9492 0.2717778 17 0.002865768 0.069523704 2.605353734 0.0687721 0.9446 0.269586633 20 0.006670849 0.068737326 2.609158815 0.068239404 0.9466 0.267498463 }\MonHM

\pgfplotstableread{ num_mon bias std mean se coverage length 0 0.009039872 0.072466711 1.381145572 0.071396392 0.952 0.279873857 1 0.027490055 0.070429752 1.399595755 0.070936667 0.9348 0.278071734 5 0.013296048 0.069357202 1.385401748 0.069994264 0.9502 0.274377515 9 0.006107393 0.069132457 1.378213093 0.069608902 0.9466 0.272866895 13 0.003929105 0.068376285 1.376034806 0.069223481 0.9496 0.271356047 17 0.003258974 0.067923819 1.375364674 0.068949095 0.9526 0.270280452 20 0.005902683 0.067105004 1.378008383 0.068780165 0.951 0.269618247 }\ConLM

\pgfplotstableread{ num_mon bias std mean se coverage length 0 0.007204591 0.073929027 2.013739087 0.073076758 0.9482 0.286460891 1 0.026496917 0.071428294 2.033031414 0.07158757 0.9344 0.280623275 5 0.010887851 0.070784223 2.017422348 0.070813949 0.9508 0.277590678 9 0.003194483 0.070317319 2.009728979 0.070465918 0.9464 0.276226399 13 0.004310014 0.069594818 2.01084451 0.070066682 0.9494 0.274661395 17 0.005456295 0.068750468 2.011990792 0.0697093 0.9518 0.273260455 20 0.007150389 0.068319208 2.013684886 0.069475436 0.952 0.272343709 }\ConMM

\pgfplotstableread{ num_mon bias std mean se coverage length 0 0.004767242 0.07633053 2.607255209 0.075294589 0.9508 0.29515479 1 0.02398886 0.072251341 2.626476827 0.072069217 0.9382 0.282511329 5 0.005773408 0.072026708 2.608261374 0.071776503 0.948 0.281363892 9 -0.002663477 0.071276426 2.599824489 0.071464746 0.9458 0.280141804 13 0.003718299 0.070855048 2.606206266 0.071112049 0.9472 0.278759232 17 0.008086438 0.069885816 2.610574405 0.071041648 0.9532 0.278483261 20 0.008900983 0.069606871 2.611388949 0.070940508 0.9522 0.27808679 }\ConHM

\pgfplotsset{ GroupCommon/.style={ group style={horizontal sep=0.8cm,vertical sep=1cm}, scaled y ticks=false, xmin=0, xmax=20, xtick={0,1,5,9,13,17,20}, xlabel={$p_{\mathrm{c}}$}, x label style={at={(axis description cs:0.95,0)},anchor=south,font=}, tick label style={/pgf/number format/fixed,font=}, yticklabel style={rotate=90,/pgf/number format/fixed,/pgf/number format/precision=3,/pgf/number format/fixed zerofill}, enlarge x limits={abs=0.5}, enlarge y limits, legend style={cells={align=center},column sep=3pt,row sep=3pt,draw=none}, cycle list={ {smooth,tension=0.5,color=Paired-B, mark=10-pointed star,mark size=1.75pt,line width=0.5pt}, {smooth,tension=0.5,color=Dark2-B, mark=halfsquare*,mark size=1.75pt,line width=0.5pt}, {smooth,tension=0.5,color=Paired-D, mark=halfcircle*,mark size=1.75pt,line width=0.5pt}, {smooth,tension=0.5,color=Paired-J, mark=triangle*,mark size=1.75pt,line width=0.5pt}, } }, }

Figures (ref) and (ref) summarize the ATE simulations for the baseline signal-to-noise ratio based on different supports for the coefficients $\mathsf a_{\mu,\mathrm{c}}$, $\mathsf b_{\mu,\mathrm{c}}$, $\mathsf a_{\tau,\mathrm{c}}$, and $\mathsf b_{\tau,\mathrm{c}}$. In both monotonicity and concavity designs, the bias remains small relative to sampling variability, coverage is close to the nominal level, and the standard deviation and confidence-interval length generally decrease as more shape restrictions are imposed. These patterns suggest that shape restrictions can serve as substantive regularization, improving estimation precision and inference accuracy while maintaining small bias.

We emphasize that a uniform tuning scheme is used across the number $p_{\mathrm{c}}$ of constrained inputs. In particular, the choices of network architectures, dropout rates, learning rates, early stopping criteria, and penalty levels are neither data-driven nor necessarily optimal. The resulting performance is therefore particularly encouraging. It would be of interest to develop principled guidance for selecting these tuning parameters that both admits theoretical guarantees and performs well empirically.

Functionals of NPIV Regression

Let $(D,Z_{\mathrm{c}}^{\intercal},Z_{\mathrm{f}}^{\intercal})^{\intercal}\in\mathbf{R}^{1+d_{\mathrm{c}}+d_{\mathrm{f}}}$ be the vector of endogenous inputs and $(S,W_{\mathrm{c}}^{\intercal},W_{\mathrm{f}}^{\intercal})^{\intercal}\in\mathbf{R}^{1+d_{\mathrm{c}}+d_{\mathrm{f}}}$ the vector of instruments such that

To build the data generating process, we let $U_0,V_0,U_{\mathrm{c}},U_{\mathrm{f}}$, and $\{V_{\mathrm{c},j},V_{\mathrm{f},j}\colon j\ge 1\}$ be independent standard normal random variables and set $U=(U_0+U_{\mathrm{c}}+U_{\mathrm{f}})/\sqrt{3}$,

$S=\Phi(V_0)$, $W_{\mathrm{c},1}=\Phi(V_{\mathrm{c},1})$, $W_{\mathrm{f},1}=\Phi(V_{\mathrm{f},1})$, and, for $j\ge2$,

Intuitively, $\eta$ controls the strength of endogeneity while $\rho$ controls the correlation between $\{D,Z_{\mathrm{c},1},Z_{\mathrm{f},1}\}$ and $\{Z_{\mathrm{c},j},Z_{\mathrm{f},j}\colon j\geq 2\}$. In turn, we set

where $\mathsf a_{\mathrm{c}}$ and $\mathsf b_{\mathrm{c}}$ have i.i.d.\ entries drawn from $\mathrm{Uni}[0.2,0.3]$, $\mathsf a_{\mathrm{f}}$ has i.i.d.\ entries drawn from $\mathrm{Uni}[-0.1,0.1]$, and $\mathsf b_{\mathrm{f}}$ has i.i.d.\ entries drawn from $\mathrm{Uni}[-0.05,0.05]$. By construction, $\gamma_0$ is monotonic in $Z_{\mathrm{c}}$, and in this NPIV design we focus on monotonicity.

The parameter of interest is the weighted average derivative: for $Z\equiv(Z_{\mathrm{c}}^{\intercal},Z_{\mathrm{f}}^{\intercal})^\intercal$,

where $w$ is a weight function and $\nabla_d$ indicates the partial derivative with respect to $D$. Intuitively, $\theta_0$ measures the average marginal/treatment effect of a small change in $D$ on the outcome $Y$ ImbensNewey2009Triangular. Suppose that $D$ is supported on $[a,b]$ with $-\infty<a<b<\infty$ and that $d\mapsto w(d)f_0(d,z)$ vanishes on $a$ and $b$ for (almost) every $z$ with $f_0$ the density of $(D,Z)$. Then we obtain under regularity conditions that

which is continuous and linear in $\gamma_0$ if $w(D)\nabla_d\ln(w(D)f_0(D,Z))\in L^2(D,Z)$.

The main comparative statics vary the endogeneity level $\eta$ from $0.1$ to $1.0$ based on $\rho=0.1$. We set $w$ to be the normalized indicator of the $[0.1,0.9]$ quantile region of $D$ so that $\theta_0=1$. Relative to the ATE designs, the additional unknown object is the projection operator $\Pi_0$ given by $\Pi_0(\delta)=E[\delta(D,Z_{\mathrm{c}},Z_{\mathrm{f}})|S,W_{\mathrm{c}},W_{\mathrm{f}}]$ for any $\delta\in L^2(D,Z_{\mathrm{c}},Z_{\mathrm{f}})$. We obtain $\hat\Pi_n(\delta)$ by a ridge regression of $\delta$ on a vector of basis functions of $(S,W_{\mathrm{c}},W_{\mathrm{f}})$. We consider two constructions of these basis functions. The first one is an adaptively learned basis inspired by Bruns2025TSML. Specifically, we fit a deep neural network to predict $Y$ from $(S,W_{\mathrm{c}},W_{\mathrm{f}})$ and employ the activations in its final hidden layer as basis functions. The second one is a specific choice of spline basis functions used in ChenChenTamer2023NPIV.

To ease computation, we consider two specifications, one with $d_{\mathrm{c}}=d_{\mathrm{f}}=1$ and one with $d_{\mathrm{c}}=d_{\mathrm{f}}=5$, and then draw 500 independent i.i.d.\ samples of size 1000 for each design. We estimate the nuisances $\gamma_0$ and $\alpha_0$ using the same partially monotonic neural network architecture as in Section (ref), varying the number of imposed constrained inputs over $p_{\mathrm{c}}=0,1,\ldots,d_{\mathrm{c}}$. Training of the networks is carried out using the same AdamW optimizer, dropout, early stopping, and $\ell_2$-penalization as in the ATE simulations. For each sample, we construct $K=5$ folds for cross-fitting. As before, we benchmark our results against DML via Riesz regression but without shape constraints, which corresponds to $p_{\mathrm{c}}=0$ in the plots below.

\pgfplotstableread{ eta pm0 pm1 0.1 0.504000000 0.954000000 0.2 0.506000000 0.948000000 0.3 0.498000000 0.948000000 0.4 0.496000000 0.942000000 0.5 0.496000000 0.952000000 0.6 0.510000000 0.940000000 0.7 0.498000000 0.928000000 0.8 0.490000000 0.942000000 0.9 0.504000000 0.942000000 1.0 0.524000000 0.950000000 }\NPIVOneCoverage

\pgfplotstableread{ eta pm0 pm1 0.1 -0.062901740 -0.029790198 0.2 -0.067800841 -0.028454479 0.3 -0.072203688 -0.025911979 0.4 -0.077364539 -0.029073808 0.5 -0.082567237 -0.027291200 0.6 -0.092458112 -0.025538096 0.7 -0.098472011 -0.023261110 0.8 -0.108009344 -0.023605807 0.9 -0.114783656 -0.022755386 1.0 -0.122451760 -0.025255426 }\NPIVOneBias

\pgfplotstableread{ eta pm0 pm1 0.1 0.029913114 0.013890082 0.2 0.030633043 0.014574314 0.3 0.033374373 0.015511749 0.4 0.036130226 0.016587384 0.5 0.040420215 0.017691700 0.6 0.045231303 0.020028298 0.7 0.050188412 0.022012582 0.8 0.055315317 0.023618394 0.9 0.059976672 0.026304121 1.0 0.064264659 0.029091482 }\NPIVOneMSE

\pgfplotstableread{ eta pm0 pm1 pm3 pm5 0.1 0.092000000 0.698000000 0.748000000 0.794000000 0.2 0.074000000 0.688000000 0.742000000 0.810000000 0.3 0.070000000 0.666000000 0.726000000 0.774000000 0.4 0.052000000 0.628000000 0.692000000 0.752000000 0.5 0.042000000 0.622000000 0.694000000 0.708000000 0.6 0.032000000 0.594000000 0.658000000 0.704000000 0.7 0.034000000 0.568000000 0.622000000 0.678000000 0.8 0.032000000 0.542000000 0.588000000 0.654000000 0.9 0.024000000 0.526000000 0.570000000 0.624000000 1.0 0.016000000 0.488000000 0.554000000 0.572000000 }\NPIVFiveCoverage

\pgfplotstableread{ eta pm0 pm1 pm3 pm5 0.1 -0.261582946 -0.057512979 -0.046756125 -0.040874387 0.2 -0.277904062 -0.064951272 -0.049404491 -0.041663732 0.3 -0.302376437 -0.070214754 -0.056650298 -0.049620320 0.4 -0.327251424 -0.075996287 -0.064755355 -0.054990715 0.5 -0.361374658 -0.085708130 -0.075923315 -0.066349879 0.6 -0.399536190 -0.099981222 -0.086775860 -0.078157370 0.7 -0.435843352 -0.115162929 -0.099078348 -0.091002018 0.8 -0.470035469 -0.132825014 -0.115458741 -0.103403743 0.9 -0.501143920 -0.154898285 -0.131658743 -0.117845826 1.0 -0.533352633 -0.178865136 -0.151951597 -0.136498229 }\NPIVFiveBias

\pgfplotstableread{ eta pm0 pm1 pm3 pm5 0.1 0.109574737 0.021298294 0.019640862 0.018855358 0.2 0.119968799 0.023154430 0.020243381 0.018749997 0.3 0.137380778 0.024828199 0.022833140 0.021257430 0.4 0.155227120 0.027946675 0.026323367 0.023673684 0.5 0.179837904 0.031125245 0.029017552 0.028058314 0.6 0.210131340 0.036028697 0.033435727 0.032539123 0.7 0.241718314 0.042074721 0.040455180 0.038274088 0.8 0.274076424 0.049555720 0.047203496 0.044008821 0.9 0.302487071 0.059079170 0.055331926 0.051948342 1.0 0.334651062 0.070926121 0.064312499 0.061423557 }\NPIVFiveMSE

\pgfplotsset{ NPIVGroupCommon/.style={ group style={horizontal sep=0.8cm,vertical sep=1.25cm}, scaled y ticks=false, xmin=0.1, xmax=1.0, xtick={0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0}, xlabel={$\eta$}, x label style={at={(axis description cs:0.95,0)},anchor=south,font=}, tick label style={/pgf/number format/fixed,font=}, yticklabel style={rotate=90,/pgf/number format/fixed,/pgf/number format/precision=3,/pgf/number format/fixed zerofill}, enlarge x limits={abs=0.02}, enlarge y limits, legend columns=-1, legend style={cells={align=center},column sep=6pt,row sep=0pt,draw=none}, cycle list={ {smooth,tension=0.5,color=Paired-B, mark=10-pointed star,mark size=1.75pt,line width=0.5pt}, {smooth,tension=0.5,color=Dark2-B, mark=halfsquare*,mark size=1.75pt,line width=0.5pt}, {smooth,tension=0.5,color=Paired-D, mark=halfcircle*,mark size=1.75pt,line width=0.5pt}, {smooth,tension=0.5,color=Paired-J, mark=triangle*,mark size=1.75pt,line width=0.5pt}, } }, }

Figure (ref) reports the simulation results based on the adaptive basis. Overall, imposing monotonicity substantially improves coverage probability while reducing bias and mean squared error across endogeneity levels relative to the unconstrained benchmark. When $d_{\mathrm{c}}=d_{\mathrm{f}}=1$, adding the single monotonicity restriction raises coverage from about $50\%$ to close to the nominal $95\%$ level. When $d_{\mathrm{c}}=d_{\mathrm{f}}=5$, coverage probabilities increase substantially as more monotonicity restrictions are imposed, rising from below $10\%$ for the unconstrained benchmark to roughly 50--80% when all five monotonicity restrictions are imposed. The bias and mean squared error also decline with higher values of $p_{\mathrm{c}}$, indicating clear gains in estimation accuracy in the high dimensional case. Figure (ref) shows that the spline basis implementation delivers similar findings. We note that achieving accurate coverage probabilities in high dimensional NPIV models is notoriously difficult---see, e.g., BennettKallusMaoNeweySyrgkanisUehara2025StrongID and Bruns2025TSML for related simulation evidence. Our results suggest that imposing shape constraints can help mitigate this challenge in these high dimensional inverse problems.

\pgfplotstableread{ eta pm0 pm1 0.1 0.392000000 0.904000000 0.2 0.392000000 0.910000000 0.3 0.376000000 0.906000000 0.4 0.366000000 0.900000000 0.5 0.372000000 0.908000000 0.6 0.348000000 0.890000000 0.7 0.360000000 0.896000000 0.8 0.362000000 0.880000000 0.9 0.362000000 0.894000000 1.0 0.356000000 0.876000000 }\NPIVSplineOneCoverage

\pgfplotstableread{ eta pm0 pm1 0.1 -0.105117650 -0.004186915 0.2 -0.110102383 -0.004988590 0.3 -0.119382818 -0.008409777 0.4 -0.127541303 -0.007055514 0.5 -0.135488984 -0.010703606 0.6 -0.146909398 -0.016768821 0.7 -0.156998746 -0.020608792 0.8 -0.169641960 -0.024061935 0.9 -0.178635922 -0.028447377 1.0 -0.188881501 -0.037457928 }\NPIVSplineOneBias

\pgfplotstableread{ eta pm0 pm1 0.1 0.039067700 0.015671467 0.2 0.041118862 0.014893703 0.3 0.046010553 0.015872737 0.4 0.051120141 0.016665281 0.5 0.054559042 0.017371498 0.6 0.061881182 0.020172121 0.7 0.067125027 0.022230449 0.8 0.073843113 0.025237452 0.9 0.079824932 0.027235375 1.0 0.085710711 0.031353640 }\NPIVSplineOneMSE

\pgfplotstableread{ eta pm0 pm1 pm3 pm5 0.1 0.092000000 0.612000000 0.642000000 0.700000000 0.2 0.080000000 0.606000000 0.650000000 0.668000000 0.3 0.062000000 0.558000000 0.638000000 0.664000000 0.4 0.044000000 0.520000000 0.586000000 0.646000000 0.5 0.042000000 0.516000000 0.558000000 0.574000000 0.6 0.028000000 0.472000000 0.516000000 0.570000000 0.7 0.020000000 0.396000000 0.480000000 0.492000000 0.8 0.016000000 0.378000000 0.446000000 0.448000000 0.9 0.016000000 0.324000000 0.398000000 0.428000000 1.0 0.010000000 0.284000000 0.350000000 0.372000000 }\NPIVSplineFiveCoverage

\pgfplotstableread{ eta pm0 pm1 pm3 pm5 0.1 -0.272919339 -0.055662661 -0.050733471 -0.048711623 0.2 -0.289564225 -0.062605444 -0.053068998 -0.058637260 0.3 -0.318010973 -0.076724572 -0.062639810 -0.062741916 0.4 -0.352754749 -0.089463980 -0.077486646 -0.072945493 0.5 -0.383779067 -0.101679530 -0.090679725 -0.088697680 0.6 -0.427628349 -0.126041521 -0.111128489 -0.108037349 0.7 -0.472427572 -0.150115667 -0.130843867 -0.130986624 0.8 -0.511157828 -0.174339298 -0.156925132 -0.156850090 0.9 -0.548284021 -0.207918637 -0.181296992 -0.180934853 1.0 -0.580650960 -0.239008996 -0.215794017 -0.217189526 }\NPIVSplineFiveBias

\pgfplotstableread{ eta pm0 pm1 pm3 pm5 0.1 0.113631430 0.020202204 0.020279492 0.020462485 0.2 0.125668490 0.021904933 0.021499798 0.022526410 0.3 0.145563844 0.025318539 0.024270672 0.024128664 0.4 0.171730738 0.029544026 0.029059087 0.027937183 0.5 0.196141021 0.033825982 0.033054951 0.033350312 0.6 0.231426659 0.041128897 0.039113649 0.039851514 0.7 0.269379741 0.052006776 0.046892182 0.048758793 0.8 0.307598293 0.063650828 0.060019823 0.059927726 0.9 0.341110076 0.079129346 0.071728581 0.068506644 1.0 0.375602957 0.094874464 0.086709819 0.089780207 }\NPIVSplineFiveMSE

Empirical Applications

In this section, we present two empirical applications. The first is a DiD analysis of Medicaid expansions and mortality, while the second studies the intertemporal relationship between working hours and wage growth. These applications illustrate the role of shape constraints as a middle ground between flexible nonparametric estimation and restrictive parametric modeling. Fully nonparametric methods can accommodate rich heterogeneity but may lack statistical guarantees in high dimensional settings. Parametric specifications are more tractable but rely on functional form assumptions that may be too restrictive. Shape constraints impose substantive restrictions that can improve precision while retaining substantial flexibility. The two applications complement one another by illustrating how shape constraints can be motivated either by institutional knowledge or by economic theory.

DiD under Shape Constraints

We proceed by introducing the DiD setup and formalizing the estimation framework as a specialization of our general theory. We then bring the model to the data by revisiting the Medicaid and mortality example in BakerCallawayCunninghamGoodmanSantAnna2025Guide.

The DiD Setup

We consider a staggered treatment regime with a binary treatment, i.e., once a unit is treated, it remains treated in all subsequent periods. Thus, an individual's treatment path is fully characterized by the period in which treatment first occurs. Specifically, suppose that there are $T$ periods indexed by $t=1,\ldots,T$ and let $G$ be the time when an individual is first treated. Then individuals may be classified into different groups based on the support $\mathcal G$ of $G$. We assume that no individuals are treated at $t=1$ and set $G=\infty$ if an individual never receives the treatment. Without loss of generality, we also assume the existence of a never-treated group by, if necessary, discarding observations during the periods when the last cohort is treated; otherwise, treatment effects in those periods would not be identified due to the absence of comparison groups. Given the absorbing nature of the treatment, the path of potential outcomes for an individual from group $G=g$ may be denoted as $\{Y_t(g)\}_{t=1}^T$. The observed outcome $Y_t$ at time $t$ relates to the potential outcomes via $Y_t=\sum_{g\in\mathcal G} Y_t(g)1\{G=g\}$. In addition, we observe a vector $Z$ of time invariant (pre-treatment) covariates.

Following CallawaySantAnna2021MultipleDID, we define the average treatment effect for group $g$ at time $t$ (GT-ATT) by

These GT-ATTs are of interest in their own right while also serving as building blocks for inference on more aggregated causal effects. For example, letting $e\equiv t-g$ be the event time, we may consider the event study parameter:

Intuitively, $\theta_{\mathrm{es}}(e)$ measures the average treatment effect among groups that have been exposed to treatment for $e$ time periods.

In order to identify $\theta_0(g,t)$, we must select a collection of comparison groups which we denote by $\mathcal C_{gt}$. A common choice is $\mathcal C_{gt}=\{s\in\mathcal G\colon s>g\vee t\}$, i.e., groups not yet treated at time $t$, although it may be desirable to choose $\mathcal C_{gt}$ as a proper subset of $\{s\in\mathcal G\colon s>g\vee t\}$. Suppose that (i) $Y_t(g)=Y_t(\infty)$ for all $t\le g-1$ and $g\in\mathcal G$ (no anticipation), (ii) $P(G=g|Z)>0$ for all $g\in\mathcal G$ (overlap), and (iii) conditional parallel trends hold in the sense that, for any $g\in\mathcal G\backslash\{\infty\}$, any $t\ge g$, and any $c\in\mathcal C_{gt}\subset \{s\in\mathcal G\colon s>g\vee t\}$,

Then it follows that, for any $g\in\mathcal G\backslash\{\infty\}$ and $t\ge g$,

where $\mu_{gt,0}(C_{gt},Z)\equiv E[Y_t-Y_{g-1}|C_{gt},Z]$ for $C_{gt}\equiv 1\{G\in\mathcal C_{gt}\}$. Equation (ref) is essentially the identification of $\theta_0(g,t)$ in CallawaySantAnna2021MultipleDID based on outcome regression using not yet treated groups as controls.

Given the identification in (ref), we may embed the setup into our framework by identifying the data vector $X$ with $(Y_{g-1},Y_t, G, C_{gt},Z)$, the first step $\gamma_0$ with $(\eta_{g,0},\mu_{gt,0})$ for $\eta_{g,0}\equiv P(G=g)$, and the moment function $g$ with

Note that (ref) is already orthogonal to $\eta_{g,0}$. By ChernozhukovNeweySingh2022Automatic, we may then orthogonalize (ref) with respect to $\gamma_0$ by setting

for any $\gamma\equiv (\eta,\mu)$ with $\eta\in(0,1)$ and $\mu\colon \{0,1\}\times\mathcal Z\to\mathbf R$ and any $\alpha\colon\{0,1\}\times\mathcal Z\to\mathbf R$, where the truth of $\alpha$ (i.e., the Riesz representer) is

To incorporate shape constraints, we note that although, in principle, they can be imposed directly on $\mu_{gt,0}$ (the conditional expectation of the trend $Y_t-Y_{g-1}$), it may be easier to motivate them at the levels. For example, if $E[Y_t|C_{gt},Z]$ and $E[Y_{g-1}|C_{gt},Z]$ are monotonic with respect to (a subset of) $Z$, then we may let the parameter space $\Gamma$ of $\mu_{gt,0}$ be $\Gamma=\Gamma^{\uparrow}-\Gamma^{\uparrow}$ where $\Gamma^{\uparrow}\subset L^2(\{0,1\}\times\mathcal Z)$ is a set of monotonic functions. In this case, our theory implies that

Alternatively, we note that we may replace $\mu(C_{gt},Z)$ with $\mu(1,Z)$ in (ref) so that we may simplify the estimation of the first step by training $\mu_{gt,0}(1,Z)$ using observations from the comparison groups. Accordingly, we set $\alpha_{gt,0}(C_{gt},Z)=C_{gt}\alpha^\flat_{gt,0}(Z)$ where $\alpha^\flat_{gt,0}(Z)\equiv -P(G=g|Z)/[P(G=g)P(C_{gt}=1|Z)]$ solves

where $\Gamma_\flat^\uparrow\subset L^2(\mathcal Z)$ is a set of monotonic functions $\alpha^\flat\colon\mathcal Z\to\mathbf R$. This is the specification we use to impose shape constraints in the empirical application below.

Given a sample $\{X_i\}_{i=1}^n$ of $X$, equation (ref) implies that we may estimate $\theta_0(g,t)$ based on observations at time $t$ consisting of units from group $g$ and $\mathcal C_{gt}$:

where $\hat\eta_{g,n}\equiv\sum_{i=1}^{n}1\{G_i=g\}/n$ is the sample analog of $\eta_{g,0}$, $\hat\mu_{gt,n}(1,\cdot)$ is an estimator of $\mu_{gt,0}(1,\cdot)$, and $\hat\alpha_{gt,n}$ is obtained by running the Riesz regression.

Medicaid and Mortality

Medicaid was enacted in 1965 as a public health insurance program providing coverage to low-income individuals. The Affordable Care Act (ACA), passed by Congress in 2010, expanded Medicaid eligibility to all adults with incomes up to $138\%$ of the federal poverty level. The ACA initially mandated expansion in all states, but a subsequent Supreme Court ruling rendered it optional, resulting in a staggered treatment design---see Table (ref) for the expansion timing. A central question is whether the Medicaid expansion causally reduces mortality. To address this, we build on the county level panel data constructed by BakerCallawayCunninghamGoodmanSantAnna2025Guide in which the outcome is the crude mortality rate (per 100,000) among adults aged 20--64 and the covariates consist of the percentages of the population that are respectively female, white, and Hispanic, as well as the unemployment rate, the poverty rate, and median income. The panel includes 2,604 counties and spans the period 2009--2019. The group index $G_i\in\{2014,2015,2016,2019,\infty\}$---recall that $G_i=\infty$ means that county $i$ did not expand the coverage by 2019. The pre-treatment periods allow us to conduct pre-trend analysis.

In exploring shape constraints in the present setting, we note that there is well documented empirical evidence supporting monotonic relationships between mortality and two key covariates: the poverty rate and income. For example, based on county level data for 1990, 2000, and 2010, CurrieSchwandt2016Mortality document that mortality rates increase with the poverty rate and decrease with median income. Moreover, both relationships are robust across gender and age groups. These findings align with GordonSommers2016Recessions, who report that higher poverty rates and lower median incomes are associated with higher mortality across multiple causes of death and adult subpopulations. Using administrative data at the individual level from 2001--2014, ChettyEtAl2016IncomeLife show that life expectancy increases continuously with income, with no evidence of a threshold beyond which higher income ceases to be associated with longer life. More recently, BradyKohlerZheng2023Novel find that mortality risk increases with poverty exposure, from a 42% higher hazard under current poverty to a 71% higher hazard for those in poverty over the past 10 years. In what follows, we incorporate this prior information by imposing that the conditional mean of the outcome is increasing in the poverty rate and decreasing in median income.

\pgfplotsset{ gtatt panel/.style={ width=0.5\textwidth, height=0.3\textwidth, xmin=2008.5, xmax=2019.5, xtick={2009,2011,2013,2015,2017,2019}, extra y ticks={0}, extra y tick labels=, extra y tick style={grid=major, grid style={solid, gray!40}}, major tick length=1pt, tick label style={/pgf/number format/fixed,/pgf/number format/1000 sep=,font=\tiny}, yticklabel style={rotate=90}, title style={yshift=-0.65em}, legend style={draw=none,legend columns=-1,column sep=6pt}, mark size=0.68pt, line width=0.5pt, error bars/error mark options={solid,rotate=90,mark size=0.8pt}, extra x tick labels=, extra x tick style={grid=major, grid style={densely dashed}}, }, right axis/.style={ytick pos=right, yticklabel pos=right}, style cs/.style={color=Dark2-A,mark=10-pointed star}, style dml/.style={color=Dark2-B,mark=halfsquare left*,error bars/error bar style={densely dotted}}, style sdml/.style={color=Dark2-C,mark=halfcircle*,every mark/.append style={rotate=270}, error bars/error bar style={densely dashdotted}}, event panel/.style={ gtatt panel, xmin=-10.5, xmax=5.5, xtick={-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5}, extra x ticks={-1}, }, }

\pgfplotstableread{ group time event cs ci_low_cs ci_high_cs dml ci_low_dml ci_high_dml sdml ci_low_sdml ci_high_sdml 2014 2009 -5 3.948646 -0.928772 8.826064 3.523899 -0.090910 7.138708 2.303480 -1.297603 5.904563 2014 2010 -4 1.747725 -4.280025 7.775476 -1.608158 -4.278177 1.061861 -2.208195 -5.237202 0.820812 2014 2011 -3 2.954795 -1.574432 7.484022 1.176764 -1.304205 3.657732 -0.166376 -2.787810 2.455058 2014 2012 -2 2.848708 -3.738555 9.435972 1.972770 -0.130555 4.076095 1.090128 -1.534811 3.715067 2014 2013 -1 0.000000 -0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 2014 2014 0 -2.131197 -8.505874 4.243479 -1.582444 -3.627916 0.463027 -1.590870 -4.028102 0.846363 2014 2015 1 -3.647074 -11.953958 4.659811 -0.384803 -2.604925 1.835318 -0.176149 -3.117000 2.764703 2014 2016 2 -4.963733 -16.472488 6.545023 1.420331 -2.216208 5.056870 0.932510 -2.672616 4.537637 2014 2017 3 -6.064612 -17.166432 5.037207 1.722951 -1.629148 5.075050 1.534260 -2.279705 5.348226 2014 2018 4 -6.370834 -19.331034 6.589367 2.127473 -1.022951 5.277897 2.002689 -1.972682 5.978059 2014 2019 5 2.477249 -13.582551 18.537050 5.376410 1.826185 8.926635 7.665728 3.451323 11.880133 2015 2009 -6 4.920485 -3.160473 13.001443 3.885150 -4.864860 12.635160 2.348260 -6.376455 11.072974 2015 2010 -5 -0.183365 -7.437234 7.070505 -1.059362 -8.754017 6.635293 -0.799739 -9.054191 7.454714 2015 2011 -4 8.479854 -1.637637 18.597344 7.545236 -3.986218 19.076689 6.619476 -5.340010 18.578962 2015 2012 -3 3.459603 -3.716142 10.635348 2.316672 -5.106735 9.740079 0.644504 -8.021409 9.310416 2015 2013 -2 3.012421 -3.650499 9.675342 1.442521 -3.756346 6.641388 0.523516 -4.616835 5.663868 2015 2014 -1 0.000000 -0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 2015 2015 0 2.703258 -4.332511 9.739028 3.261534 -3.338786 9.861855 5.271275 -1.560370 12.102921 2015 2016 1 12.148210 5.284887 19.011533 12.374765 5.041894 19.707637 14.384918 7.330039 21.439796 2015 2017 2 17.677326 9.950412 25.404241 18.348153 9.385488 27.310818 18.711310 9.490248 27.932371 2015 2018 3 4.595611 -2.448458 11.639680 4.871064 -1.883436 11.625564 7.865058 0.558312 15.171804 2015 2019 4 7.016006 -0.709274 14.741287 6.148231 -0.733801 13.030263 10.657663 3.173123 18.142202 2016 2009 -7 14.048370 -9.286057 37.382796 1.214199 -14.551690 16.980088 1.152305 -13.893301 16.197911 2016 2010 -6 40.983577 -48.429000 130.396153 -7.623372 -20.089579 4.842835 -7.141116 -18.891828 4.609596 2016 2011 -5 9.289281 -37.210715 55.789277 -14.680419 -31.665132 2.304293 -14.095416 -29.237747 1.046916 2016 2012 -4 2.551677 -12.092603 17.195958 -5.603523 -20.466134 9.259087 -4.673502 -17.510532 8.163528 2016 2013 -3 5.143678 -8.033429 18.320785 0.817462 -11.277045 12.911970 1.415289 -10.454732 13.285309 2016 2014 -2 3.683277 -10.974922 18.341476 0.866514 -14.093323 15.826352 0.943132 -14.278174 16.164439 2016 2015 -1 0.000000 -0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 2016 2016 0 -8.697067 -23.608093 6.213958 -9.067751 -23.949157 5.813655 -9.088724 -24.113680 5.936232 2016 2017 1 -9.916357 -22.257083 2.424368 -8.862719 -20.500548 2.775111 -7.126508 -18.576569 4.323552 2016 2018 2 -18.922757 -31.700941 -6.144574 -17.601591 -29.699477 -5.503704 -14.166636 -26.137776 -2.195496 2016 2019 3 -13.986775 -27.535852 -0.437697 -10.306329 -22.918923 2.306265 -9.414144 -23.282549 4.454261 2019 2009 -10 -15.350323 -33.067353 2.366706 -4.483858 -17.087003 8.119286 -11.639466 -26.602817 3.323885 2019 2010 -9 -25.791807 -42.620676 -8.962938 -12.418756 -25.535209 0.697697 -15.168911 -29.609977 -0.727846 2019 2011 -8 -17.258425 -34.714323 0.197472 -8.090339 -21.380892 5.200214 -13.943306 -26.425235 -1.461376 2019 2012 -7 -12.467257 -28.174203 3.239689 -5.231264 -18.218487 7.755960 -15.405140 -29.171861 -1.638420 2019 2013 -6 -15.205226 -28.510416 -1.900036 -10.250854 -21.778406 1.276699 -16.603800 -27.936032 -5.271568 2019 2014 -5 -13.879145 -28.047122 0.288833 -9.493391 -20.655045 1.668263 -12.877358 -23.461946 -2.292770 2019 2015 -4 -5.506914 -16.912073 5.898245 -3.143221 -12.861272 6.574829 -10.407823 -20.032382 -0.783263 2019 2016 -3 -1.100543 -10.292344 8.091258 -0.373400 -8.412398 7.665598 -5.264003 -13.615630 3.087624 2019 2017 -2 3.201367 -5.774498 12.177232 4.199000 -3.632495 12.030495 0.769034 -6.930262 8.468329 2019 2018 -1 0.000000 -0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 2019 2019 0 3.308994 -5.763475 12.381463 1.806603 -6.992260 10.605467 2.671930 -6.729099 12.072958 }\GTATT

\pgfplotstableread{ event cs ci_low_cs ci_high_cs dml ci_low_dml ci_high_dml sdml ci_low_sdml ci_high_sdml -10 -15.350323 -33.067353 2.366706 -4.483858 -17.084838 8.117121 -11.639466 -26.600880 3.321948 -9 -25.791807 -42.620676 -8.962938 -12.418756 -25.533155 0.695643 -15.168911 -29.607958 -0.729864 -8 -17.258425 -34.714323 0.197472 -8.090339 -21.378802 5.198125 -13.943306 -26.423059 -1.463552 -7 -2.808912 -15.750252 10.132429 -3.750111 -13.716750 6.216527 -8.630821 -19.539082 2.277441 -6 5.068632 -10.275735 20.412999 -2.187134 -9.047862 4.673593 -4.448304 -11.360484 2.463876 -5 2.668480 -1.679572 7.016532 1.765245 -1.474290 5.004780 0.628801 -2.693216 3.950818 -4 2.133383 -2.924812 7.191578 -0.999102 -3.583249 1.585044 -1.703218 -4.719892 1.313455 -3 2.857439 -0.934498 6.649376 1.132792 -1.003286 3.268871 -0.060983 -2.338991 2.217026 -2 2.912967 -2.394760 8.220695 1.924381 0.019453 3.829310 1.071209 -1.197136 3.339554 -1 0.000000 -0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0 -1.492955 -6.716460 3.730549 -1.077010 -3.004868 0.850848 -0.524763 -2.715184 1.665658 1 -1.969392 -9.001513 5.062729 0.857938 -1.353896 3.069771 1.359657 -1.503312 4.222627 2 -2.725066 -12.488982 7.038849 2.859542 -0.656536 6.375619 2.533242 -0.942054 6.008538 3 -5.055663 -14.535884 4.424558 1.723326 -1.221129 4.667781 1.730479 -1.605308 5.066267 4 -4.716163 -16.157557 6.725231 2.725095 -0.174855 5.625045 3.072480 -0.559026 6.703987 5 2.477249 -13.582551 18.537050 5.376410 1.826298 8.926522 7.665728 3.451455 11.880001 }\EventATT

We implement three estimation methods: the Callaway-Sant'Anna (parametric) method employed in BakerCallawayCunninghamGoodmanSantAnna2025Guide as the benchmark, debiased machine learning (via Riesz regression) without shape constraints, and our shape constrained debiased machine learning, respectively labeled as CS, DML, and SDML. The implementation details of DML and SDML are aligned with those in Section (ref). For simplicity, we only report pointwise (in time) confidence intervals at the 5% significance level.

Figures (ref) and (ref) present GT-ATTs over calendar time and event time respectively for each expansion group. For the 2014 cohort, the CS estimates are negative in most post-expansion years, whereas the DML and SDML estimates become positive after 2015 and are significantly positive by 2019. For the 2015 cohort, all three methods produce similar positive effects during the first two post-expansion years, but the SDML estimates are systematically larger and remain significant in 2018 and 2019, when the CS and DML confidence intervals include zero. For the 2016 cohort, all three methods produce negative estimates that are significantly negative in 2018; by 2019, only the CS estimate remains significant. The contemporaneous estimates for the 2019 cohort are similar and insignificant across methods, but their pre-treatment estimates differ substantially, with SDML generally producing more negative estimates than DML. These differences carry over to Figure (ref): the CS estimates are negative from event time zero through four, whereas the DML and SDML estimates are positive from event time one onward, although all corresponding confidence intervals include zero. The aggregate pre-treatment estimates are not uniformly close to zero: CS and SDML are significantly negative at event time $-9$, SDML remains significantly negative at event time $-8$, and DML is narrowly positive and significant at event time $-2$. At event time five, the DML and SDML estimates are positive and significant, while the CS estimate is smaller and considerably less precise. Overall, relative to DML, imposing shape constraints produces meaningful differences in the magnitude and statistical significance of the estimated effects across several cohorts and time periods.

Working Hours and Wage Growth

How working hours translate into wage growth is central to understanding career progression, labor supply incentives, and the accumulation of earnings differences across workers. The theoretical case for allowing this relationship to be nonlinear goes back at least to Barzel1973Wage, who treats hours and wages as jointly determined components of an employment arrangement rather than separate objects linked by a fixed hourly wage. Consistent with this perspective, subsequent empirical work finds that hourly earnings vary systematically with hours worked, pointing to nonlinear wage-hours relationships BiddleZarkin1989Wage,BickBlandinRogerson2022Wage. Building on the dynamic complete information framework of GibbonsWaldman1999Wage, Gicheva2013Working develops a model that predicts a convex intertemporal relationship between working hours and wage growth, which is then informally validated based on a data set from the panel survey of registrants in the United States for the Graduate Management Admission Test (GMAT) between June 1990 and March 1991---see Fang2021Unifying for a formal test.

To map this problem into our framework, we adopt a flexible nonparametric model (instead of the partially linear model in Gicheva2013Working):

where $Y$ is the annual wage growth rate, $D$ is weekly working hours, $Z$ is a vector of demographic control variables (e.g., experience, education, gender, age, race, and family characteristics), and $E[U|D,Z]=0$. In the GMAT sample, there are $19$ control variables including experience, education, gender, age, race, and family characteristics. The parameter of interest is the weighted average derivative:

where $w$ is a weight function.

In the implementation below, we employ the same GMAT dataset as Gicheva2013Working and set $w$ to be the normalized indicator of the empirical $[0.1,0.9]$ quantile region of $D$. The first step $\gamma_0$ is estimated by a two-branch neural network subject to convexity with respect to $D$. The Riesz representer $\alpha_0$ is characterized as the solution to

where $\Gamma^\smallsmile\subset L^2(D,Z)$ is the class of functions $(d,z)\mapsto\gamma(d,z)$ convex in $d$. We keep the tuning of architectures consistent with the simulations for convexity in Section (ref). The final cleaned sample consists of $1,911$ survey respondents. We also consider men and women separately as in Gicheva2013Working.

Table (ref) reports estimates of the weighted average derivative. The unconstrained DML estimates are positive and significant at the $5\%$ level for the full sample and for men, while the estimate for women is smaller and statistically indistinguishable from zero. The SDML point estimates are also positive for the full sample and for men but are estimated less precisely, and the SDML estimate for women is slightly negative; all three SDML confidence intervals include zero. Overall, the table suggests a modest average marginal association, if any, between weekly working hours and subsequent wage growth. This is consistent with the broader wage-hours literature, which emphasizes that the relationship between hours and earnings is nonlinear and heterogeneous rather than uniformly positive at the margin. For example, DenningJacobLefgrenLehn2022Wage show that marginal returns to hours can be small within labor-market settings, even when hours are important for broader wage differences.

Conclusion

In this paper, we present a general framework of identification and estimation for DML where the parameter of interest $\theta_0$ is identified by a GMM model in the presence of a first step nuisance $\gamma_0$. The parameter $\gamma_0$ may be a high dimensional parameter defined by a model with endogeneity and possibly subject to nonlinear shape constraints. In particular, we establish identification of the Riesz representer $\alpha_0$ under mild conditions and generalize the Riesz regression to accommodate a generic first step $\gamma_0$. Since machine learning alone is incapable of overcoming the curse of dimensionality, we incorporate shape constraints on $\gamma_0$ as a vehicle for bringing economic structures into DML. By exploiting linearity of the first step influence function with respect to $\alpha_0$, we are able to develop a theory that is conceptually simple, transparent, and unifying. We believe our framework benefits statistical inference on many causal and structural problems, in both high dimensional and classical semiparametric settings.