On Gaussian Process Priors in Conditional Moment Restriction Models

On Gaussian Process Priors in Nonparametric Conditional Moment Restriction Models

Introduction

This paper considers a general conditional moment restriction model where an unknown function $h_0$ is identified from the restriction

Here, $Y \in \mathbb{R}^{d_y}$ is a vector of observables, $X \in \mathbb{R}^{d}$ is a vector of regressors and $W \in \mathbb{R}^{d_w}$ is a vector of conditioning (or instrumental) variables. The vector $\rho(.) = [ \rho_{1}(.) , \dots , \rho_{d_{\rho}}(.) ]' $ is a known $d_{\rho}$ dimensional vector of generalized residual functions and $h_0(.)$ is an unknown function of interest. The model in ((ref)) has a long history in econometrics and statistics (e.g. liao2011posterior; chen2012estimation). As a special case, it nests the nonparametric instrumental variable model (NPIV) studied in newey2003instrumental; hall2005nonparametric; \citet*{blundell2007semi} and the nonparametric quantile instrumental variable (NPQIV) model studied in \citet*{chernozhukov2007instrumental,horowitz2007nonparametric}.

The conditional moment restriction model in ((ref)) represents a general class of ill-posed inverse problems with an unknown and possibly nonlinear operator. In this setup, the (possibly nonlinear) operator $h \rightarrow m(W,h) = \mathbb{E}[\rho(Y,h(X))|W]$ smoothes out features of $h$. Moreover, the operator is not known as it depends on the true data generating process through the population conditional expectation operator $\mathbb{E}(.|W)$. One consequence of this, from a statistical perspective, is that feasible procedures which replace $m(W,h)$ with a finite sample analog $\widehat{m}(W,h)$ are sensitive to features of the generalized residual $\rho(.)$, that would otherwise be smoothed out in the case of a fully known operator. In particular, if $\rho(.)$ is pointwise nonsmooth or discontinuous in its arguments, as in the case of a NPQIV model, a feasible estimator $\widehat{m}(W,h)$ may inherit similar properties in finite samples.

