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.
104,587 characters · 20 sections · 79 citation commands
Generalized Bayes in Conditional Moment Restriction Models
Conditional moment restrictions are widely used to identify structural parameters in complex economic models. In many applications, the object of interest is an unknown nonparametric structural function $h_0(\cdot)$ that satisfies
where $Y \in \mathbb{R}^{d_y}$ is a vector of outcomes, $X \in \mathbb{R}^d$ is a vector of endogenous regressors, $W \in \mathbb{R}^{d_w}$ is a vector of conditioning (or instrumental) variables, and the conditional distribution of $(Y,X) \mid W$ is left unrestricted. Here, $\rho(.) = [ \rho_{1}(.) , \dots , \rho_{d_{\rho}}(.) ] $ is a $d_{\rho}$ dimensional vector of generalized residual functions, whose functional forms are assumed to be fully known. Common applications of this framework include consumer demand \citep*{blundell2007semi}, firm productivity (doraszelski2013r), differentiated product markets berry2024nonparametric, production functions (\citealp*{ackerberg2015identification}), international trade \citep*{adao2017nonparametric}, treatment effects chernozhukov2005iv and asset pricing bansal1993no,chen2009land.
A common challenge for practitioners is that, although these restrictions are informative in the population, their finite-sample information content can be quite limited. In parametric models, this issue is typically attributed to weak instruments \citep*{stock2002survey}, whereas in nonparametric endogenous settings it reflects an “ill-posed inverse” problem \citep*{chen2012estimation}. As a result, classical nonparametric estimators often display undesirable properties such as high finite-sample variability, irregular behavior, and extreme sensitivity to small data perturbations. These difficulties are particularly evident in applications with multivariate endogenous regressors or when closed-form solutions are unavailable.
Motivated by these concerns, this paper proposes a class of nonparametric estimators and confidence sets obtained as solutions to generalized (quasi-) Bayes decision rules. In this framework, the conditional restrictions are interpreted as a quasi-likelihood which, when combined with a prior, yields a generalized Bayesian nonlinear inverse problem for the structural parameter. To fix ideas, let $\widehat{m}(\cdot)$ denote a feasible first-stage estimator of $m(W,h) = \mathbb{E}[\rho(Y,h)\mid W]$, $\widehat{\Sigma}(\cdot)$ a positive semi-definite weighting matrix, and $d\mu(\cdot)$ a prior on structural functions. We then study the generalized Bayes posterior distribution: \[ \mu(\cdot \mid \mathcal{D}_n) \;=\; \frac{\exp\!\big(-\frac{n}{2}\,\mathbb{E}_n\!\big[\widehat{m}(W,\cdot)^{\prime}\,\widehat{\Sigma}(W)\,\widehat{m}(W,\cdot)\big]\big)\, d\mu(\cdot)} {\int \exp\!\big(-\frac{n}{2}\,\mathbb{E}_n\!\big[\widehat{m}(W,h)^{\prime}\,\widehat{\Sigma}(W)\,\widehat{m}(W,h)\big]\big)\, d\mu(h)} \, . \] In the nonparametric endogenous models considered here, this framework provides a powerful form of data-driven regularization. Importantly, it also allows researchers to incorporate auxiliary information that strengthens the finite sample information content of the moments. Such information may range from weakly informative features, such as smoothness, to restrictions informed by application-specific microfoundations.
Over the past two decades, parametric quasi-Bayes procedures have found a variety of applications in econometrics, from models with nonsmooth objectives chernozhukov2005iv to settings with nonstandard identification (\citealp*{chen2018monte}; andrews2022optimal). Most of the literature has focused on the properties of quasi-posteriors in parametric models. By contrast, relatively little is known about the behavior of quasi-Bayes in settings with a nonparametric structural parameter. This article helps bridge that gap by providing a unified treatment of quasi-Bayes for the broad class of nonparametric conditional moment restriction models commonly encountered in applied work. As we illustrate, when paired with a suitable nonparametric prior, quasi-Bayes naturally functions as a powerful form of data-driven regularization in endogenous models.
The main theoretical contributions of this paper are as follows. First, we introduce a theoretically motivated class of Gaussian process priors to model the nonparametric structural parameter. Together with the conditional restrictions, this induces a generalized (quasi-) Bayes posterior for the parameter. Second, we derive posterior contraction rates for the quasi-Bayes posterior in classical $L^2$ metrics. Third, we establish conditions under which a nonparametric Bernstein–von Mises (BvM) theorem holds for the quasi-Bayes posterior. We use this to provide frequentist guarantees for certain optimally weighted quasi-Bayes credible sets that are centered around the posterior mean. In particular, we show that such credible sets achieve asymptotically exact frequentist coverage. This provides the first nonparametric quasi-Bayes inferential guarantee in the literature, extending classical results (e.g. chernozhukov2003mcmc) for parametric GMM models.
We demonstrate the viability of our procedures across a broad class of models, including classical linear nonparametric IV, conditional quantile restrictions, and general nonlinear conditional restrictions. We complement this with extensive simulation evidence, replicating all univariate benchmark designs from the literature and extending them to settings with multivariate endogenous regressors. To highlight the flexibility of our approach, we additionally estimate models under alternative sets of restrictions whenever such alternatives are available. Overall, we expect our generalized Bayes procedures and accompanying implementation toolkit to be broadly useful for nonlinear conditional moment restrictions, particularly in ill-posed problems or when closed-form solutions are unavailable.
The paper is organized as follows. Section (ref) introduces the class of conditional moment restriction models and develops the generalized (quasi-) Bayes framework. Section (ref) discusses our motivation for generalized Bayes procedures and relates it to the broader econometric literature. Section (ref) presents the assumptions and develops the main results. Sections (ref) and (ref) provide simulation evidence on the performance of generalized Bayes estimators relative to common alternatives. In Section (ref), we apply our methodology to nonlinear restrictions that arise in the nonparametric estimation of production functions. Section (ref) provides additional remarks and concludes. Appendices (ref), (ref), (ref), and (ref) provide additional details on simulations, implementation, theory, and proofs, respectively.
There is a large literature on nonparametric sieve-based frequentist estimation and inference for conditional moment restriction models. As part of our general analysis, we review a subset of this literature in Sections (ref)–(ref). For a more comprehensive survey, particularly on early contributions, see chen2016methods.
In econometrics, our work is most closely related to chen2012estimation, chen2015sieve, who developed the foundational frequentist sieve-based analysis of general conditional moment restriction models. At a high level, our procedures provide a generalized Bayes counterpart to their theory for infinite-dimensional sieves. However, instead of relying on traditional sieves and penalization, we develop procedures that are built around a class of infinite dimensional Gaussian process priors.
chernozhukov2003mcmc developed the quasi-Bayes limit theory for parametric models strongly identified by a collection of moments. For finite-dimensional structural parameters, several alternative approaches have been proposed, including exponentially tilted empirical likelihoods \citep*{schennach2005bayesian,chib2018bayesian, chib2022bayesian} and methods that project a posterior on the data-generating distribution onto the parameter of interest chamberlain2003nonparametric, walker2024semiparametric. By contrast, our focus is on endogenous models in which the parameters of interest are nonparametric structural functions. Importantly, in this setting, the structural parameter is infinite-dimensional, and its recovery is a challenging statistical ill-posed inverse problem.
In the statistical literature, early extensions of chernozhukov2003mcmc to nonparametric models focused on slowly growing uninformative flat sieve priors. This line of work includes conditions for basic consistency liao2011posterior and convergence rates in the special case of linear nonparametric IV models kato2013quasi. These approaches parallel classical frequentist analysis (e.g. ai2003efficient,newey2003instrumental), where regularization is achieved by restricting estimation to a sequence of slowly expanding sieve spaces. By contrast, we study generalized Bayes procedures with infinite dimensional Gaussian process priors and develop statistical guarantees for general nonlinear conditional moment restrictions.
As we illustrate in Sections (ref) and (ref), the regularizing properties of the Gaussian process priors we study make them particularly well-suited to nonparametric endogenous models identified via general conditional moment restrictions. This motivation connects to early econometric work on the consistency of Gaussian priors in conjugate linear models with a known operator \citep*{florens2012nonparametric}.\footnote{For related work in statistics, see also knapik2011bayesian, gugushvili2020bayesian.} Our setting allows for general nonlinear and possibly nonsmooth restrictions with an unknown operator, leading to a non-conjugate quasi-Bayes posterior based on an estimated first-stage likelihood. Addressing this general case is necessary to cover the wide range of conditional moment restrictions commonly encountered in applied work, and our analysis develops both estimation and inferential guarantees in this setting.
Finally, in the special case of regression with exogenous covariates, our procedures relate to a growing literature in applied mathematics that examines Gaussian priors for nonlinear regression models with homoscedastic Gaussian noise dashti2015bayesian,monard2021statistical, nickl2023bayesian. Our framework can be seen as complementary to this line of work, providing a generalized Bayes analogue that accomodates certain forms of heteroskedasticity and non-Gaussianity.
Let $(Y,X,W)$ denote random vectors, where $Y \in \mathbb{R}^{d_y}$ is the outcome, $X \in \mathbb{R}^{d}$ the regressors, and $W \in \mathbb{R}^{d_w}$ the conditioning (instrumental) variables. We are interested in an unknown structural function $h_0$ that satisfies the conditional moment restriction
Here, $\rho(.) = [ \rho_{1}(.) , \dots , \rho_{d_{\rho}}(.) ] $ is a $d_{\rho}$ dimensional vector of generalized residual functions, whose functional forms are assumed to be fully known. Components of $X$ that are exogenous may, without loss of generality, be included in $W$. As is standard in applications, the conditional distribution of $(Y,X)$ given $W$ is not assumed to be known.
This framework is very general. By varying the choice of $\rho(\cdot)$, we can recover a large class of structural models commonly encountered in applied work. The form of the conditional restrictions, or equivalently the choice of generalized residual $\rho(\cdot)$, typically varies significantly across applications. The following examples illustrate some of these restrictions in further detail.
For intuition and as a guide to our general analysis, we will frequently refer to Examples (ref) and (ref). We view these two examples as useful benchmark models for the following reason. In Example (ref), the residual \( \rho(.) \) is a smooth linear function of $h$, whereas in Example (ref), it is highly nonlinear and nonsmooth in $h$. In particular, they exemplify two distinct classes of models, distinguished by the regularity of the residual function. Although the restrictions encountered in empirical applications often appear more complex, their analysis and limiting structure can typically be characterized between these two extremes.
Given a function \( h(X) \), we denote the conditional mean of the generalized residual by $$ m(W,h) = \mathbb{E} \big[ \rho(Y,h(X)) \mid W \big].$$ The restriction \(m(W,h_0)=\mathbf{0}\) implies that \(h_0\) is the minimizer of the population criterion \[ Q(h)=\mathbb{E}\!\left[m(W,h)^{\prime}\,\Sigma(W)\,m(W,h)\right], \] where \(\Sigma(W)\in\mathbb{R}^{d_{\rho}\times d_{\rho}}\) is a positive-definite weighting matrix.
As the distributional structure of the data is not assumed to be known, working with \( Q(h) \) directly is infeasible. The standard approach (e.g. ai2003efficient,newey2003instrumental,chen2012estimation) replaces \( m(W,h) \) and $\Sigma (\cdot)$ with suitable empirical analogs. Specifically, let \( \widehat{m}(W,h) \) and \( \widehat{\Sigma}(W) \) be “first-stage” estimators of \( m(W,h) \) and \( \Sigma(W) \), respectively. Then, a feasible finite-sample objective function is
The classical approach to estimating \( h_0 \) involves a “second stage”, where \( Q_n(\cdot) \) is minimized over a suitable parameter space $\mathcal{H}_n$ to obtain an estimator \( \widehat{h} \). As noted in the literature (e.g. \citealp*{chetverikov2017nonparametric}), these solutions often exhibit substantial finite-sample variability and are highly sensitive to small perturbations in the data and user-selected tuning parameters such as the complexity of $\mathcal{H}_n$. Intuitively, the second stage is “ill-posed” and the large finite-sample variability of these estimators arises from their representation as the inverse of an ill-posed objective.
To stabilize the inverse problem and more efficiently utilize the information content in the conditional moments, we examines a class of nonparametric estimators that arise as solutions to generalized Bayes decision rules. Specifically, we view the conditional moment restriction as a nonlinear inverse problem for the infinite dimensional structural parameter $h_0$. The restriction $ m(W,h_0) = \mathbf{0} $ then motivates a quasi-Bayes likelihood of the form
Denote the observed data by $\mathcal{D}_n = \{ (X_1,Y_1,W_1), \dots , (X_n,Y_n,W_n) \}$. By combining the likelihood $L(.)$ with a (possibly data dependent) prior $\mu$ over structural functions, we obtain the generalized (quasi-) Bayes posterior:
Related to this construction, liao2011posterior transformed the conditional moment restrictions into a growing set of unconditional moments and proved the asymptotic consistency of a classical quasi-Bayes GMM criterion chernozhukov2003mcmc under slowly growing flat sieve priors. In contrast, we follow the conventional frequentist approach, in which the first-stage functional $\widehat{m}(\cdot)$ is estimated directly, and we then treat the objective function $L(\cdot)$ in ((ref)) as a quasi-likelihood for the model.
In this paper, we focus on a class of infinite dimensional Gaussian process priors for $d \mu (\cdot)$. When the structural function $h_0(\cdot)$ is defined over a bounded smooth domain $\mathcal{X} \subset \mathbb{R}^d$, a common choice is the family of Whittle--Matérn Gaussian process priors williams2006gaussian.
Gaussian process priors are widely employed in Bayesian nonlinear inverse problems, especially in applications arising within applied mathematics nickl2023bayesian. To fix ideas, consider a mean-zero Gaussian process \( G \) with realizations in a Hilbert space \( \mathcal{H} \) and covariance operator \( \Lambda \). By the spectral theorem, there exists an orthonormal basis of eigenfunctions \( (e_i)_{i=1}^{\infty} \subset \mathcal{H} \) that diagonalizes \( \Lambda \). If \( \lambda_i \) denotes the non-negative eigenvalue associated with \( e_i \), then \( G \) admits a unique Karhunen-Loève expansion of the form:
Intuitively, the rate at which \( \lambda_i \to 0 \) serves as a measure of the process's smoothness relative to the eigenbasis. If \( (e_i)_{i=1}^{\infty} \) denotes the standard Fourier basis, this corresponds to classical Sobolev smoothness.
Similar to the analysis in knapik2011bayesian, we consider a family of Gaussian process priors $\{G_{\alpha} : \alpha \in \mathcal{L} \}$ that are indexed by a regularity hyperparameter $\alpha \in \mathcal{L} \subset \mathbb{R}_+$. In this setting, each process $G_{\alpha}$ admits an expansion of the form\footnote{If the mapping $\alpha \mapsto \lambda_{i,\alpha}$ influences the exponent in a different way, the results can also be stated in terms of the induced exponent $s(\alpha)$, i.e., $\lambda_{i,\alpha} \asymp i^{-s(\alpha)}$.}
where $\lambda_{i,\alpha} \asymp i^{-(1 + 2\alpha/d)}$ and $(e_i)_{i=1}^{\infty}$ is an orthonormal basis of $L^2(\mathcal{X})$.
While we do not impose any restrictions on the eigenbasis $(e_i)_{i=1}^{\infty}$ directly, we will typically require the sample paths of the Gaussian process $G_{\alpha}$ (for $\alpha \in \mathcal{L}$) to satisfy some minimum regularity (see Condition (ref) below). In most cases, this can be satisfied by restricting the regularity index set to $ \alpha \in \mathcal{L} \subseteq [\underline{\alpha}, \infty)$ for some minimum regularity $\underline{\alpha} > 0$. The following example illustrates the general idea for a widely used family of Gaussian process priors.
Intuitively, larger values of $\alpha$ correspond to smoother sample paths. In certain applications, suitable smoothness levels can be motivated by prior studies or application-specific microfoundations. In settings where such guidance is unavailable, $\alpha = 3/2$ and $\alpha = 5/2$ are widely used as standard defaults williams2006gaussian, offering a balance between regularity and flexibility to accommodate irregular variation.
Researchers have considerable flexibility in the choice of the first-stage estimator for the conditional mean $ m(W,h) = \mathbb{E}[\rho(Y,h(X)) \mid W] $. This can accomodate a broad range of regression and machine learning methods. In practice, however, it will be convenient to focus on estimators that are computationally efficient, as this ensures that the quasi-likelihood $L ( \cdot) $ in ((ref)) can be evaluated efficiently.
A common and efficient choice is to consider sieve-based first stages, defined as linear projections onto a set of basis functions. Let $b^K(W) = [ b_1(W), \dots , b_K(W) ]'$ denote a vector of first stage approximating functions. Then, for a given function $h(X)$, we estimate the conditional mean by the least squares projection:
In low dimensions, approximating functions can be formed from tensor products of standard univariate bases (e.g. Fourier series, splines), eigenfunction expansions and indicator functions to accommodate discrete instruments. In higher dimensions, common alternatives are bases constructed using randomized features (e.g. rahimi2007random).
To facilitate detailed analysis and clarity of exposition, we focus on a classical first stage defined by a linear projection onto approximating functions.\footnote{In Appendix (ref), we provide some theory for contraction with generic first-stage estimators.} Although our main results extend to other first-stage estimators, the conditions required to obtain statistical guarantees will generally depend on the specific choice of estimator. By concentrating on the sieve case, we keep the first-stage analysis self-contained and directly comparable to the classical frequentist analysis of conditional moment restriction models.
In the classical frequentist literature (e.g., \citealp*{blundell2007semi, chen2012estimation}), the choice of first stage estimator is typically not viewed as a “key tuning parameter.” Intuitively, estimating the smooth conditional mean $ \mathbb{E}[\rho(Y,h(X)) \mid W] $ is a well-posed regression problem and is far less sensitive to tuning than a classical ill-posed inverse problem. This is also true in our setting. Specifically, if $\Theta_n$ denotes a suitable collection of high probability regular sample paths of the Gaussian process, the first stage is best viewed as providing an efficient approximation to the conditional mean operator $ \Theta_n \ni h \mapsto \mathbb{E}[\rho(Y,h(X)) \mid W] $.
In this section, we discuss the econometric and practical motivation for quasi-Bayes procedures, with emphasis on their application to nonparametric endogenous models. We begin with the econometric motivation, particularly in comparison with fully Bayesian and classical frequentist approaches.
A fully Bayes approach to this problem would typically require explicit modeling of the conditional distribution $(Y,X)\mid W$. Since our primary object of interest is the structural parameter, this distribution is a complex nuisance, and modeling it may be undesirable in many settings. Analogous to the econometric motivation underlying classical GMM hansen1982generalized, it is often preferable to target the structural parameter directly, particularly when the parameter itself is a complex nonparametric object.\footnote{For finite dimensional structural parameters, a similar point was made by chernozhukov2003mcmc.}
Beyond modeling challenges, the analysis in \citet*{bornn2019moment,florens2021gaussian} also highlight that, even with parametric structural parameters, there are subtle probabilistic difficulties in specifying a joint prior on the nuisance law $F_{(Y,X)\mid W}$ and structural parameter.\footnote{Constructing a reasonable prior on the low dimensional manifold $\Theta=\{(h,F): \mathbb{E}_{F}[\rho(Y,h(X)) \mid W]=0\}$ is challenging: for any fixed $h$, classical priors typically assign probability zero to the fiber $\mathcal{F}_h=\{F:\mathbb{E}_{F}[\rho(Y,h(X)) \mid W]=0\}$. This difficulty arises even in simpler settings with unconditional moments and finite-dimensional structural parameters.} In our setting with an infinite dimensional structural function, this becomes considerably more challenging. Although it may be possible, in theory, to proceed without a prior on the structural function, this is ill-advised for the nonparametric endogenous models we study, as it forgoes the regularization, interpretability, and flexibility gained by placing the prior directly on the structural function.
From the preceding discussion, it follows that quasi-Bayes can be viewed as a convenient hybrid between frequentist and fully Bayes methods. Similar to classical frequentist procedures, it utilizes the efficient first stage as a sufficient statistic for the nuisance. In the second stage, the difficult, ill-posed recovery of the structural function is formulated as a generalized Bayesian nonlinear inverse problem nickl2023bayesian. In this setting, the prior on the structural function provides a powerful form of data driven regularization, while also allowing the researcher to incorporate domain-specific knowledge.
To illustrate some of our motivation in greater detail, we make use of all the benchmark designs previously employed in the nonparametric instrumental variable (NPIV) literature. Specifically, we consider the designs from newey2003instrumental, santos2012inference, \citet*{chernozhukov2015constrained}, chetverikov2017nonparametric, and \citet*{chen2025adaptive}, which we refer to as NP, S, CNS, CW and CCK, respectively. In all of these designs, the regressor is univariate and the structural function is estimated under a nonparametric instrumental variable (NPIV) restriction (Example (ref)). Details on all the designs are contained in Appendix (ref).
Let $\mathcal{D}_n$ denote the observed data, and let $X'$ be an independent draw from the distribution of $X$. Given an estimator $\widehat{h} = \widehat{h}(\mathcal{D}_n)$, we define the expected out-of-sample root mean squared risk: \[ \mathcal{R}(\widehat{h}, h_0) = \left\{ \mathbb{E}_{\mathcal{D}_n,\,X'} \left[ \big( \widehat{h}(X') - h_0(X') \big)^2 \right] \right\}^{1/2}. \] Let 2SLS denote the two-stage least squares estimator, where the first stage uses thin-plate splines and the structural function uses natural splines, both of dimension $J$.\footnote{Natural splines provide some regularization by enforcing $h''(x) = 0$ at the data boundary, implying linearity beyond. For larger $J$, results appeared more unstable with alternative bases.}
As Table (ref) illustrates, in endogenous models, classical estimators are highly sensitive to tuning parameters that determine the complexity of the parameter space. In some univariate settings (e.g. NPIV, \citealp*{chen2025adaptive}), this complexity can be tuned in a data driven way. However, in models with generalized nonlinear restrictions, multivariate regressors, or no closed-form solutions, effective tuning becomes substantially more challenging. Indeed, to the best of our knowledge, no regularization mechanism has yet been demonstrated to perform successfully across the broad range of models, restrictions and data generating processes encountered in theoretical and empirical work.
It is well known that Bayes procedures regularize naturally via the prior, albeit at the cost of potential finite-sample bias. In endogenous settings, the resulting variance reduction can be substantial. In nonparametric Bayes procedures, this bias typically takes the form of a preference for well-behaved or regular functions. We argue that this property is particularly valuable as a regularization mechanism in nonparametric endogenous models, where structural function regularity is typically already a prerequisite for any meaningful analysis. Indeed, this feature is evident in all the designs reported in Table (ref) and all other designs considered in the broader literature.
To further illustrate the preceding point, consider all the designs in Table (ref). They can be estimated using either of the following generalized residuals: \[
\]
In general, the NPQIV restriction is considered more challenging, as it involves a nonlinear and nonsmooth residual. Let QB denote the quasi-Bayes posterior mean, based on a first-stage thin plate spline of dimension $K$ and a classical Whittle–Matérn Gaussian process prior. We use the same prior and implementation algorithm across all designs and both sets of restrictions. Further details are provided in Appendix (ref).
Table (ref) reports the quasi-Bayes risk for all designs in Table (ref), under both NPIV and NPQIV restrictions. The estimates appear remarkably accurate and stable across both restrictions. A natural question is how far these findings extend. For example, can they generalize to more challenging settings with multivariate regressors? In Section (ref), we provide additional evidence by examining multivariate extensions of the designs in Table (ref).
As a final remark, we note that these procedures differ from classical frequentist regularization in two key ways. First, as noted earlier, devising a broadly effective data-driven regularization scheme that works across all models and restrictions is highly challenging. By contrast, in our quasi-Bayes framework, the priors we employ induce a nontrivial form of regularization that has proven effective in a wide range of applications, particularly in nonlinear inverse problems.\footnote{See ghosal2017fundamentals,nickl2023bayesian for an overview of applications.} Second, quasi-Bayes procedures are inherently data-driven through the interplay between the prior and the information in the conditional moments. This interaction is precisely what allows the information content in the moments to dominate in settings with strong identification and enables a single prior specification to yield reasonable results across all the designs and restrictions in Table (ref).
In this section, we develop the limit theory for the generalized (quasi-) Bayes posterior in ((ref)). Specifically, we examine the following questions in detail: (i) What are the minimal conditions on the model and prior that ensure quasi-Bayes consistency? (ii) How do convergence rates depend on the smoothness of the structural function \( h_0 \)? (iii) When do nonparametric quasi-Bayes credible sets achieve exact frequentist coverage?
To begin with, we state our main conditions on the generalized residual function $\rho(\cdot)$ that defines the conditional moment restriction in ((ref)). We assume that the endogenous regressor \(X\) is supported on a smooth bounded domain \(\mathcal X\subset\mathbb{R}^{d}\), and the instrument \(W\) is supported on a domain \(\mathcal W\subseteq\mathbb{R}^{d_w}\). This is standard in the literature and, if necessary, can always be satisfied by applying an appropriate transformation of the regressors.\footnote{In practice, apart from basic standardization, no transformations are used in our implementation.}
For any $t > 0$, let $(\mathbf{H}^t, \| \cdot \|_{\mathbf{H}^t})$ denote the usual Sobolev space of order $t$ over $\mathcal{X}$. The Sobolev ball of radius $M$ is denoted by $ \mathbf{H}^t(M) = \{ h : \| h \|_{\mathbf{H}^t} \leq M \}$.
\begin{condition2}[Local $L^2$ continuity] For some $ \kappa \in (0,1]$, $ t > d/ 2\kappa $ and any $M < \infty$, there exists $C_1 = C_1(M) < \infty$ such that
holds for all $h' \in \mathbf{H}^{t}(M)$ and $\xi > 0$ small enough. \end{condition2} In Condition (ref), the two expectations differ in the metrics they employ. The first expectation is over the the supremum with respect to the stronger $\| \cdot \|_{\infty}$ norm, whereas the outer supremum of the second expectation is taken under the weaker $ \|\cdot \|_{L^2(\mathbb{P})}$ norm. Intuitively, because the expected supremum is more difficult to control, it is taken over functions that are closer in a stronger metric.
Condition (ref) is analogous to conditions that are frequently used in the analysis of non-smooth objectives \citep*{chen2003estimation}. In particular, it permits a pointwise discontinuous residual function (e.g. NPQIV models) provided that $\rho(\cdot)$ is suitably uniformly continuous in $L^{2}(\mathbb{P})$ expectation. The parameter $\kappa$ is typically referred to as the local continuity exponent. It holds with $\kappa=1$ for the NPIV model (Example (ref)) and $\kappa=1/2$ for the NPQIV model (Example (ref)). \begin{condition2}[Residual moments] There exists $\epsilon, \delta > 0 $ and $t > d/2\kappa$ such that for any $M > 0 $, there exists finite constants $C_2(M) , C_3(M), C_4(M) < \infty$ that satisfy
\end{condition2} Condition (ref) imposes mild moment restrictions on the residual function: the bounds only need to hold over any fixed Sobolev ball. The assumption is trivially satisfied with bounded residual functions (e.g. NPQIV). More generally, if $t>d/2$, the Sobolev embedding theorem evans2022partial implies that $\mathbf H^{t}$ embeds continuously into a Hölder space, so functions in $\mathbf H^{t}(M)$ are uniformly bounded in the $\|\cdot\|_{\infty}$ norm. In most settings, this observation makes it straightforward to verify Condition (ref). For example, in the NPIV model, Condition (ref) holds if the unobserved error $u$ satisfies $\mathbb E\big(|u|^{2+\epsilon}\big)<\infty$ and $\mathbb E[u^{2}\big| W]\le\bar{\sigma}^{2}$ for some $\bar{\sigma}^{2} < \infty$.
While the generalized residual may be non-smooth, we assume (as is standard) that its smoothed conditional mean $m(W,h) = \mathbb{E}[\rho(Y,h(X)) \mid W]$ is sufficiently regular, in the sense that it satisfies a local Lipschitz property. This is formalized below in Condition (ref). \begin{condition2}[Locally Lipschitz conditional mean] For some $t > d/(2\kappa)$, the map $h \mapsto m(W,h)$ from $(\mathbf{H}^t, \| \cdot \|_{L^2(\mathcal{X})})$ to $(L^2(W),\| \cdot \|_{L^2(\mathbb{P})})$ is continuous. Furthermore, for every $M > 0 $, there exists a constant $C_5(M) < \infty$ such that $\| m(W,h) - m(W,h_0) \|_{L^2(\mathbb{P})} \leq C_5 \| h - h_0 \|_{L^2(\mathcal{X})}$ for every $h \in \mathbf{H}^t(M)$. \end{condition2}
In this section, we establish the consistency of general quasi-Bayes posteriors arising from suitably rescaled Gaussian process priors. As discussed in Section (ref), we consider a classical first stage based on projecting onto a set of basis (approximating) functions $ b^K(W) = \big[b_1(W), \dots, b_K(W)\big].$ Denote by $\Pi_K(\cdot)$, the $L^2(\mathbb{P})$ projection operator onto the span of these functions.
Following the discussion in Section (ref), let $G_{\alpha}$ denote a mean-zero Gaussian process with regularity parameter $\alpha > 0$. Let $(e_i)_{i=1}^{\infty}$ be the orthonormal eigenfunction basis of its covariance operator $\Lambda_{\alpha}$. Similar to the analysis in \citet*{knapik2011bayesian}, it will be convenient to measure regularity directly with respect to this basis.\footnote{When $G_{\alpha}$ is a Whittle--Matérn Gaussian process, or when $(e_i)_{i=1}^{\infty}$ is a standard Fourier basis, this reduces to classical Sobolev regularity.} To that end, for any $p > 0$, we define the associated $p$-regularity class as
Given $G_{\alpha}$ and first stage sieve dimension $K$, we consider the rescaled prior:
Rescaled Gaussian process priors are frequently employed in the analysis of Bayesian nonlinear inverse problems monard2021statistical,nickl2024posterior,nickl2025bayesian. In our conditional moment restriction framework, the scaling provides additional regularization that is crucial both for $(i)$ controlling the nonlinear ill-posedness of the inverse problem and $(ii)$ obtaining high-probability guarantees on the behavior of the first-stage estimator $\widehat{m}(W,h)$ used to approximate the conditional mean $h \mapsto m(W,h)$.
Intuitively, the posterior limit theory is determined by the interplay between the prior and the quasi-Bayes likelihood $ h \mapsto \mathbb{E}_n\!\big[ \widehat{m}(W,h)' \, \widehat{\Sigma}(W) \, \widehat{m}(W,h) \big] $. To formalize this interplay, we impose low-level conditions on three components: the prior, the weighting matrix $\widehat{\Sigma}(\cdot)$, and the first-stage basis functions $\{b_1(W), \dots, b_K(W)\}$ used to construct the conditional mean estimator $\widehat{m}(W,h)$. Our main requirements on these objects are summarized in the following two conditions.
\begin{condition2}[Regularity] $(i)$ The density of $X$ with respect to the Lebesgue measure is bounded away from $0$ and $\infty$ on $\mathcal{X}$. $(ii)$ $G_{\alpha}$ is a Gaussian random element on a separable subspace of the Sobolev space $\mathbf{H}^t$ for some $t > d/2\kappa$. \end{condition2} Condition (ref)$(i)$ is imposed for convenience, as it ensures the equivalence of the norms $\| \cdot \|_{L^2(\mathbb{P})}$ and $\| \cdot \|_{L^2(\mathcal{X})}$, where the latter is taken with respect to the Lebesgue measure. Condition (ref)$(ii)$ can be interpreted as a minimum regularity requirement in that it ensures the Gaussian process $G_{\alpha}$ has continuous and bounded sample paths.\footnote{This is a consequence of the Sobolev inequality evans2022partial, since $\mathbf{H}^t$ (for $t > d/2$) embeds into a H\"{o}lder space $C^{\beta}$ for some $\beta > 0$.}
\begin{condition2}[First stage approximation] $(i)$ The matrix $G_{b,K} = \mathbb{E} \big( [b^K(W)] [ b^K(W) ]' \big)$ is positive definite for all $K$ and $\zeta_{b,K} = \sup_{w \in \mathcal{W} } \| G_{b,K}^{-1/2} b^K(w) \|_{\ell^2} \lessapprox \sqrt{K}$. $(ii)$ The eigenvalues of $\widehat{\Sigma}(W)$ are asymptotically bounded above and below: $\mathbb{P}\big( c \leq \lambda_{\min}( \widehat{\Sigma}(W) ) \leq \lambda_{\max}( \widehat{\Sigma}(W) ) \leq C \big) \rightarrow 1$ for some $0 < c \leq C < \infty$. $(iii)$ For any fixed $ M > 0 $, the first stage is uniformly consistent over the Sobolev ball $\mathbf{H}^t(M)$: $ \sup_{h \in \mathbf{H}^t(M) } \|(\Pi_K - I) m(W,h) \|_{L^2(\mathbb{P})} \rightarrow 0 $ as $K \rightarrow \infty$. \end{condition2} Both Condition (ref)$(i)$, which restricts the growth of the $\| \cdot \|_{\ell^2}$ norm, and Condition (ref)$(iii)$, which requires uniform consistency over bounded regularity classes, are mild assumptions. They are satisfied by many standard bases, including splines, CDV wavelets, and Fourier series (see, e.g., chen2015optimal,belloni2015some).
Theorem (ref) establishes that the quasi-Bayes posterior is consistent provided that the regularity of the true function exceeds that of the Gaussian process by a factor of $d/2$. The upper bound constraint on $K_n$ is very mild: it guarantees that the first stage estimator $\widehat{m}(w,h)$ is well defined and uniformly approximates its population analog $\Pi_K m(w,h)$. By contrast, the theorem imposes a strict lower bound on the growth rate of the first-stage basis. Intuitively, larger values of $K_n$ increase sampling variability but simultaneously act as a form of regularization by shrinking the Gaussian process prior in ((ref)). This regularization is essential for controlling the nonlinear ill-posedness in the model. The lower bound on $(K_n)_{n=1}^{\infty}$ can be further relaxed in settings where the conditional mean function $ m(W,h) = \mathbb{E}\!\left[\rho(Y,h(X)) \,\middle|\, W\right] $ is known to smooth features of $h$ in a neighborhood of $h_0$.
Theorem (ref) can be extended in several directions. One possibility is to consider a continuously updated version of the quasi-Bayes posterior. In this case, the data-dependent weighting matrix $\widehat{\Sigma}$ may depend pointwise on both $W$ and the prior realization $h$, i.e.\ $\widehat{\Sigma} = \widehat{\Sigma}(W,h)$. The continuously updated quasi-Bayes posterior is then given by
For example, a natural choice is a feasible estimate of the optimal continuously updated weighting matrix: \[ \Sigma(W,h) = \big\{ \, \mathbb{E}[ \rho(Y,h(X)) \rho(Y,h(X))' \mid W ] \, \big\}^{-1}. \] Another possible extension is to generalize the contraction result in Theorem (ref) to settings where the unknown function $h_0$ is not uniquely identified from the data. In this case, the identified set is given by $ \Theta_0 = \big\{ h : \| m(W,h) \|_{L^2(\mathbb{P})} = 0 \big\}.$ Intuitively, regardless of point identification, samples from the quasi-Bayes posterior should concentrate in regions where the quasi-Bayes objective function is minimized, i.e.\ around the identified set $\Theta_0$. Below, we state a version of Theorem (ref) that accommodates both of the preceding extensions. To this end, we impose the following condition on the weighting matrix.
Theorem (ref) establishes the consistency of the continuously updated quasi-Bayes posterior, provided that at least one element of the identified set possesses sufficient regularity relative to the Gaussian process sample paths.
For the remainder of Section (ref), we focus on the case with a standard weighting matrix and a uniquely identified structural function. Extensions to continuously updated weighting and partial identification can be addressed analogously to Theorem (ref).
In this section, we establish contraction rates for the quasi-Bayes posterior. Although Theorem (ref) established consistency, it did not quantify the rate of convergence. In the following analysis, we provide explicit posterior contraction rates.
In our setting, as we illustrate below, the posterior contraction rate is determined by the interplay among $(i)$ the sample path properties of the Gaussian process prior, $(ii)$ the local curvature of the objective function that defines the quasi-Bayes posterior, $(iii)$ the smoothing properties of the $h \mapsto m(W,h)$ locally around $h_0$, and $(iv)$ the basis functions $b^K(W) = (b_1(W), \dotsc, b_K(W))'$ used to construct a first-stage estimate of $m(W,h)$.
The behavior of the nonlinear map $h \mapsto m(W,h)$ can be locally approximated around $h_0$ by a suitable linearization. Depending on the model and the assumptions on the data $\mathcal{D} = (Y,X,W)$, there may be multiple candidates for such a linearization. If the map $h \mapsto m(W,h)$ is sufficiently regular in a neighborhood of $h_0$, the natural choice is the Fréchet derivative at $h_0$, i.e. the unique continuous linear operator $D_{h_0}: L^2(X)\to L^2(W)$ such that \[ \| m(W,h_0+h) - m(W,h_0) - D_{h_0}[h] \|_{L^2(\mathbb{P})} = o(\|h\|_{L^2(\mathbb{P})}) \quad \text{as } \|h\|_{L^2(\mathbb{P})}\to 0. \] Intuitively, if $D_{h_0}[h]$ provides a good local approximation to $m(W,h)$ around $h_0$, then the smoothing properties of the nonlinear map $h \mapsto m(W,h)$ can be studied through the simpler linear operator $h \mapsto D_{h_0}[h]$. In what follows, we relate the smoothing behavior of $D_{h_0}$ to changes in regularity with respect to the orthonormal basis $(e_i)_{i=1}^{\infty}$ defining the Gaussian process in ((ref)). Since the smoothness of $h_0$ is also defined relative to this basis through membership in the Sobolev ball ((ref)), this allows us to analyze the action of $D_{h_0}(\cdot)$ on $(G_{\alpha}, h_0)$ under a common regularity scale. To this end, it will be convenient to define a family of weak norms on $L^2(\mathcal{X})$, obtained by shrinking the Fourier coefficients of a function relative to the basis $(e_i)_{i=1}^{\infty}$. We introduce the following definition:
For $\gamma > 0$ and $\epsilon > 0$, we denote a bounded smooth local neighborhood of $h_0$ by
The following two conditions quantify the smoothing properties of the map $h \rightarrow m(W,h)$ in a local neighborhood of $h_0$ by relating it to a suitable weak norm. \begin{condition2}[Smoothing Link] There exists $\epsilon > 0$ sufficiently small, $\gamma > 0$ and a sequence $\sigma_i \to 0$ such that, for any $M > 0$, there are constants $C_1(M), C_2(M) < \infty$ satisfying $ \| D_{h_0}[h - h_0] \|_{L^2(\mathbb{P})} \leq C_1(M) \| h - h_0 \|_{w,\sigma} $ and $ \| h - h_0 \|_{w,\sigma} \leq C_2(M) \| D_{h_0}[h - h_0] \|_{L^2(\mathbb{P})} $ for every $h \in \Omega(M, \epsilon, \gamma)$. \end{condition2} \begin{condition2}[Local Curvature] There exists $\epsilon > 0$ sufficiently small and $ \gamma > 0$ such that, for any $M > 0$, there exists a constant $B = B(M) < \infty$ satisfying $ \| m(W,h) \|_{L^2(\mathbb{P})} \leq B \| D_{h_0}[h - h_0] \|_{L^2(\mathbb{P})} $ and $ \| D_{h_0}[h - h_0] \|_{L^2(\mathbb{P})} \leq B \| m(W,h) \|_{L^2(\mathbb{P})} $ for every $h \in \Omega(M, \epsilon, \gamma)$.
\end{condition2}
\begin{condition2}[First Stage] Let $\alpha > \gamma$ denote the regularity of the Gaussian process $G_{\alpha}$. There exist sufficiently small $\epsilon, \delta > 0$, a non-increasing function $\varphi: \mathbb{R}_{+} \to \mathbb{R}_{+}$ and a constant $D > 0$ such that, for any $M > 0$,
for all sufficiently large $K$ and $\zeta \in (\alpha - \delta, \alpha)$. \end{condition2} Conditions (ref)--(ref), albeit in varied formulations, are standard in the literature.\footnote{Our conditions are equivalent to the assumptions in chen2012estimation; see, for example, Corollary 5.3 therein. For further discussion on alternative formulations, see also Remark (ref) below.} These conditions can be further weakened to hold with a sequence $\epsilon = \epsilon_n \to 0$ sufficiently slowly. Condition (ref) holds trivially when $h \mapsto m(W,h)$ is linear, as in the NPIV model. If $D_{h_0}^*$ denotes the adjoint, a sufficient (but not necessary) assumption for Condition (ref) is that the self-adjoint operator $D_{h_0}^* D_{h_0}$ diagonalizes in the eigenbasis $(e_i)_{i=1}^\infty$ of the Gaussian process in ((ref)). Stronger versions of Condition (ref) are often imposed in the literature on linear inverse problems with a known operator (e.g. \citealp*{knapik2011bayesian}; \citealp*{gugushvili2020bayesian}).
The intuition behind Condition (ref), following chen2012estimation, is that locally around $h_0$, the map $(h,h_0) \mapsto m(W,h) - m(W,h_0)$ exhibits smoothing properties that are comparable to those of its local linear approximation $(h,h_0) \mapsto D_{h_0}[h - h_0]$. Thus, it is expected that the decay rate of $\varphi(K)$ is of the same order as the sequence $\sigma_K$ in Condition (ref), while $K^{-\zeta/d}$ represents the usual sieve approximation error for bounded smoothness classes $\mathcal{H}^{\zeta}(M)$.
Following standard practice in the literature, we distinguish two regimes of estimation difficulty. The model is said to be mildly ill-posed if $\sigma_K$ and $\varphi(K)$ decay at a polynomial rate, and severely ill-posed if they decay at an exponential rate. The following result establishes contraction rates for the generalized Bayes posterior.
In the literature (e.g. \citealp*{chen2012estimation,chernozhukov2023constrained}), the assumption $\varphi(K) \asymp \sigma_K$ is often imposed, as it corresponds, in a certain sense, to an optimal choice of first-stage approximating functions. Theorem (ref) allows for some degree of misspecification in this choice, with the rates simplifying under the conventional hypothesis (see Corollary (ref) below). For clarity and simplicity of notation, we proceed under the conventional hypothesis for the remainder of the paper.
As a point estimator for $h_0$, we consider the posterior mean
Given the posterior contraction rate in Theorem (ref), the posterior mean, as a point estimator, is expected to converge at a comparable rate. Intuitively, this follows if the posterior probability of the set where contraction fails decays sufficiently quickly. The next result formalizes this intuition.
In this section, we study the limiting quasi-posterior distribution for a class of linear functionals. Let $\mathbf{L}(h_0)$ denote a linear functional of interest—for example, the average value of $h_0(\cdot)$ over an interval or its average derivative. Our analysis focuses on two main questions: $(i)$ What is the limiting quasi-Bayes posterior distribution of $\mathbf{L}(h)$? $(ii)$ Under what conditions do quasi-Bayes credible sets for $\mathbf{L}(h_0)$ attain valid frequentist coverage?
To begin our analysis, we view the linear functional as a map $\mathbf{L} : L^2(\mathcal{X}) \rightarrow \mathbb{R}$. Then, by the Riesz representation theorem, there exists a function $\Phi \in L^2(\mathcal{X})$ such that
The advantage of this representation is that properties of $\mathbf{L}(\cdot)$ (e.g. regularity) can be analyzed through its representer function $\Phi(X)$.
In the preceding sections, the choice of the weighting matrix $\widehat{\Sigma}(\cdot)$ in the quasi-Bayes posterior ((ref)) did not affect the limit theory, provided that the eigenvalues of $\widehat{\Sigma}(\cdot)$ remained asymptotically bounded away from $0$ and $\infty$. Intuitively, under this condition, the rates of convergence can be characterized by analyzing a quasi-Bayes posterior based on the identity weighted objective $h \mapsto \mathbb{E}_n \big( \| \widehat{m}(W,h) \|_{\ell^2}^2 \big)$. To characterize finer aspects of the posterior, it will be necessary to account for the limiting behavior of $\widehat{\Sigma}(\cdot)$ in the analysis. We impose the following low level condition on the limiting behavior of the weights. \begin{condition2}[Limiting Weights] There exists a limit symmetric matrix $\Sigma_0(\cdot)$ such that $ \sup_{w\in\mathcal{W}} \bigl\|\widehat{\Sigma}(w)-\Sigma_0(w)\bigr\|_{op} =O_{\mathbb{P}}(\gamma_n), $ where $ (\gamma_n)_{n=1}^{\infty} $ satisfies \(\gamma_nK_n\to0\). Furthermore, the eigenvalues of \(\Sigma_0(W)\) are uniformly bounded away from zero and infinity: \[ \mathbb{P}\Bigl(c\le\lambda_{\min}\bigl(\Sigma_0(W)\bigr) \le\lambda_{\max}\bigl(\Sigma_0(W)\bigr)\le C\Bigr) =1 \] for some universal constants $c,C > 0$. \end{condition2} We are primarily interested in the setting where \(\Sigma_0(\cdot)\) is an efficient weighting matrix for the conditional moment restriction, so that \(\widehat\Sigma(\cdot)\) may be viewed as a preliminary first‑step estimate of the optimal weighting matrix. In finite‑dimensional GMM models, a celebrated result by chernozhukov2003mcmc establishes the frequentist validity of optimally weighted quasi‑Bayes credible sets. In this section, we provide a nonparametric extension to their results by studying the frequentist coverage of quasi-Bayes credible sets for the functional \(\mathbf L(h_0)\).
As in Section (ref), let \( D_{h_0}(\cdot) \) denote the Fréchet derivative of the map \( h \mapsto m(W,h) \) at \( h_0 \). We denote its adjoint by \( D_{h_0}^* \).\footnote{In defining $D_{h_0}^*$, we view $D_{h_0}$ as a map $(L^2(X) , \| . \|_{L^2(\mathbb{P})}) \mapsto (L^2(W , \| . \|_{L_{\Sigma_0}^2(\mathbb{P})}) $, where $\| . \|_{L^2_{\Sigma_0}(\mathbb{P})}$ denotes the optimal weighted norm $ \| D_{h_0}(h) \|_{L^2_{\Sigma_0}(\mathbb{P})}^2 = \mathbb{E} \left[ D_{h_0}(h)' \Sigma_0(W) D_{h_0}(h) \right]$.} Let $\mathbb{H}$ denote the reproducing kernel Hilbert space (RKHS) of the Gaussian process $G_{\alpha}$. The following condition specifies our main regularity requirements on the representer function $\Phi(\cdot)$.
\begin{condition2}[Regular Functional] There exists $\tilde{\Phi} \in \mathbb{H}$ such that $\Phi = D_{h_0}^* D_{h_0} \tilde{\Phi}$. The first-stage approximation biases of $D_{h_0}[\tilde{\Phi}]$ and $\Sigma_0(W) D_{h_0}[\tilde{\Phi}]$ satisfy:
\end{condition2} The requirement that $\Phi$ lie in a suitable range of the adjoint is a well-known necessary condition for $\sqrt{n}$ estimation of linear functionals, appearing in a variety of settings. For exogenous nonlinear regression models, see \citet*{monard2021statistical}; for NPIV models, see \citet*{severini2012efficiency}, bennett2022inference, \citet*{deaner12trade}; and for NPQIV models, see \citet*{chen2019penalized}. This condition implicitly imposes regularity constraints on $\Phi$. Although extending to more general settings, such as irregular functionals, would be desirable, we view our analysis as an important first step toward a comprehensive nonparametric quasi-Bayes inferential theory.
Given the posterior contraction rate established in Theorem (ref), it suffices, for deriving the distributional limit theory, to restrict our analysis to a quasi-Bayes posterior whose support is contained within local neighborhoods of $h_0$. Specifically, if $\Theta_n$ denotes a sequence of shrinking local neighborhoods around $h_0$, it suffices to focus on the localized posterior:
Let $\delta_n$ denote the posterior contraction rate established in Theorem (ref). In our analysis, we will also make use of the contraction rate $\xi_n$, obtained with the weaker metric $d_w(h,h_0) = \| m(W,h) - m(W,h_0) \|_{L^2(\mathbb{P})}$. As a byproduct of our earlier analysis, it is straightforward to verify that this contraction rate is given by
If $\gamma > 0$ is as in Condition (ref)-(ref), we consider the localized distribution $\mu^{\star}(\cdot \:|\:\mathcal{D}_n)$ obtained through the sequence of smooth local neighborhoods:
where $D,M > 0$ are sufficiently large universal constants.
To connect with the usual linear distributional theory, we quantify the discrepancy between $m(W,h)$ and its linear approximation $D_{h_0}[h-h_0]$ locally around $h_0$. To that end, given any function $h : \mathcal{X} \rightarrow \mathbb{R}$, we denote the remainder obtained from linearizing the map $h \rightarrow m(W,h)$ locally around $h_0$ by
For linear problems such as NPIV (Example (ref)), we have $R_{h_0}(h,W) = 0$ for every $h$. As such, including ((ref)) in the analysis is only relevant for nonlinear models. Analogous to the finite dimensional Euclidean case, the remainder vanishes as $\| h- h_0 \|_{L^2(\mathbb{P})} \rightarrow 0$. The precise rate at which this occurs depends on (among other factors) $(i)$ the ill-posedness in the model, $(ii)$ the regularity of $h$ and $(iii)$ the convergence rate of $ \| h- h_0 \|_{L^2(\mathbb{P})}$.
Let $ \mathcal{M}_n = \{ m(\cdot,h) : h \in \Theta_n \} $ denote the image of $\Theta_n$ under the first stage map $h \mapsto m(W,h)$. As is standard, we quantify the complexity of $\mathcal{M}_n$ through its entropy integral:
where $N(\mathcal{S},d,\delta)$ denotes the usual $\delta-$covering number of a set $\mathcal{S}$ with respect to the metric $d$. The following condition specifies our requirements on the localized support $\Theta_n$, its image $\mathcal{M}_n$ and nonlinear remainder $\{ R_{h_0}(h,W) : h \in \Theta_n \}$.
\begin{condition2} Let $\kappa$ and $t$ denote the local $L^2$ continuity parameters of the generalized residual $\rho(\cdot)$, as defined in Condition (ref). Suppose that:
\end{condition2} Conditions (ref)$(i)$--$(ii)$ arise primarily from empirical process techniques used to control the uniform empirical deviation: \[ \chi_n = \sup_{h \in \Theta_n} \left| \mathbb{E}_n \left[ \widehat{m}(W,h)' \Sigma(W) \widehat{m}(W) \right] - \mathbb{E} \left[ \Pi_{K} m(W,h)' \Sigma(W) \Pi_{K} m(W,h) \right] \right|. \] If we substitute the posterior contraction rate $\delta_n$ and the optimal first-stage sieve dimension sequence $K_n$ from Theorem (ref), Condition (ref) can be reduced to minimum smoothness requirements on the structural function $h_0$ and prior. The dependence on $\kappa$ and $t$ arises because the generalized residual function $\rho(\cdot)$ may be nonlinear and pointwise discontinuous in $h$. Accordingly, our analysis relies on the weaker $L^2(\mathbb{P})$ continuity condition specified in Condition (ref).
The following result establishes that the quasi-Bayes posterior distribution of a regular functional $\mathbf{L}(.) = \langle \cdot , \Phi \rangle_{L^2(\mathbb{P})}$ is well approximated by a suitable Gaussian measure.
The two variances in Theorem (ref) coincide if and only if the quasi-Bayes posterior is optimally weighted. That is, when the weighting matrix is
An important implication of Theorem (ref) is that optimally weighted quasi-Bayes credible sets, centered around the posterior mean, attain asymptotically exact frequentist coverage. Specifically, given a linear functional \( \mathbf{L}(\cdot) \) and a significance level \( \gamma \in (0,1) \), define \[ c_{1-\gamma} = (1-\gamma) \text{ quantile of } \left| \mathbf{L}(h) - \mathbf{L}\left( \mathbb{E}[h \mid \mathcal{D}_n] \right) \right|, \quad h \sim \mu(\cdot \mid \mathcal{D}_n). \] The quasi-Bayes credible set at level \( \gamma \) is defined as: \[ C_n(\gamma) = \left\{ t \in \mathbb{R} : \left| t - \mathbf{L}\left( \mathbb{E}[h \mid \mathcal{D}_n] \right) \right| \leq c_{1-\gamma} \right\}. \]
To the best of our knowledge, Theorem (ref) and Corollary (ref) provide the first nonparametric quasi-Bayes inferential guarantees in the literature. These results extend classical quasi-Bayes inferential results for parametric GMM chernozhukov2003mcmc to nonparametric conditional moment restriction models.
In this section, we present additional simulation evidence on the finite-sample performance of quasi-Bayes posteriors. Whereas Section (ref) focused on structural functions with a univariate regressor, here we consider settings with multivariate regressors.
Specifically, we examine multivariate generalizations of the designs in newey2003instrumental, santos2012inference, \citet*{chernozhukov2015constrained}, chetverikov2017nonparametric, and \citet*{chen2025adaptive}, which we denote as NP, S, CNS, CW, and CCK, respectively. These generalizations are constructed to mimic the endogeneity structure and ill-posedness of the original univariate designs.\footnote{GPT-5 assisted in the construction of these generalizations.} The structural functions are: \( \textstyle
\) In these designs, the endogenous regressor is five-dimensional, \(X \in \mathbb{R}^5\), and the instrument is two-dimensional, \(W \in \mathbb{R}^2\). The structural functions extend those used in the original univariate designs, and collectively span a reasonable spectrum of functional complexity. Beyond maintaining a similar endogeneity structure, we also scaled up the variance of the disturbances to ensure that the signal-to-noise ratios remain comparable to, or smaller than, those in the original univariate designs. All details are provided in Appendix (ref).
In endogenous models with multivariate regressors, it is very challenging to estimate the structural function using classical methods. Indeed, with a five dimensional endogenous regressor, even a minimal tensor-product sieve with three terms per coordinate yields $J = 3^5 = 243$ basis functions. In all designs, 2SLS estimation based on this tensor product produced an extremely large and unstable risk. This mirrors the univariate behavior in Table (ref), except that in higher dimensions the minimal feasible $J$ is already prohibitively large.
Let QB denote the quasi-Bayes posterior mean, based on a first-stage thin-plate spline with dimension $K=15$ and a Whittle–Matérn Gaussian process prior. The same prior and implementation algorithm are used across all designs and both sets of restrictions (see Appendix (ref) for details). For comparison, we also report nonparametric regression estimates using random forests (RF), implemented via the ranger package in R.
Random forests (RF) are a reliable supervised learning method for high-dimensional regression and are expected to capture much of the variation in the structural functions. However, because of the non-trivial endogeneity in the designs, it exhibits substantial bias. The designs in Table (ref) span a wide range of structural function complexities and endogeneity patterns, with some expected to serve as relatively challenging stress tests. In practice, we expect our methods to perform considerably better in more conventional settings.
The results in Table (ref) demonstrate that the quasi-Bayes estimators perform well and are viable in higher dimensions. In particular, the estimators are accurate and stable across both restrictions. This is especially noteworthy since nonparametric quantile IV (NPQIV) estimation is often regarded as a substantially more difficult problem due to its nonlinear and discontinuous generalized residual. Together with the simulation evidence in Section (ref), our findings suggest that quasi-Bayes estimators may provide a broadly useful toolkit for the large class of nonlinear restrictions frequently encountered in applied work.
Figure (ref) plots a sample realization of quasi-Bayes predicted vs true values on a generated test data. The predictions closely follow the trajectory of the true values, concentrating around the 45-degree line of equality. Figure (ref) plots the associated fit for the biased OLS predictions.
As a final remark, it would be desirable to compare the quasi-Bayes estimators with other nonparametric alternatives. However, we are not aware of any reliable implementations for general conditional moment models with multivariate regressors. To the best of our knowledge, our simulation study also provides the first nonparametric risk estimates for quantile IV models with multivariate regressors.
In this section, we apply our methodology to estimate firm-level production functions in Chile, using data from the national census of manufacturing plants conducted by Chile’s Instituto Nacional de Estadística. This dataset is frequently employed in studies of firm-level production functions (e.g. \citealp*{levinsohn2003estimating,gandhi2020identification}). Our analysis focuses on the food products industry, one of the country’s largest manufacturing sectors. We use firms with more than 10 employees and complete observations for the years 1979–1996.
Let $y_{it}, k_{it}, l_{it}$ denote the logarithms of gross output, capital, and labor, respectively, and let $m_{it}$ denote intermediate inputs (fuels, materials, and electricity). All variables are in real terms. Consider the structural value-added production model \[ y_{it} = F(l_{it}, k_{it}) + \omega_{it} + \varepsilon_{it}, \] where $F(\cdot)$ is the production function in inputs $(l,k)$, $\varepsilon_{it}$ are exogenous shocks unobserved by the firm, and $\omega_{it}$ are first-order Markov shocks observed (or predictable) by the firm prior to its input decisions at time $t$. We assume $\omega_{it}$ is a deterministic function of inputs, $\omega_{it} = \tilde f_t(k_{it}, l_{it}, m_{it})$, for some function $\tilde f_t$. One interpretation of this specification, following ackerberg2015identification, is that the gross-output production function is Leontief in the intermediate input. Define the conditional means \[ g(\omega_{it-1}) = \mathbb{E}[\omega_{it} \mid \omega_{it-1}] \quad, \quad \Phi_t(l_{it}, k_{it}, m_{it}) = \mathbb{E}[y_{it} \mid l_{it}, k_{it}, m_{it}]. \] Note that, since $\varepsilon_{it}$ is exogenous noise, the function $g(\cdot)$ can be interpreted as the conditional mean regression of $\Phi_{t}(l_{it}, k_{it}, m_{it}) - F(l_{it}, k_{it})$ on $\Phi_{t-1}(l_{it-1}, k_{it-1}, m_{it-1}) - F(l_{it-1}, k_{it-1})$. If $\mathcal{I}_{t}$ denotes the firm’s information set at time $t$, it is shown in ackerberg2015identification that $ F(\cdot)$ satisfies the conditional moment restriction:
In most industries, it is assumed that firms choose labor $l_{it}$ after period $t-1$. Under this timing assumption, the natural information set, as in ackerberg2015identification, is $ \mathcal{I}_{t-1} = \{ k_{it}, \, l_{it-1}, \, \Phi_{t-1} \} $. We use the same information set in our analysis.
The functions $g(\cdot)$ and $\Phi_{t-1}(\cdot)$ are smooth, low-dimensional regressions and can therefore be estimated accurately with standard nonparametric methods. In practice, $\Phi_{t-1}(\cdot)$ is typically estimated using a flexible sieve regression (e.g. splines). Similarly, for any input function $\tilde{F}$, the output of $g(\cdot)$ in the restriction is obtained from a one-dimensional conditional mean regression, typically implemented with a flexible polynomial. We adopt this approach and thus treat both functions as known for the restriction in ((ref)). Further implementation details are provided in Appendix (ref).
We aim to estimate the production function that satisfies the conditional moment restriction in ((ref)). This is a particularly challenging problem, as the restriction defines a complex and highly nonlinear inverse problem in $F(\cdot)$.
Figure (ref) shows the posterior mean estimator $\widehat{F}(k,l) = \mathbb{E}[F(k,l) \mid \mathcal{D}_n]$ as a function of log capital $k$, with labor fixed at selected quantiles. For each labor quantile, the production function displays the familiar S-shape: convex at low $k$, where additional capital raises productivity at an increasing rate, and concave at higher $k$, where diminishing returns set in. Consequently, the marginal product in Figure (ref) first increases with capital but eventually declines, yielding the classical inverted-U pattern.
Figure (ref) shows the estimated production function $\widehat{F}(k,l)$ as a function of log labor $l$, holding capital fixed at selected quantiles. At low to moderate capital quantiles, the function is roughly linear for small values of $l$, becomes convex at intermediate levels, and turns concave at higher levels. By contrast, at very high capital quantiles, the function begins at a higher level of output and maintains an almost linear trajectory with a steep slope over most of the range of $l$, turning concave only at higher values. Figure (ref) illustrates these patterns via the corresponding marginal product curves.
In the data, real capital at the 0.5, 0.75, and 0.95 quantiles equals 740.96, 3656.74, and 24,325.68, respectively, indicating a sharp increase at the upper end of the distribution. One interpretation of these patterns is that they reflect how labor interacts with available capital. With low to moderate capital, complementarities cause output to expand more rapidly as labor increases before diminishing returns set in, yielding convexity followed by concavity. With abundant capital, each worker is already highly productive, so output rises almost linearly with a steep slope in labor until very high levels, where diminishing returns set in.
As a final remark, we note that the identifying restriction for $F(\cdot)$ in ((ref)) is complex and highly non-linear. It is therefore noteworthy that our procedures are still able to recover reasonable and meaningful features of $F(\cdot)$ from this restriction alone. To our knowledge, this represents the first fully nonparametric estimate of $F(\cdot)$, obtained without imposing any predetermined parametric structure. Beyond serving as a valuable nonparametric benchmark, these estimates may also provide guidance for the empirical design of approximating parametric specifications. In particular, our findings suggest a preference for specifications that can capture flexible variation in marginal products across input levels.
This paper develops a generalized Bayes framework for a broad class of nonparametric conditional moment restriction models. Simulations demonstrate that the proposed procedures are viable and perform well. We expect these methods to be broadly useful, particularly in ill-posed settings or when closed-form solutions are unavailable. As an empirical illustration, we apply the methodology to estimate nonparametric production functions using Chilean plant-level data. We conclude with a few remarks and outline possible extensions.
In Section (ref), we motivated quasi-Bayes procedures as an attractive form of data-driven regularization for endogenous nonparametric inverse problems. An additional advantage is in their flexibility to incorporate application specific information. For instance, extending Remark (ref), one may specify informative priors centered at a fixed structural function $\widetilde{h}(\cdot)$. In many applications (e.g., \citealp*{adao2017nonparametric, bergquist2020competition}), researchers may have strong microfounded preferences for a parametrically estimated $\widetilde{h}(\cdot)$, yet still wish to accommodate potential misspecification.
As with all nonparametric methods, some degree of finite-sample tuning can often improve performance. In our setting, following Remark (ref), partial tuning of the Gaussian process covariance hyperparameter $\theta = (\sigma, \ell)$ can be beneficial. When the regressors are normalized, a reasonable default is to set $\ell = 1$ and choose $\sigma$ near the scale of the observables. In nonparametric regression with Gaussian errors, it is standard practice (e.g. williams2006gaussian) to empirically select $\theta$ by maximizing the Bayesian marginal likelihood. Writing the prior dependence on $\theta$ as $d\mu(h \mid \theta)$, the natural analogue in our framework is to choose $\theta$ by maximizing the quasi-Bayes marginal likelihood:
In practice, evaluating this normalizing factor over a large grid can be computationally challenging. An intermediate strategy is to place a weakly informative prior on $\theta$, run a short exploration phase in which we sample from the full posterior over $(h, \theta)$, and then fix $\theta$ at $\hat{\theta}$—the posterior mean computed from the latter part of this exploration phase. Then, proceed with full posterior sampling from the quasi-Bayes posterior $d\mu(h \mid \mathcal{D}_n, \hat{\theta})$. This is the approach we adopt in our implementation. In high-dimensional settings, a common approach for updating $\theta$ during the exploration phase is via slice sampling steps murray2010slice.
The first-stage regression in our procedures can use any available source of variation, including both continuous and discrete instruments. Furthermore, there is no requirement that the number of functions in the first stage exceed a fixed threshold. This is in contrast to classical IV 2SLS, which requires at least $K \geq J$ functions in the first stage to estimate a $J$-dimensional second-stage parameter. This flexibility should be particularly valuable in empirical settings where researchers have mixed sources of variation and substantially fewer instruments than endogenous regressors.
We use the same implementation algorithm across all settings considered in this paper, discussed further in Appendix (ref). Briefly, the approach consists of preconditioned Crank–Nicolson (pCN) steps applied to a suitable non-centered parametrization of the Gaussian process sample paths.\footnote{pCN proposals are frequently employed to target infinite-dimensional posteriors that arise in inverse problems with Gaussian process priors cotter2013mcmc,nickl2023bayesian.} We view this as an attractive feature, as it suggests that the same algorithm, perhaps with only minor modifications, can be applied broadly.
For ease of exposition, we focused on a single structural function $h_0(\cdot)$ that depends on the entire endogenous vector $X$. Adapting the framework to settings with multiple structural functions and restrictions defined on different subcomponents of the observables is straightforward, though notationally more cumbersome.
Our limit theory is developed for a class of infinite-dimensional Gaussian process (GP) priors. Extending the results to other widely used prior classes (e.g., chipman2012bart) or to priors that directly impose specific shape restrictions would be valuable. For GP priors in particular, there is already a substantial literature on enforcing such constraints in regression models (e.g. lin2014bayesian).
Section (ref) develops, to our knowledge, the first inferential results for a nonparametric quasi-Bayes framework, extending classical parametric GMM results chernozhukov2003mcmc. The analysis focused on regular, $\sqrt{n}$-estimable functionals. A natural direction for future work is to broaden the framework to irregular functionals that are slower than $\sqrt{n}$-estimable, similar to the frequentist analysis in chen2015sieve.