In this paper, we study the quasi-Bayesian posterior distribution that arises from the conditional moment restriction in ((ref)). If we denote the conditional mean of the residual at a function $h$ by $m(W,h) = \mathbb{E}[\rho(Y,h(X))|W]$, the unknown function is identified by the restriction $ \mathbb{E} \big( \| m(W,h_0) \|_{\ell^2}^2 \big) = 0 $. As such, estimation can be based on the finite sample objective function $h \rightarrow \mathbb{E}_n \big[ \widehat{m}(W,h) ' \widehat{\Sigma}(W) \widehat{m}(W,h) \big] $, where $\widehat{m}(W,h)$ denotes a feasible estimator of $m(W,h)$ and $\widehat{\Sigma}(W)$ is a positive semi-definite weighting matrix. In a quasi-Bayes framework, we view this objective function as a pseudo-likelihood. When combined with a prior $\mu$, this leads to the quasi-posterior distribution

The quasi-posterior in ((ref)) was proposed by liao2011posterior. Using sieve based priors, they established posterior consistency. Posterior contraction rates, also using sieve based priors, were obtained in kato2013quasi for the special case of a nonparametric instrumental variable model. The main results of this paper develop the limit theory for the quasi-Bayes posterior when $\mu$ belongs to a class of sufficiently smooth Gaussian process priors. As a first step, our analysis extends the preceding work on contraction rates for the quasi-posterior in that $(i)$ we establish contraction rates for the general conditional moment restriction model in ((ref)) and $(ii)$ we do so using a class of Gaussian process priors. Beyond contraction rates, we also provide conditions under which a Bernstein–von Mises (BvM) theorem holds for the quasi-posterior distribution. To be specific, if $h \rightarrow \mathbf{L}(h)$ is a sufficiently smooth linear functional, we show that the induced quasi-posterior distribution of $\mathbf{L}(h)$ can be well approximated (in the sense of weak convergence in probability) by a suitable Gaussian measure. We use this to provide frequentist guarantees for quasi-Bayesian credible sets that are centered around the posterior mean. In particular, we show that such credible sets have asymptotically exact frequentist coverage, provided that the quasi-Bayes posterior in ((ref)) is optimally weighted. This extends classical results (e.g. chernozhukov2003mcmc) on the frequentist validity of optimally weighted quasi-Bayes credible sets for parametric generalized method of moments (GMM) models.

Our approach to obtaining contraction rates is similar to the general strategy employed in inverse problems (e.g. knapik2018general) in that we first obtain contraction rates for a suitable direct problem, i.e $ d_{w}(h,h_0) = \| \Pi_n \big[m(W,h) - m(W,h_0) \big] \|_{L^2(\mathbb{P})} $ where $\Pi_n : L^2(\mathbb{P}) \rightarrow L^2(\mathbb{P})$ is a sample size dependent orthogonal projection operator. Similar to the approach in giordano2020consistency; \citet*{monard2021consistent,monard2021statistical}, we consider rescaled Gaussian process priors which concentrate (with high probability) on bounded subsets of a sufficiently smooth function class. By combining this with contraction rates for the direct problem, we obtain contraction rates for stronger metrics such as $d(h,h_0) = \| h - h_0 \|_{L^2(\mathbb{P})}$. To obtain asymptotic Bernstein-von Mises Gaussian approximations, we make use of posterior local fluctuation and change of parametrization arguments that are frequently used in the analysis of Gaussian process priors in density or regression frameworks (e.g. castillo2015bernstein; monard2021statistical). In the setup considered here, these arguments are suitably modified to account for $(i)$ a quasi-posterior as opposed to a traditional likelihood; $(ii)$ a general conditional moment restriction as in ((ref)); $(iii)$ an unknown and possibly nonlinear operator $h \rightarrow m(W,h)$ that must be estimated as a first step and $(iv)$ a generalized residual function $\rho(.)$ that may be nonlinear and pointwise nonsmooth in its arguments.

This setting considered in this paper is closely related to several strands of literature that study statistical inverse problems in distinct setups. Rates of convergence and confidence sets for the model in ((ref)) was established in chen2012estimation, chen2015sieve under a frequentist sieve based framework. There is a large literature (e.g. \citealp*{knapik2011bayesian}; \citealp*{agapiou2013posterior}; florens2016regularizing; \citealp*{gugushvili2020bayesian}) that studies Gaussian priors within the context of linear inverse regression models with Gaussian noise and a known linear operator. In these cases, the model is conjugate with a known Gaussian posterior distribution and statistical properties can be analyzed directly. In the nonlinear case, monard2021statistical consider Bayesian inference with Gaussian process priors on a class of nonlinear inverse regression models with Gaussian noise and a known nonlinear operator. Similar to their analysis, we consider additional posterior regularization obtained through scaling the Gaussian process. In our conditional moment setting, this scaling is crucial to $(i)$ control the nonlinearity and ill-posedness of the inverse problem and $(ii)$ obtain high probability guarantees on the behavior of the first stage estimator $\widehat{m}(W,h)$ used to approximate the unknown nonlinear operator $h \rightarrow m(W,h)$.

The paper is organized as follows. Section (ref) provides a brief review of Gaussian process priors. Section (ref) introduces the quasi-Bayes framework and main assumptions. Section (ref) develops the quasi-Bayes limit theory and main theoretical results. Section (ref) contains proofs and auxiliary results for all the statements in the main text.

Notation

Given positive sequences $(x_n)_{n=1}^{\infty} $and $ (y_n)_{n=1}^{\infty}$, we write $x_n \lessapprox y_n$ if $\limsup_{n \rightarrow \infty} x_n/y_n < \infty $ and $x_n \asymp y_n$ if $x_n \lessapprox y_n \lessapprox x_n$. Given a positive definite matrix $\Sigma \in \mathbb{R}^{k \times k }$, the induced inner product and norm on $\mathbb{R}^k$ is denoted by $ \langle ., . \rangle_{\Sigma} $ and $ \| . \|_{\Sigma}$, respectively. That is, $ \langle u , v \rangle_{\Sigma} = u' \Sigma v$. The Euclidean norm $(\Sigma = I_k)$ is denoted by $\| .\|_{\ell^2}$. Let $\| . \|_{\infty}$ denote the usual supremum norm on functions and vectors. We use $\mathbb{E}$ and $\mathbb{P}$ to denote the usual expectation and probability operators. Let $\mathbb{E}_n$ and $\mathbb{P}_n$ denote the empirical analog of $\mathbb{E}$ and $\mathbb{P}$, respectively. Given a random vector $Z$, let $L^2(Z)$ denote the usual $L^2$ space of $\mathbb{R}$ valued functions that are measurable with respect to the $\sigma$ algebra generated by $Z$. Similarly, the $L^2$ space corresponding to $\mathbb{R}^d$ valued functions is denoted by $L^2(Z,\mathbb{R}^d)$. In both cases, the $L^2$ norm is denoted by $\| f \|_{L^2(\mathbb{P})}^2 = \int \| f(Z) \|_{\ell^2}^2 d \mathbb{P}$. Given a cube $\mathcal{X} \subset \mathbb{R}^d$, we use $ \mathbf{H}^{p} = ( \mathbf{H}^p(\mathcal{X}), \| . \|_{\mathbf{H}^p}) $ to denote the usual $p$-Sobolev space of functions on $\mathcal{X}$. The $L^2$ space with respect to the Lebesgue measure on $\mathcal{X}$ is denoted by $L^2(\mathcal{X})$. We denote a $p$-Sobolev ball of radius $M > 0 $ by $\mathbf{H}^p(M) = \{ h \in \mathbf{H}^p : \| h \|_{\mathbf{H}^p} \leq M \}$.

Review

In this section, we briefly review Gaussian random elements and the related topic of generating covariance operators and Hilbert scales through self-adjoint operators. For further details on Gaussian process priors, we refer to ghosal2017fundamentals.

A Gaussian process over an index set $T$ is a stochastic process $\{ G_t : t \in T \}$ such that the vector $(G_{t_1} , \dots , G_{t_k})$ is multivariate normal distributed, for every $t_1,\dots,t_k \in T$ and $k \in \mathbb{N}$. Furthermore, if $\{G_t :t \in T \}$ is a Borel measurable random element with sample paths in a separable subset of $\ell^{\infty}(T)$, then $G$ is a Gaussian random element in this space. The covariance of a Gaussian process $\{ G_t : t \in T \}$ may either refer to the operator $\Lambda$ in Definition (ref) (when the sample paths lie in a Hilbert space) or to the covariance function $C(s,t) = \text{Cov}(G_s,G_t) $.

Consider a mean-zero Gaussian process $G$ with realizations in a separable Hilbert space $\mathcal{H}$ with covariance operator $\Lambda$. By the spectral theorem, there exists an orthonormal basis of eigenfunctions $(e_i)_{i=1}^{\infty} \subset \mathcal{H}$ that diagonalizes the operator $\Lambda$. Furthermore, if $\lambda_i $ denotes the non-negative eigenvalue associated to $e_i$, the sequence is of trace class $(\sum_{i=1}^{\infty} \lambda_i < \infty)$ and $G$ admits an expansion of the form

Conversely, given any non-negative sequence $(\lambda_i)_{i=1}^{\infty}$ with $\sum_{i=1}^{\infty} \lambda_i < \infty$, the representation in ((ref)) defines a Gaussian random element on $(\mathcal{H}, \langle . \: , . \: \rangle_{\mathcal{H}})$. As the following definition illustrates, this leads to a family of Gaussian random elements that differ only by a smoothness scale.

The family $\{ \mathcal{H}^{\beta} : \beta \in \mathbb{R} \}$ in Definition (ref) is an example of a Hilbert or smoothness scale. More generally, Hilbert scales can be generated through a densely defined self-adjoint operator $L$. The following definition clarifies this connection.

If the operator $L$ in the preceding definition admits a compact self-adjoint inverse $L^{-1}$, we can define a Gaussian Process by viewing $L^{-s}$ for $s > 0$ as a covariance operator. In this case, the Gaussian Process can be expressed in a similar form to ((ref)), where $(e_i)_{i=1}^{\infty}$ are the eigenfunctions of $L^{-1}$. The following example illustrates the essential idea for a commonly used class of priors.

Framework and Assumptions

In this section, we introduce the conditional moment restriction (CMR) framework and state our main assumptions on the model.

Model

Suppose we observe a random sample $\mathcal{Z}_n = \{ (X_i,Y_i,W_i) \}_{i=1}^n $, where $Y \in \mathbb{R}^{d_y}$ is a vector of observable variables, $X \in \mathbb{R}^{d}$ is a vector of regressors with support $\mathcal{X} \subset \mathbb{R}^d$ and $W \in \mathbb{R}^{d_w}$ is a vector of conditioning (or instrumental) variables with support $\mathcal{W} \subset \mathbb{R}^{d_w}$. We are interested in a structural function $h_0(X) : \mathcal{X} \rightarrow \mathbb{R}$ that is assumed to satisfy the conditional moment restriction:

Here, $\rho(.) = [ \rho_{1}(.) , \dots , \rho_{d_{\rho}}(.) ] $ is a vector of generalized residual functions with functional forms that are assumed known up to the structural function $h_0(.)$.

Many commonly used statistical and econometric models can be reformulated as a conditional moment restriction of the form in $(\ref{cmr})$. In Example (ref), the generalized residual function is a linear function of $h_0$, while in Example (ref), the residual function is nonlinear and nonsmooth in $h_0$. Thus, Examples (ref) and (ref) may be seen as representative of two distinct classes of statistical models, characterized by regularity of the generalized residual function.

Given a function $h: \mathcal{X} \rightarrow \mathbb{R}$, we denote the conditional mean of the residual at $h$ by

As the true function is assumed to satisfy $m(W,h_0) = \mathbf{0}$, the function can be identified as a minimizer of the objective function

where $\Sigma(W) \in \mathbb{R}^{d_{\rho} \times d_{\rho}} $ is a suitable positive definite matrix. As the true data generating process is not assumed to be known, it is infeasible to work with $Q(h)$ directly. Denote by $\widehat{m}(W,h)$ any feasible estimator of $m(W,h)$ in $ (\ref{condm}) $. The usual conditional moment restriction objective function is given by

where $\widehat{\Sigma}(W) = \widehat{\Sigma}_n(W) $ is a (possibly data dependent) positive definite weighting matrix. In a quasi-Bayes framework, ((ref)) is viewed as a psuedo-likelihood for the model. Given a prior probability measure $\mu$, the quasi-Bayes posterior induced from $\mu$ and ((ref)) is denoted by

Assumptions

In this section, we state our main assumptions on the model and data generating process.

Assumption (ref) is standard. We can always ensure $X,W$ are contained in cubes by transforming them through suitable homeomorphisms $ \tilde{X} = \Phi(X) $ and $ \tilde{W} = \Psi(W)$. The conditional moment restriction in ((ref)) then holds with $\tilde{X} , \tilde{W}$ and $ \tilde{h}_0 = h_0 \circ \Phi^{-1} $. The bounded density condition facilitates the analysis as it implies, among other things, that $ \| f \|_{L^2(\mathbb{P})}^2$ and $f \rightarrow \int_{\mathcal{X}} \left| f(x) \right|^2 dx $ are equivalent metrics over $L^2(X)$.

Assumption (ref) is similar to conditions frequently imposed in the literature (e.g. \citealp*{chen2003estimation}) to facilitate analysis involving non smooth objective functions. In particular, it allows for a pointwise discontinuous residual function $\rho(.)$. However, it requires that the residual function be uniformly continuous in $L^2(\mathbb{P})$ expectation. The parameter $\kappa$ determines the modulus of continuity. It holds with $\kappa = 1$ for the NPIV model (Example (ref)) and $\kappa=1/2$ for the NPQIV model (Example (ref)).

Before stating the remaining assumptions, we fix any sufficiently large $t \geq d/\kappa$ that satisfies Assumption (ref).

Assumption $\ref{residuals2}$ imposes weak moment bounds on the residual function. The assumption is trivially satisfied with bounded residual functions, such as in a NPQIV model. For more general cases, observe that if $t > d/2$, the Sobolev inequality evans2022partial implies that $\mathbf{H}^t$ embeds into a H\"{o}lder space. In particular, the functions in $\mathbf{H}^t(M)$ are bounded in $\| . \|_{\infty}$ norm. In most cases, this can be used to verify Assumption $\ref{residuals2}$ directly.

Assumption (ref)$(i)$ is a standard identification condition for the conditional moment restriction model. We consider a relaxation of this assumption in Section (ref). Assumption (ref)$(ii)$ imposes that the conditional mean function $h \rightarrow m(W,h)$ is Lipschitz over any fixed Sobolev ball $\mathbf{H}^t(M)$ of radius $M$. This is made for convenience and can be relaxed further.

Assumption (ref) is a mild differentiability restriction. As we illustrate in Section (ref), this allows us to study the behavior of the map $h \rightarrow m(W,h)$ in a local neighborhood around $h_0$ through properties of its linearization $h \rightarrow D_{h_0}[h]$.

Main Results

Setup

We denote by $b^K(W) = ( b_1(W), \dots , b_K(W) )'$ a vector of first stage basis functions. Here, $K \in \mathbb{N}$ denotes the dimension of the basis. Denote by $\mathcal{V}_K$, the linear subspace of $L^2(W)$ spanned by the basis. Let $\Pi_K : L^2(W) \rightarrow \mathcal{V}_K $ denote the $L^2(\mathbb{P})$ orthogonal projection onto $\mathcal{V}_K$. Given a function $h: \mathcal{X} \rightarrow \mathbb{R}$, we estimate the conditional mean function $m(W,h) = \mathbb{E}[\rho(Y,h(X))|W]$ using the empirical analog of $\Pi_K$. Given the observed data $\mathcal{Z}_n = \{ (X_i,Y_i,W_i) \}_{i=1}^n $, we estimate the conditional expectation using a least squares projection onto $\mathcal{V}_K$. That is,

This choice is not crucial towards obtaining contraction rates but, as a closed form expression, greatly facilitates the analysis in obtaining inferential results.

As a first step towards defining our prior probability measure $\mu$, we consider a family of Gaussian process priors $G_{\alpha} = \{ G_{\alpha}(x) : x \in \mathcal{X} \}$ that are indexed by a regularity hyperparameter $\alpha \in \mathcal{L} \subset \mathbb{R}_+$. As in Section (ref), we assume that the regularity hyperparameter $\alpha$ influences the family through the coefficients on a series expansion with respect to particular basis or through the exponent of an unbounded self-adjoint operator (e.g. the Laplacian) on $\mathcal{X}$. In either case, without loss of generality\footnote{If the mapping $\alpha \rightarrow \lambda_{i,\alpha} $ influences the exponent in a different way, the results can be stated in terms of the induced exponent $ s(\alpha) $, i.e $\lambda_{i,\alpha} \asymp i^{-s(\alpha)}$.}, we can express $G_{\alpha}$ as

where $\lambda_{i,\alpha} \asymp i^{-(1+ 2 \alpha/d)} $, $ Z_i \stackrel{i.i.d}{\sim} N(0,1)$ 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 require that the sample paths of the Gaussian process $G_{\alpha}$ (for $\alpha \in \mathcal{L}$) belong to a separable linear subspace of the Sobolev space $\mathbf{H}^t$ for some $t > d/2$.\footnote{By the Sobolev inequality evans2022partial, $t > d/2$ implies that same paths of $G$ belong, almost surely, to a H\"{o}lder space. In particular, this ensures that the sample paths are continuous. As the H\"{o}lder space is separable with respect to the $\| . \|_{\infty}$ norm, $G$ can be viewed as a Gaussian random element on $(C(\mathcal{X}), \|. \|_{\infty})$.} In fact, due to the possible non smoothness of the residual function $\rho$, we will generally require more stringent conditions (see Assumption (ref)) on the minimal such $t$. In most cases, this can be viewed as a restriction on the hyperparameter set by considering $ \mathcal{L} \subseteq [\underline{\alpha}, \infty)$ for some minimum regularity $\underline{\alpha} > 0$. Given a Gaussian process $ G_{\alpha}$ of regularity $\alpha$ and first stage sieve dimension $K$, we consider the prior distribution

We scale the Gaussian process by the first stage sieve dimension (up to a log term) to provide additional regularization. Let $\widehat{\Sigma}(W) $ denote a (possibly data dependent) positive semidefinite weighting matrix. The quasi-Bayes posterior induced from $(\mu, \widehat{\Sigma} )$ and the conditional moment restriction in ((ref)) is given by

In the analysis that follows, we will frequently measure regularity with respect to the orthonormal basis $(e_i)_{i=1}^{\infty}$ in ((ref)). To that end, we define the p-Sobolev space and p-Sobolev ball relative to the orthonormal basis $(e_i)_{i=1}^{\infty}$ by

Consistency

In this section, we establish consistency of the quasi-Bayes posterior in ((ref)). From the prior representation in ((ref)), it is expected that posterior limit theory depends on some interplay between sample path realiziations of the Gaussian process $G_{\alpha}$ and the quasi-Bayes objective function. Intuitively, given a function class $\mathcal{D}$, if the conditional moment restriction in ((ref)) arises from a nonlinear and/or nonsmooth residual function $\rho(.)$, further (smoothness) restrictions on $\mathcal{D} $ are necessary to quantify the uniform sampling uncertainty of the map $ \mathcal{D} \ni h \rightarrow \mathbb{E}_n\big[ \widehat{m}(W,h)' \widehat{\Sigma}(W) \widehat{m}(W,h) \big]$. In our setting, it suffices for the Gaussian process to be contained (with high probability) in a sufficiently regular function class $\mathcal{D}$. This is formalized in the following condition.

\begin{condition2} The Gaussian process $\{G_{\alpha}(x) : x \in \mathcal{X} \}$ in ((ref)) is a Gaussian random element (in the sense of Definition (ref)) on a separable subspace of the Sobolev space $\mathbf{H}^t$, where $t$ is as in Assumption (ref).

Beyond sampling uncertainty, consistency also depends on the weighting matrix $\widehat{\Sigma}(.)$ and the first stage basis functions $b^K(W) = (b_1(W),\dots,b_K(W))'$ used to construct the estimator $\widehat{m}$. At a minimum, we will impose the following structure. \end{condition2}

\begin{condition2} $(i)$ The matrix $G_{b,K} = \mathbb{E} \big( [b^K(W)] [ b^K(W) ]' \big)$ is positive definite for every $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 away from $0$ and $\infty$: $\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 $, we have $ \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} If Assumption (ref)$(i)$ holds, Condition (ref)$(i)$ is satisfied by splines, Cohen–Daubechies–Vial (CDV) wavelets and Fourier series (see e.g. chen2015optimal; belloni2015some).

Theorem (ref) shows 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 weak, it ensures that the least squares estimator $\widehat{m}(w,h)$ is well defined and that it uniformly approximates its population analog $\Pi_K m(w,h)$. On the other end, the theorem does enforce a strict lower bound on how slowly the first stage basis can grow.\footnote{The restriction on $(K_n)_{n=1}^{\infty}$ can be weakened even further in settings where the conditional mean function $m(W,h) = \mathbb{E}[\rho(Y,h(X))|W]$ is known to smooth out features of $h$ in a neighborhood of $h_0$. In particular, Theorem (ref) requires no conditions (except the weak contraction requirement of Assumption (ref)$(ii)$) on the smoothing properties (or ill-posedness) of the mapping $h \rightarrow m(W,h)$ in a local neighborhood around $h_0$.} Intuitively, large values of $K_n$ induce greater sampling uncertainty but also act as a form of regularization by shrinking the Gaussian process in ((ref)). This regularization is crucial to control the ill-posedness in the model.

Theorem (ref) can be extended in several ways. One possibility is to consider a continuously updated version of the quasi-Bayes objective function. 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)$.\footnote{The usual continuously updated objective function takes $\widehat{\Sigma}(W,h)$ to be a suitable estimator of $$ \Sigma(W,h) = \{ \mathbb{E}[ \rho(Y,h(X)) \rho(Y,h(X))'|W ] \}^{-1} .$$} In this case, the associated quasi-Bayes posterior is

Another possible avenue, as in liao2011posterior, is to generalize the contraction in Theorem (ref) to settings where the unknown function $h_0$ is not uniquely identified from the data. In this setting, the identified set is given by $\Theta_0 = \{ h : \| m(W,h) \|_{L^2(\mathbb{P})} = 0\}$. Intuitively, regardless of point identification, draws from the quasi-Bayes posterior should concentrate in areas 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 extensions discussed above. To that end, we impose the following analog of Condition (ref).

Theorem (ref) shows that a continuously updated quasi-Bayes posterior is consistent, provided that at least one element of the identified set has sufficient regularity relative to the sample paths of the Gaussian process.

For the remainder of Section (ref), we focus on the setting where the true structural function $h_0$ is uniquely identified. Furthermore, unless otherwise specified, all further limit theory is developed with the quasi-Bayes posterior in ((ref)).\footnote{Extensions to the continuously updated version in ((ref)) can be handled in a similar manner to Theorem (ref). }

Contraction Rates

In this section, we further develop the quasi-Bayes limit theory. From the concluding statements in the proof of Theorem (ref), we can deduce that the quasi-Bayes posterior concentrates on bounded sample paths that lie in shrinking local neighborhoods around $h_0$. To be specific, there exists a sequence $\delta_n \rightarrow 0$ and $M > 0$ sufficiently large such that

As the focus in the preceding section was on general consistency, the analysis did not lead to an explicit form for $\delta_n$. In this section, we improve the preceding consistency results by explicitly quantifying the posterior contraction rate $\delta_n$. As we illustrate below, the posterior contraction rate depends on the interplay between $(i)$ Sample path realizations of the Gaussian process prior, $(ii)$ The local curvature of the pseudo-likelihood that defines the quasi-Bayes posterior, $(iii)$ the precise smoothing properties of the operator obtained from linearizing $h \rightarrow m(W,h)$ in a sufficiently regular local neighborhood around the true structural function $h_0$ and $(iv)$ the sequence of basis functions $b^K(W) = (b_1(W) , \dots , b_K(W))'$ used to compute a feasible estimate of the residual conditional mean function $m(W,h) = \mathbb{E}[ \rho(Y,h(X))|W]$.

Given the consistency result in ((ref)), to determine the rate of convergence, it suffices to restrict our attention to sample path realizations that are sufficiently regular and lie in a local neighborhood around $h_0$. Locally around $h_0$, we approximate the behavior of the map $h \rightarrow m(W,h)$ through a linearized version of it. Depending on the model and assumed hypothesis on the data generating process for $ \mathcal{Z} = (Y,X,W)$, there may be several distinct maps that serve as a reasonable linearization.\footnote{As the existence of such a map primarily serves as a proof technique, the precise choice is not crucial.} If the map $h \rightarrow m(W,h)$ is sufficiently regular around $h_0$, the natural linearization to consider is the Fr\'echet differential at $h_0$. This is the unique continuous linear operator $D_{h_0} : L^2(X) \rightarrow L^2(W)$ that satisfies

Intuitively, if $D_{h_0}[h]$ closely approximated $m(W,h)$ in local neighborhood around $h_0$, the smoothing properties of the map $h \rightarrow m(W,h)$ can be analyzed through the simpler linear map $h \rightarrow D_{h_0}[h]$. In the analysis that follows, we relate the smoothing property of $h \rightarrow D_{h_0}[h]$ to a change in regularity relative to the orthonormal basis $(e_i)_{i=1}^{\infty}$ that defines the Gaussian process in ((ref)). As the smoothness of $h_0$ is also defined relative to this basis through membership in the Sobolev ball ((ref)), this will allow us to study the action of $D_{h_0}$ on $(G_{\alpha},h_0)$ under a common regularity scale.

It will be convenient in our analysis to define a family of weak norms on $L^2(\mathcal{X})$, all obtained by shrinking the Fourier coefficients of a general function with respect to the basis $(e_i)_{i=1}^{\infty}$. To that end, we employ the following definition.

The following two conditions quantify the smoothing action of the map $h \rightarrow m(W,h)$ in a local neighborhood around $h_0$ by relating it to a weak norm in ((ref)). For $\alpha \in \mathcal{L}$ and $\gamma \in (0,\alpha)$ sufficiently small, denote a smooth local neighborhood around $h_0$ by

where $\mathcal{H}^{p}$ is defined as in ((ref)) for every $ p \geq 0$. The $\gamma$ in ((ref)) is used to account for the fact that a Gaussian process $G= G_{\alpha}$ as in ((ref)) does not possess exact regularity $\alpha$. The process has nearly exact regularity $\alpha$ in the sense that $ \mathbb{E} \big (\| G_{\alpha} \|_{\mathcal{H}^{\alpha}} \big) = \infty$ but $ \mathbb{E} \big (\| G_{\alpha} \|_{\mathcal{H}^{\beta}} \big) < \infty $ for every $\beta < \alpha$.

\begin{condition2}[Smoothing Link Condition] $(i)$ There exists $\epsilon \in (0,1), \gamma \in (0, \alpha)$ and a non-negative sequence of constants $ \sigma = (\sigma_i)_{i=1}^{\infty}$ such that for any $ M > 0 $, there exists constants $C_1(M),C_2(M) < \infty$ that satisfy $ \| D_{h_0}[h- h_0] \|_{L^2(\mathbb{P})} \leq C_1 \| h - h_0 \|_{w,\sigma} $ and $ \| h - h_0 \|_{w,\sigma} \leq C_2 \| D_{h_0}[h- h_0] \|_{L^2(\mathbb{P})} $ for every $h \in \Omega(M,\epsilon,\gamma) $. $(ii)$ The model is either mildly or severely ill-posed in the sense that $$ \sigma_{i} \asymp

$$ for some $R, \zeta \geq 0$. \end{condition2}

\begin{condition2}[Local Curvature] There exists $\epsilon \in (0,1), \gamma \in (0, \alpha)$ such that for any $ M > 0$, there exists a constant $ B = B(M) < \infty $ that satisfies $\| 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} Condition (ref) and (ref) are similar to Assumption $4.1$ and $5.2$ in chen2012estimation. Condition $\ref{lcurv}$ is trivially satisfied if $h \rightarrow m(W,h)$ is a linear map, as in NPIV (Example (ref)). If $D_{h_0}^*$ denotes the adjoint of $D_{h_0}$, Condition (ref) is satisfied if the self-adjoint operator $D_{h_0}^* D_{h_0}$ diagonalizes in the same eigenbasis $(e_i)_{i=1}^{\infty}$ that defines the Gaussian process in ((ref)). In this special case, Condition (ref) coincides with assumptions frequently used in the literature (e.g. \citealp*{knapik2011bayesian}, Assumption 3.1) for linear inverse problems. A stronger version of Condition (ref) is also used in \citet*{gugushvili2020bayesian} within the context of a white noise model with a known linear operator.

The quasi-Bayes objective function uses the estimator $\widehat{m}(w,h)$ in ((ref)) as a feasible analog for the true conditional mean function $m(W,h) = \mathbb{E}[ \rho(Y,h(X))|W]$. The discrepancy between $\widehat{m}(W,h)$ and $m(W,h)$ depends on two factors: $(i)$ the stochastic error between $\widehat{m}(W,h)$ and the population projection $\Pi_K m(W,h)$ and $(ii)$ the distance between the projection $\Pi_K m(W,h)$ and the true $m(W,h)$. While the effect of $(i)$ can be quantified\footnote{By Lemma (ref), under Condition (ref)$(i)$, the stochastic error has order $\sqrt{K}/ \sqrt{n}$ uniformly over any Sobolev ball $\mathbf{H}^t(M)$.} under very weak conditions for a large set of basis functions, the effect of $(ii)$ is intrinsic to the choice of first stage basis functions $b^K(W) = (b_1(W) , \dots , b_K(W))'$ used in the projection. Our final condition quantifies the approximation properties of the basis in a smooth local neighborhood around $h_0$.

\begin{condition2}[Basis Approximation] Let $\sigma_K$ be as in Condition (ref). There exists $\epsilon \in (0,1)$ and a sequence of (possibly diverging) constants $ \lambda_K < \infty$ such that for any $ M > 0 $, there exists $ B = B(M) < \infty $ that satisfies

for all $K $ sufficiently large and $\alpha \in \mathcal{L}$. \end{condition2} Condition (ref) is similar to assumptions commonly imposed in the literature.\footnote{For example, Corollary 5.1 of chen2012estimation imposes that the bias in the left side of ((ref)) is upper bounded by the stochastic projection error, having order $\sqrt{K}/ \sqrt{n}$. As the optimal $K$ there is chosen to balance the stochastic projection error and the projection bias $\sigma_{K+1} K^{-\alpha/d}$, this is equivalent to assuming $\lambda_K \lessapprox 1$.} Locally around $h_0$, the map $h \rightarrow m(W,h)$ acts as a smoothing operator that is similar to $D_{h_0}[h-h_0]$. The quantity $\sigma_{K+1} K^{-\alpha/d}$ represents the projection bias of the smoothed function.\footnote{This is the bias of $\| (\Pi_K - I) D_{h_0}[h-h_0] \| _{L^2(\mathbb{P})}$ if $b^K(W) = (b_1(W) , \dots , b_K(W))'$ are chosen to be the eigenfunctions generated by the singular value decomposition of $D_{h_0}$.} The sequence $\lambda_K$ in Condition (ref) represents a slack factor to account for rates that may be slightly larger than the usual local linear bias $\| (\Pi_K - I) D_{h_0}[h-h_0] \| _{L^2(\mathbb{P})}$, possibily due to a non optimal choice of first stage basis functions.

Given any Gaussian process $G_{\alpha}$ of regularity $\alpha$, let $ K_n (\alpha)$ denote a sequence of sieve dimensions defined by

Intuitively, $(K_n)_{n=1}^{\infty}$ denotes the optimal sieve dimension that balances the stochastic projection error and the local (around $h_0$) projection bias. In particular, if the model is mildly ill-posed, we have $ K_n \asymp n^{d/[2(\alpha + \zeta) + d]} $. If the model is severely ill-posed, we have $K_n \asymp (\log n)^{d/\zeta}$.

The following result provides contraction rates for the quasi-Bayes posterior in ((ref)).

As a point estimator for $h_0$, we consider the posterior mean

Given the posterior contraction rate in Theorem (ref), we expect that the posterior mean, as a point estimator, achieves a similar rate of convergence. Intuitively, this follows from the preceding result if the posterior tail probabilities appearing in ((ref)) and ((ref)) decay sufficiently fast. The following result verifies this.

The contraction rates for the mildly ill-posed case can be improved further if one includes the slack factor $\lambda_K$ from Condition (ref) as part of the projection bias in the definition of $K_n$ in ((ref)). If $ \nu_{n} = \lambda_{K_n}$ is asymptotically negligible relative to the other factors, the resulting rates will not change significantly from those stated (without $\nu_{n}$) in Theorem (ref) and Corollary (ref). For ease of exposition and notation, in the remaining sections we focus on the case where $\nu_{n} \lessapprox 1$.\footnote{It is straightforward to incorporate $\nu_{n} \uparrow \infty$ into the analysis, as in the preceding results.}

Inference

In this section, we establish the limiting quasi-posterior distribution for a class of linear functionals of $h$. We verify that, under suitable regularity conditions, the quasi-posterior distribution of a fixed linear functional is asymptotically Gaussian. Moreover, quasi-Bayesian credible sets have asymptotically exact frequentist coverage, provided that the quasi-Bayes posterior in ((ref)) is optimally weighted.

Given the posterior contraction rate in Theorem (ref), to derive the distributional limit theory, it suffices to restrict our analysis to a quasi-Bayes posterior with support contained on shrinking local neighborhoods around $h_0$. Let $\delta_n$ denote the posterior contraction rate in Theorem (ref). On route to deriving the posterior contraction rate, the previous results also provide contraction rates with respect to weaker metrics such as $ d_{w}(h,h_0) = \| m(W,h_{}) - m(W,h) \big] \|_{L^2(\mathbb{P})} $. It will be convenient in our analysis to emphasize this when defining the posterior support. To that end, define

Fix any $\gamma \in (0,\alpha)$ sufficiently small.\footnote{It is fine to let $\gamma \rightarrow 0$ slowly here. In fact, the posterior concentrates on subsets of $ \mathbf{H}^t(M) \cap \mathcal{H}^{\alpha - \gamma_n} (M) $, where $$ \gamma_n \asymp

$$} For constants $M,D > 0$ sufficiently large, we denote the localized posterior support and its image under the map $h \rightarrow m(W,h)$ by

We denote the entropy integral of the image of the localized posterior support by

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 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, the regularity of $h$ and $(iii)$ the convergence rate of $\| h- h_0 \|_{L^2(\mathbb{P})}$. Our main conditions on the localized support and remainder are as follows.

\begin{condition2} Let $\kappa,t$ be as in Assumption (ref). Suppose that

\end{condition2} Condition (ref)$(i-ii)$ arise as a consequence of empirical process tools that are used to control the deviation of $ \mathbb{E}_n \big[ \widehat{m}(W,h) ' \Sigma(W) \widehat{m}(W) \big]$ from its population analog $ \mathbb{E}[ \Pi_{K_n} m(W,h)' \Sigma(W) \Pi_{K_n} m(W,h)] $, uniformly over $h \in \Theta_n$. The dependence on $\kappa$ arises because the generalized residual function $\rho(.)$ in ((ref)) may be nonlinear and pointwise discontinuous in the input function $h$. As such, our argument instead utilizes the weaker $L^2(\mathbb{P})$ uniform continuity in Assumption (ref).

The Reproducing Kernel Hilbert Space (RKHS) of the Gaussian Process $G_{\alpha}$ in ((ref)) can be expressed as

Let $\mathbf{L} : L^2(X) \rightarrow \mathbb{R}$ denote a linear functional of interest. By duality (Riesz representation) there exists a function $\Phi(X) \in L^2(X)$ such that the linear functional can be expressed as

As such, to derive the limit theory with $\mathbf{L}(.)$, it suffices to work with $\Phi(.)$ directly.

In the preceding sections, the choice of weighting matrix $\widehat{\Sigma}(.)$ in the quasi-Bayes posterior ((ref)) did not influence the limit theory, provided that the eigenvalues of $\widehat{\Sigma}(.)$ are asymptotically bounded away from $0$ and $\infty$. Intuitively, if such a condition holds, rates of convergence can be determined by studying a quasi-Bayes posterior based on the simpler objective function $ h \rightarrow \mathbb{E}_n \big( \| \widehat{m}(W,h) \|_{\ell^2}^2 \big) $. However, to determine finer aspects of the posterior (such as a precise limiting distribution for functionals), it will be necessary to include the limiting behavior of $\widehat{\Sigma}(.)$ in the analysis. The following condition imposes that $\widehat\Sigma(W)$ converges to a limiting positive definite matrix $\Sigma(W)$. Properties of the quasi-posterior can then be analyzed through the deterministic matrix $\Sigma(.)$.

\begin{condition2} $(i)$There exists a positive definite matrix $\Sigma(.) $ such that $ \sup_{w \in \mathcal{W}} \| \widehat{\Sigma}(W) - \Sigma(W) \|_{op} = O_{\mathbb{P}}(\gamma_n) $ for some sequence $(\gamma_n)_{n=1}^{\infty}$ satisfying $\gamma_n \log(n) K_n \rightarrow 0$. $(ii)$ The eigenvalues of $\Sigma(W)$ are bounded away from $0$ and $\infty$: $\mathbb{P}(c \leq \lambda_{\min} (\Sigma(W)) \leq \lambda_{\max}(\Sigma(W)) \leq C ) = 1$ for some constants $ 0 < c \leq C < \infty$. \end{condition2} Given a positive definite matrix $\Sigma(.)$ as in Condition (ref), we can endow $L^2(W,\mathbb{R}^{d_{\rho}})$ with a norm (and inner product) by

This is the weighted norm induced by $\Sigma(.)$. As the quasi-Bayes objective function is based on a feasible version of ((ref)), it will be convenient in our analysis to view ((ref)) as the natural norm for $m(W,h)$.\footnote{Observe that, under Condition (ref)$(ii)$, this norm is equivalent to the usual unweighted norm $\|m \|_{L^2(\mathbb{P})} = \mathbb{E}[ \| m(W) \|_{\ell^2}^2]$. However, since certain operators such as adjoints are defined based on the specific form of the inner product, it is convenient to view $L^2(W,\mathbb{R}^{d_{\rho}})$ as a Hilbert space with natural norm as in ((ref)) directly.} Unless otherwise stated, for the remainder of this section, we view $L^2(W,\mathbb{R}^{d_{\rho}})$ as a Hilbert space with norm as in $(\ref{innerprod-induce})$. Denote the adjoint of $D_{h_0} : L^2(X) \rightarrow L^2(W,\mathbb{R}^{d_{\rho}})$ by $D_{h_0}^*$. We impose the following condition on the Riesz representer $\Phi(.)$ in ((ref)) that determines the linear functional $\mathbf{L}(.)$.

\begin{condition2} $(i)$ There exists $ \tilde{\Phi} \in \mathbb{H} $ such that $ \Phi = D_{h_0}^* D_{h_0} \tilde{\Phi}$. $(ii)$ $ D_{h_0}[\tilde{\Phi}]$ and $\Sigma(.)$ satisfy $ \sqrt{K_n} \sqrt{\log n} \|(\Pi_{K_n} - I)D_{h_0}[\tilde{\Phi}] \|_{L^2(\mathbb{P})} \rightarrow 0$ and $ \sqrt{K_n} \sqrt{\log n} \|(\Pi_{K_n} - I) \Sigma(W) D_{h_0}[\tilde{\Phi}] \|_{L^2(\mathbb{P})} \rightarrow 0$. \end{condition2} Condition (ref)$(i)$ is a source type condition on $\Phi$. This is similar to Condition 3.3 in monard2021statistical. We note that the condition for $\Phi$ to be in the range of the adjoint operator is a well known necessary condition for $\sqrt{n}$ estimability of linear functionals.\footnote{ For classical semiparametric models, this follows from van1991differentiable. For specific applications to NPIV and NPQIV setups, see severini2012efficiency and \citet*{chen2019penalized}.} Condition (ref)$(ii)$ imposes weak smoothness requirements on $D_{h_0}[\tilde{\Phi}]$ and $\Sigma(.) D_{h_0}[\tilde{\Phi}] $.

Let $ \rho $ denote any metric which metrizes weak convergence of probability measures on $\mathbb{R}$. If $(\mu_{n})_{n=1}^{\infty}$ is a random (data dependent) sequence of measures on $\mathbb{R}$, we say that $\mu_{n} \overset{\mathbb{P}}{\rightsquigarrow} \nu $ for some non-random limit measure $\nu$ if $\rho(\mu_n , \nu) \xrightarrow{\mathbb{P}} 0$. The following theorem shows that the quasi-posterior distribution of a sufficiently smooth (in the sense of Condition (ref)) linear functional $\mathbf{L}(.)$ can be well approximated by a suitable Gaussian measure.

Observe that the two variances in Theorem (ref) agree if and only if the limiting weighting matrix of the quasi-Bayes objective is the optimal weighting matrix. That is, when

In this case, the common limiting variance is given by

One implication of an optimally weighted quasi-Bayes objective function, or equivalently the equality of variances in Theorem (ref), is that quasi-Bayesian credible bands centered around the posterior mean have asymptotically exact frequentist coverage. To be specific, given a linear functional $\mathbf{L}(.)$ and significance level $\gamma \in (0,1)$, define

Denote the quasi-Bayesian credible band at significance level $\alpha$ by

An immediate consequence of Theorem (ref) is that optimally-weighted quasi-Bayes credible sets have asymptotically exact frequentist coverage.

Proofs

In this section, we provide proofs for all the main results. We first state some notation that appears throughout the proofs. We denote by $\widehat{G}_{b,K}^{o}$ the matrix

where $G_{b,K} = \mathbb{E} \big( [b^K(W)] [ b^K(W) ]' \big)$. We denote by $\zeta_{b,K}$ the quantity $ \zeta_{b,K} = \sup_{w \in \mathcal{W} } \| G_{b,K}^{-1/2} b^K(w) \|_{\ell^2} $. By Condition (ref), we have $\zeta_{b,K} \lessapprox \sqrt{K}$. For ease of notation, in several of the limit results, we often suppress the dependence of $K=K_n$ on $n$.