arXiv:1607.07343v1 [math.ST] 25 Jul 2016

Gaussian processes and Bayesian moment
estimation∗
Jean-Pierre Florens†

Anna Simoni‡

Toulouse School of Economics

CREST and Ensae

This version: July 2016
Abstract
Given a set of estimating equations (EE) that characterize a parameter θ,
we investigate a semiparametric Bayesian approach for inference on θ that does
not restrict the data distribution F apart from the EE. As main contribution,
we construct a degenerate Gaussian process prior that, conditionally on θ, restricts the F generated by this prior to satisfy the EE with probability one.
Our prior works even in the more involved case where the number of EE is
larger than the dimension of θ. We show that this prior is computationally
convenient. Since the likelihood function is not specified by the model, we
approximate it based on a linear functional transformation of F that has an
asymptotically Gaussian empirical counterpart. This likelihood is used to construct the posterior distribution. We provide a frequentist validation of our
procedure by showing consistency and asymptotic normality of the posterior
distribution of θ.
Key words: Estimating Equations, Gaussian processes, overidentification, invariance, posterior consistency, functional equation.
∗

First version: February 2012. The authors especially thank Yuichi Kitamura for insightful discussions. The authors are grateful to seminars and conferences in: Berlin, Boston College, Bristol,
Carlos III, CREST, Northwestern, SBIES 2015, ICEEE 2015, NASM 2012, Toulouse. We thank
financial support from ANR-13-BSH1-0004 (IPANEMA). Anna Simoni gratefully acknowledges financial support from Labex ECODEC (ANR - 11-LABEX-0047), SFB-884 and hospitality from the
University of Mannheim.
†
Toulouse School of Economics - 21, allée de Brienne - 31000 Toulouse (France). Email: jeanpierre.florens@tse-fr.eu
‡
CREST, 15 Boulevard Gabriel Péri, 92240 Malakoff (France). Email: simoni.anna@gmail.com
(corresponding author ).

1

1

Introduction
Statistical models are often formulated via estimating equations (EE) - or moment

restrictions - of the form EF [h(θ, X)] = 0, where h(θ, X) is a vector-valued function of an
observable random element X with distribution F and a parameter vector θ. These EE
provide the only information available about θ and the data distribution. Given a set of
EE, this paper builds a semiparametric Bayesian inference procedure for θ that imposes
these moment conditions in the nonparametric prior distribution for the data distribution
F and that is computationally convenient. Apart from these moment restrictions, F is left
unrestricted.
A main advantage of Bayesian inference consists in providing a well-defined posterior
distribution that is important for many decision problems and for predictive analysis. On
the other hand, constructing Bayesian inference procedures for moment condition models
presents two difficulties. A first difficulty is due to the fact that a likelihood is not available.
A second difficulty arises when the number of EE exceeds the dimension of θ (overidentification) because imposing overidentifying moment restrictions on the prior distribution for
the nonparametric F is challenging. The contribution of this paper is to propose an elegant
approach that allows to deal with these two difficulties. As a by-product we show that
the quasi-likelihood of some Laplace-type procedures arises as the limit of our Bayesian
procedure when our prior becomes diffuse.
The model we consider is as follows. Let X be an observable random element in Rm
with distribution F and X1 , . . . , Xn be an i.i.d. sample of X. The parameter θ ∈ Θ ⊂ Rp
is linked to the data generating process (DGP) F through the EE
EF [h(θ, X)] = 0,

(1.1)

where h(θ, x) := (h1 (θ, x), . . . , hd (θ, x))T and the functions hj (θ, x), j = 1, . . . , d are realvalued and known up to θ. We assume d ≥ p and our main interest is the case where
d > p, which is in general more challenging than the case d = p. Examples of the case
d > p can be found in Broniatowski and Keziou [2012]. Apart from (1.1), F is completely
unrestricted. The Bayesian procedure proposed in this paper constructs a nonparametric
prior distribution for F with support equal to the subset of distributions that satisfy the
EE for a given θ.
Imposing moment restrictions via semiparametric priors may be challenging depending
on the relationship existing between θ and F . More precisely, when the model is exactly
identified (i.e. p = d) and (1.1) characterizes θ as an explicit function of F , say θ = b(F )
for some function b(·), then (1.1) does not restrict F . Therefore, one places an unrestricted
nonparametric prior on F and recovers the prior of θ via the transformation θ = b(F ). The
(θ, F )s generated by this prior automatically satisfy the constraints.

2

On the contrary, when d > p (overidentified model), θ cannot be expressed as an explicit function of F . Indeed, (1.1) imposes constraints on F and the existence of a solution
θ to (1.1) is guaranteed only for a subset of distributions F . Therefore, a restricted nonparametric prior on F must be specified conditionally on θ and the support of this prior is
a proper subset of the set of probability distributions. It turns out that incorporating overidentifying moment restrictions in a semiparametric prior for (θ, F ) is not straightforward.
In this paper we propose a way to construct a semiparametric prior that incorporates the
overidentifying EE.
Our strategy is based on a degenerate Gaussian process (GP) prior with restricted support which is easy to deal with and that works as follows. The DGP F is assumed to admit
a density function f with respect to some positive measure Π chosen by the researcher
(for instance the Lebesgue measure). Then, we endow f with a GP prior conditional on
θ. The d > p EE are incorporated by constraining the prior mean and prior covariance

of this GP in an appropriate way. Because this prior imposes the moment restrictions, it

is degenerate on a proper subset of the set of probability density functions. The reason

for the appropriateness of a GP prior in such a framework is due to the fact that the EE
in (1.1) are linear in f and the linearity of the model matches extremely well with a GP

prior. A remarkable feature of our method is that the EE are imposed directly through the
GP prior of f given θ without requiring a second step projection over the set of density
functions satisfying the moment restrictions. To the best of our knowledge a GP prior has
not been used yet in the EE framework.

Our Bayesian procedure, that we call the GP-approach, is constructed as follows. In

the overidentified case (d > p) we first specify a prior on θ and then a GP prior on f conditional on θ. In the case d = p we may either proceed as in the overidentified case or specify
an unrestricted GP prior on f and then deduce from it the prior for θ through the explicit
relationship θ = b(f ). We circumvent the difficulty of the likelihood function specification,

which is not available, by constructing a linear functional transformation of the DGP F
such that its empirical counterpart, say rn , has an asymptotic Gaussian distribution. This
will be used as the sampling model. Therefore, our model is approximately conjugate and
allows easy computations while being nonparametric in F .
We provide a closed-form expression for the marginal posterior distribution of θ (obtained by integrating out f ) and propose the maximum of this distribution as an estimator
for θ. The maximum a posteriori of θ is usually not available in closed-form but can be
easily computed via drawn from the marginal posterior. We show that the quasi-likelihood
function – also called limited information likelihood (LIL) – used, among others, by Kim
[2002] and Chernozhukov and Hong [2003], can be obtained as the limit of our marginal
likelihood when the GP prior for f is allowed to become diffuse. In addition, when the prior

for f becomes noninformative, the marginal posterior distribution for θ becomes the same

3

(up to constants) as the GEL objective function with quadratic criterion and is a monotonic
transformation of the continuous updating GMM objective function (Hansen et al. [1996]).
Finally, we provide a frequentist validation of our method by showing: (i) posterior
consistency, (ii) frequentist consistency of the maximum a posteriori estimator, and (iii)
asymptotic normality of the posterior distribution of θ.

Related literature. Inference in an EE framework has received very much attention
in the past literature. Among the most popular statistical methods in this framework is
the Empirical Likelihood (EL)-based inference (e.g. Owen [1988], Qin and Lawless [1994],
Imbens, Owen [2001], Chen and Van Keilegom [2009]) and its variations like the Exponential Tilting (e.g. Kitamura and Stutzer [1997] and Kitamura [1997]), the Exponentially
Tilted EL (ETEL) proposed by Schennach [2007], and the Generalized EL (GEL) (e.g.
Smith [1997], Newey and Smith [2004]). EL and ETEL can also be used in a Bayesian
framework for valid inference as established by Lazar [2003], Grendar and Judge [2009],
Schennach [2005], Rao and Wu [2010], Chib et al. [2016] and Chaudhuri et al. [2017].
An alternative frequentist inference procedure, which is very popular especially in
econometrics, is the Generalized method of Moments (GMM) (Hansen [1982] and Hansen and Singleton
[1982]). From a Bayesian inference point of view, alternative procedures that have been
proposed in the literature can be classified in two types depending on whether the moment
conditions (1.1) are imposed in the likelihood or in the nonparametric prior. Procedures
of the first type construct a quasi-likelihood by exponentiating the quadratic criterion (or
the GEL criterion) associated with the empirical counterpart of (1.1) and include Kwan
[1999], Kim [2002], Chernozhukov and Hong [2003] and Liao and Jiang [2011] among others. Our paper shows that the quasi-likelihood used in this type of approach arises as the
limit of our GP prior as it becomes diffuse. We provide thus a fully Bayesian justification

to this approach. Procedures of the second type impose the moment conditions in the prior
for (θ, F ) while leaving the likelihood completely unrestricted and include, in addition to
the Bayesian ET and ETEL discussed above, Chamberlain and Imbens [2003] who use a
Dirichlet prior, Kitamura and Otsu [2011] and Shin [2014] who propose a two-step procedure based on a projection of a Dirichlet process mixture prior and of a mixture of Dirichlet
Process prior, respectively, and Bornn et al. [2015] who use Hausdorff measures to build
probability tools for dealing with moment estimation.
The paper is organized as follows. The GP-approach is described in section 2, which

contains our main contribution. Here, we also show the link existing between our approach
and some frequentist approaches in an EE framework. In section 3 we analyze asymptotic
properties of the posterior distribution of θ and of the maximum a posteriori estimator.
In section 4 we detail how to implement our method through simulation studies. All the
proofs are gathered in the Appendix which is contained in the Supplementary Material.

4

2

The Gaussian Process (GP)-approach
Let X be a continuous random element in S ⊆ Rm with distribution F and X1 , . . . , Xn

be an i.i.d. sample of X. Assume that F is absolutely continuous with respect to some
positive measure Π (e.g. the Lebesgue measure) with density function f . In other words,
conditionally on f the data are drawn from F : X1 , . . . , Xn |f ∼ F . The set of probability
density functions on S with respect to Π is denoted by M .

Let θ ∈ Θ ⊆ Rp be the parameter of interest characterized by (1.1). By adopting a

frequentist point of view, we denote, throughout the paper, the true value of θ by θ∗ , the

true DGP by F∗ and its density with respect to Π by f∗ . The model is assumed to be
well-specified, that is, there exists θ∗ ∈ Θ such that EF∗ (h(θ∗ , X)) = 0 holds. We endow

S ⊆ Rm with the trace of the Borelian σ-field BS and specify Π as a positive measure on

this subset. We denote by E = L2 (S, BS , Π) the Hilbert space of square integrable functions

on S with respect to Π and by BE the Borel σ-field generated by the open sets of E. The

scalar product and norm on this space are defined in the usual way and denoted by h·, ·i
and || · ||, respectively.

The parameters of the model are (θ, f ), where f is the nuisance parameter, and the

parameter space is
Λ=



(θ, f ) ∈ Θ × EM ;

Z


h(θ, x)f (x)Π(dx) = 0 ,

EM := E ∩ M,

where h : Θ × Rm → Rd is a known function. In the following of the paper we maintain the
following assumption.

Assumption 2.1. (i) The true f∗ satisfies f∗ ∈ EM := E ∩ M ; (ii) the moment function
h(θ, ·) is such that hi (θ, ·) ∈ E for every i = 1, . . . , d and for every θ ∈ Θ, where hi denotes

the i-th component of h; (iii) d ≥ p.

Assumption 2.1 (i) restricts f∗ to be square integrable with respect to Π and is for

instance verified if f∗ is bounded and Π is a bounded measure. The model is made up of

three elements that we detail in the next two subsections: a prior on θ, denoted by µ(θ), a
conditional prior on f given θ, denoted by µ(f |θ), and the sampling model. In the following,

we shorten “almost surely” by “a.s.” and omit the probability which “a.s.” refers to. We
denote by EF the expectation taken with respect to F and by E∗ the expectation taken
with respect to F∗ .

2.1

Prior distribution

We specify a prior probability measure µ for (θ, f ) of the form µ(θ, f ) = µ(θ)µ(f |θ). By

abuse of notation, µ(θ) will also denote the Lebesgue density of the prior distribution of θ

5

in the case it admits it. The prior µ(θ) may either be flat (non-informative) or incorporate
any additional information available to the econometrician about θ. In any case, it is tacitly
assumed that µ(θ) is such that the posterior of θ exists.
Given a value for θ, the conditional prior µ(f |θ) is specified such that its support equals

the subset of functions in E that integrate to one and satisfy (1.1) for this particular value
of θ. At the best of our knowledge, the construction of such a conditional prior µ(f |θ) is

new in the literature and we now explain it in detail.

Construction of the conditional prior µ(f |θ). We construct the conditional prior

distribution µ(f |θ) of f , given θ, as a GP on BE with mean function f0θ ∈ EM and covariance

operator Ω0θ : E → E. We restrict f0θ and Ω0θ to guarantee that the trajectories f generated

by µ(f |θ) are such that the corresponding F (which is such that dF = f dΠ) integrates to 1
and satisfies equation (1.1) with probability 1. The two sets of restrictions that we impose
are the following (one on f0θ and one on Ω0θ ):
Restriction 1 (Restriction on f0θ ). The prior mean function f0θ ∈ EM is chosen such that
Z

h(θ, x)f0θ (x)Π(dx) = 0.

(2.1)

Restriction 2 (Restriction on Ω0θ ). The prior covariance operator Ω0θ : E → E is chosen

such that

(

1/2

Ω0θ h(θ, x) = 0
1/2

Ω0θ 1

1/2

(2.2)

= 0
1/2

1/2

where Ω0θ : E → E denotes the positive square root of Ω0θ : Ω0θ = Ω0θ Ω0θ .
The covariance operator Ω0θ is linear, self-adjoint and trace-class.1 Due to Restriction
2, Ω0θ is not injective. In fact, the null space of Ω0θ , denoted by N(Ω0θ ), is not trivial and
contains effectively the constant 1 – which implies that the trajectories f generated by the
prior integrate to 1 a.s. (with respect to Π) – and the functions h(θ, x) – which implies that
the trajectories f satisfy the moment conditions a.s. This means that Ω0θ is degenerate
in the directions along which we want that the corresponding projections of f and f0θ are
equal. Therefore, the support of µ(f |θ) is a proper subset of E. This is the meaning of the

next lemma.

Lemma 2.1. The conditional GP prior µ(f |θ), with mean function f0θ and covariance

operator Ω0θ satisfying Restrictions 1 and 2, generates trajectories f that satisfy µ(f |θ)a.s. the conditions
1

A trace-class operator is a compact operator with eigenvalues that
R are summable. Remark that
this guarantees that the trajectories f generated by µ(f |θ) satisfy f 2 dΠ < ∞ a.s.

6

Z

f (x)Π(dx) = 1

Z

and

h(θ, x)f (x)Π(dx) = 0.

Remark 2.1. Restrictions 1 and 2 imply that the trajectories generated by µ(f |θ) in-

tegrates to 1 (with respect to Π) and satisfy (1.1) a.s. but they do not guarantee nonnegativity of the trajectories. Thus, the support of µ(f |θ) is smaller than E but bigger

than EM . To impose non-negativity one could: (i) either project the prior on the space of
R
non-negative functions, or (ii) write f = g2 or f = eg / eg , g ∈ E, and specify a conditional
prior distribution for g, given θ, instead of for f . Nonetheless, it is important to notice

that: the projected prior in (i) is not Gaussian anymore, and in (ii) we cannot use our
restricted GP prior because if it is specified on g it does not work to impose the moment

restrictions on f so that a different restricted prior should be constructed. Moreover, in
both (i) and (ii) the resulting posterior for θ is not available in closed form which is instead

one of the main advantages of our procedure. Therefore, it is not possible to impose the
non-negativity constraint if one wants to use our restricted GP prior. However, because

our goal is to make inference on θ while f is a nuisance parameter, failing to impose the
non-negativity constraint is not an issue as long as our procedure is shown to be consistent
for θ (which we show in section 3).
From a practical implementation point of view, a covariance operator satisfying Restriction 2 and a f0θ satisfying Restriction 1 may be constructed as follows.
Construction of Ω0θ . Let (λj )j∈N be a decreasing sequence of non-negative numbers such
P
that
j λj < ∞, and (ϕj )j∈N be an orthonormal basis (o.n.b.) for E. Then, ∀φ ∈ E:
P
Ω0θ φ = ∞
j=0 λj hφ, ϕj iϕj . Remark that (λj )j∈N and (ϕj )j∈N correspond to the eigenvalues

and eigenfunctions of Ω0θ , respectively. Since the null space N(Ω0θ ) ⊂ E is spanned by
{1, h1 (θ, ·), . . . , hd (θ, ·)}, we can set the first eigenfunctions of Ω0θ equal to the elements of
any o.n.b. of N(Ω0θ ). Restriction 2 is then fulfilled by setting the corresponding eigenvalues

equal to 0. For instance, if {1, h1 (θ, ·), . . . , hd (θ, ·)} are orthonormal as elements of E, then

N(Ω0θ ) has dimension d + 1, the first eigenfunctions are (ϕ0 , ϕ1 , . . . , ϕd )T = (1, hT )T and
the corresponding eigenvalues are λj = 0, ∀j = 0, 1, . . . , d. Remark that in this case, necesR
sarily, Π(dx) = 1. If {1, h1 (θ, ·), . . . , hd (θ, ·)} are not orthonormal then one can use their

orthonormalized counterparts as the first eigenfunctions of Ω0θ . The latter is the method we
use to implement our procedure. The remaining components (ϕj )j>d are chosen such that
P
(ϕj )j≥0 forms an orthonormal basis of E and (λj )j>d are chosen such that j>d λj < ∞.

Hence,

∀φ ∈ E,

Ω0θ φ =

∞
X

j=d+1

λj hφ, ϕj iϕj

where we suppress the dependence of ϕj on θ for simplicity. Examples of choices for (λj )j>d
are, for some constant c > 0: (i) λj = cj −a with a > 1, (ii) λj = ce−j . In section 4 we

7

clarify the construction of Ω0θ .
Construction of f0θ . By definition, f0θ (·) ∈ E is orthogonal to h(θ, ·). Therefore, if

the functions {ϕj }j>d used to construct Ω0θ are bounded in absolute value on S with

bounds {Bj }j>d , respectively, then we can select a sequence of constants {aj }j>d that
P
P
satisfy j>d |aj |Bj ≤ 1 and construct f0θ as f0θ (x) = 1 + j>d aj ϕj (x). Alternatively, if

not all the {ϕj }j>d are bounded we can set all but a finite number of the {aj }j>d to zero.

Remark 2.2. In the just-identified case where d = p and the EE (1.1) can be solved
explicitly for θ (that is, θ = b(f ), for some functional b), then the prior for (θ, f ) may be
constructed in an alternative way: one can first specify an unrestricted prior for f and then
recover from it the prior for θ. When b is a linear functional and θ can take any value in Rp ,
one can specify a GP prior µ(f ) for f (independent of θ) with a mean function f0 restricted
1/2

only to be a pdf and a covariance operator Ω0 restricted only to satisfy Ω0 1 = 0. Then, the
prior for θ is obtained through the transformation b(·) and will be Gaussian. Because the
support of this prior is Rp , then this approach is feasible if every value in Rp is plausible
for θ. For example, if θ = EF (x) and the support of x is Rp , then b(f ) = hf, ιi and

µ(θ) = N (hf0 , ιi, hΩ0 ι, ιi), where ι ∈ E denotes the identity function ι(x) = x.

2.2

The sampling model

Given the observed i.i.d. sample (x1 , . . . , xn ), the likelihood function is

Qn

i=1 f (xi ).

While apparently simple, using this likelihood for Bayesian inference on θ makes the analysis of the posterior distribution complicated. This is because to compute the posterior
for θ one has to marginalize out f . Since a GP prior is not a natural conjugate of the
i.i.d. model then, marginalization of f has to be carried out through numerical, or Monte
Carlo, integration on a functional space, which may be computationally costly. To avoid
this difficulty, we propose an alternative and original way to construct the sampling model
that allows for a conjugate analysis and prevents from numerical integration. Our approach
is based on a functional transformation rn of the sample X1 , . . . , Xn .
This transformation rn is chosen by the researcher and must have the following characteristics: (I) rn is an observable element of an infinite-dimensional Hilbert space F (to be
defined below), for instance a L2 -space; (II) rn converges weakly towards a Gaussian process in F; (III) the expectation of rn , conditional on f , defines a linear operator K : E → F

such that EF (rn ) = Kf ; (IV) rn is a one-to-one transformation of some sufficient statistic. Moreover, rn ∈ F is a Hilbert space-valued random variable (H-r.v.). We recall that,

for a complete probability space (Z, Z, P), rn is a H-r.v. if it defines a measurable map

rn : (Z, Z, P) → (F, BF ), where BF denotes the Borel σ-field generated by the open sets

of F.

8

Construction of rn . Let T ⊆ Rl , l > 0. To construct rn we first select a function

k(t, x) : T × S → R (or in C) that is measurable in x for every t ∈ T and that is nonconstant in (t, x). The transformation rn is then taken to be the expectation of k(t, ·)

under the empirical measure:

n

rn (t) =

1X
k(t, xi ),
n
i=1

∀t ∈ T.

Define F = L2 (T, BT , ρ) where ρ is a measure on T and BT denotes the Borel σ-field

generated by the open sets of T. The scalar product and norm on F are defined in the

usual way and denoted by h·, ·i and k·k, respectively, with the same notation as for the inner

product and norm in E. The function k(t, x) defines also a bounded operator K : E → F

and must be such that, for every ϕ ∈ E, Kϕ ∈ F and rn is an H-r.v. with realizations in
F. Hence,

K:E

→ F
R
ϕ 7→
k(t, x)ϕ(x)Π(dx).

(2.3)

For every f ∈ EM , Kf is the expectation of k(t, X) under F : (Kf )(t) = EF (k(t, X)).

Under the true distribution F ∗ the expectation of rn is Kf∗ and the covariance function of

rn is: ∀s, t ∈ T,

1
1
σ(t, s) = E∗ rn (t)rn (s) = [E∗ (k(t, X)k(s, X)) − E∗ (k(t, X))E∗ (k(s, X))] .
n
n
If the class of functions {k(t, ·), t ∈ T} is Donsker then, as n → ∞, the conditional distri√
bution of n(rn − Kf∗ ) weakly converges to a GP with covariance operator Σ : F → F
defined as

∀ψ ∈ F,

(Σψ)(t) =

Z

σ(t, s)ψ(s)ρ(ds)

(2.4)

which is one-to-one, linear, positive definite, self-adjoint and trace-class. In the following
we assume that {k(t, ·), t ∈ T} is Donsker such that rn is approximately Gaussian: rn ∼

GP(Kf∗ , Σn ) where Σn =

1
n Σ.

Among all the functions k(t, x) that satisfy the previous

assumptions, one should keep only the ones such that the corresponding rn is a one-to-one
transformation of some sufficient statistic, as required in (IV) above. This, which will be
tacitly assumed in the following of the paper, guarantees that the posterior distribution
computed by using rn does not depend on the particular choice of k(t, x).
In our analysis we treat f∗ as a realization of the random parameter f and Σn as

known. Therefore, the sampling distribution of rn |f is P f = GP(Kf, Σn ) and we construct

the posterior distribution based on it. In practice, Σn must be replaced by its empirical

9

counterpart. In finite sample, P f is an approximation of the true sampling distribution
but the approximation error vanishes as n → ∞. Moreover, the approximating sampling

distribution P f is only used to construct the posterior distribution and the proofs of our
asymptotic results do not rely on it. We give in the following two examples where rn is a

one-to-one transformation of sufficient statistics and K is injective.
Example 2.1 (Empirical cumulative distribution function (cdf)). Let (X1 , . . . , Xn ) be an
i.i.d. sample of X ∈ R. A possible choice for k(t, x) is k(t, x) = 1{x ≤ t}, where 1{A}
P
denotes the indicator function of the event A. In this case, rn (t) = Fn (t) := n1 ni=1 1{Xi ≤
R
t} is the empirical cdf and the operator K is (Kϕ)(t) = S 1{s ≤ t}ϕ(s)Π(ds), ∀ϕ ∈ E. By

the Donsker’s theorem, Fn (·) is asymptotically Gaussian with mean the true cdf F∗ (·) and

covariance operator characterized by the kernel:

1
n (F∗ (s

∧ t) − F∗ (s)F∗ (t)).

Example 2.2 (Empirical characteristic function). Let (x1 , . . . , xn ) be an i.i.d. sample of
P
x ∈ R. Let k(t, x) = eitx , so that rn (t) = cn (t) := n1 nj=1 eitxj is the empirical charR
acteristic function. In this case, the operator K is (Kϕ)(t) = S eits ϕ(s)Π(ds), ∀ϕ ∈ E.

By the Donsker’s theorem, cn (·) is asymptotically a Gaussian process with mean the true

characteristic function c(·) := E∗ [eitx ] and covariance operator characterized by the kernel:
1
n (c(s

+ t) − c(s)c(t)).

The following lemma gives an useful characterization of the operator Σ in terms of
K and its adjoint K ∗ . We recall that the adjoint K ∗ of a bounded and linear operator
K : E → F is defined as the operator from F to E that satisfies hKϕ, ψi = hϕ, K ∗ ψi,

∀ϕ ∈ E and ∀ψ ∈ F. In our case, an elementary computation shows that (K ∗ ψ)(t) =
R
T k(t, x)ψ(t)ρ(dt), ∀ψ ∈ F.

Lemma 2.2. Let K : E → F be a bounded and linear operator defined as in (2.3) and
R
K ∗ : F → E be its adjoint, that is, (K ∗ ψ)(t) = T k(t, x)ψ(t)ρ(dt), ∀ψ ∈ F. The operator

Σn = n1 Σ, with Σ : F → F defined in (2.4) takes the form
∀ψ ∈ F,

Σψ = KMf K ∗ ψ − (KMf 1)hMf , K ∗ ψi

(2.5)

where Mf : E → E is the multiplication operator Mf ϕ = f∗ ϕ, ∀ϕ ∈ E.
We denote by D the subset of E whose elements integrate to 0 with respect to Π:


Z
g(x)Π(dx) = 0
D := g ∈ E;
S

and denote by K|D the operator K restricted to D ⊂ E. Remark that D contains the
1/2

functions in E that are the difference of pdf s of F with respect to Π. Moreover, R(Ω0θ ) ⊂
D, where R(·) denotes the range of an operator, and

10

1/2
R(Ω0θ )



Z
Z
ϕ(x)Π(dx) = 0 .
ϕ(h)h(θ, x)Π(dx) = 0 and
⊆ ϕ ∈ E;
S

S

Remark that the equality holds when the sequence (λj )j>d , used to construct Ω0θ , is strictly
positive. We also denote by D∗ the subset of L2 (S, BS , F∗ ) whose elements integrate to 0

R
with respect to F∗ : D∗ := g ∈ L2 (S, BS , F∗ ); S g(x)F∗ (dx) = 0 . It follows from Lemma

2.2 that we can write Σ = HH ∗ where

H∗ : F

H : D∗ → F
ϕ

7→ Kf∗ ϕ

ψ

→ D∗

7→ K ∗ ψ − hKf∗ , ψi.

(2.6)

Remark that, H ∗ is the adjoint of H and is the projection of the adjoint of Kf∗ (determined

by using the scalar product in L2 (S, BS , F∗ )) onto D∗ : H ∗ ψ = PK ∗ ψ, where Pϕ = ϕ −
R
2
S ϕ(x)F∗ (dx) is the projection operator of L (S, BS , F∗ ) on D∗ .

2.3

Posterior distribution

The Bayesian model defines a joint distribution on (the Borel σ-field) of Λ and can be
summarized in the following way:
θ ∼ µ(θ)
f |θ ∼ µ(f |θ) = GP(f0θ , Ω0θ ),

Z

1

h(θ, x)f0θ (x)Π(dx) = 0 and

2
(1, h(θ, ·)T )T = 0
Ω0θ

rn |f, θ ∼ rn |f ∼ P f = GP(Kf, Σn )

(2.7)

where we use the GP approximation P f . Theorem 1 in Florens and Simoni [2012] shows

that the joint distribution of (f, rn ), conditional on θ, is:
f
rn

!

θ ∼ GP

f0θ
Kf0θ

!

,

Ω0θ ,

Ω0θ K ∗

KΩ0θ , Σn + KΩ0θ K ∗

!!

(2.8)

where (Σn + KΩ0θ K ∗ ) : F → F, Ω0θ K ∗ : F → E and KΩ0θ : E → F. The marginal

sampling distribution of rn conditional on θ, obtained by integrating out f , is:
rn |θ ∼ Pnθ := GP(Kf0θ , Σn + KΩ0θ K ∗ ).

(2.9)

We now discuss the posterior distribution of the parameter of interest θ, denoted by µ(θ|rn ).
The conditional posterior distribution µ(f |rn , θ) of the nuisance parameter f , given θ, and

its properties can be obtained by using results in Florens and Simoni [2014] and we briefly
discuss it in Appendix ??.

11

2.3.1

Posterior distribution of θ

The marginal posterior for θ, denoted by µ(θ|rn ), is obtained by using the marginal
sampling distribution Pnθ given in (2.9). We first have to characterize the likelihood of Pnθ
with respect to an appropriate common dominating measure that will be denoted by Pn0 .
The following theorem characterizes a probability measure Pn0 which is equivalent to Pnθ as
well as the corresponding likelihood of Pnθ with respect to Pn0 . Denote by H −1 : R(H) → D∗

(resp. (H ∗ )−1 : D∗ ∩ R(K ∗ ) → F) the left inverse of H (resp. right inverse of H ∗ ) which

takes the form H −1 · = f∗−1 K −1 · (resp. (H ∗ )−1 · = (K ∗ )−1 ·). Denote by (ljθ , ρj (θ), ψj (θ))j≥0
1/2

the singular value decomposition of the operator H −1 KΩ0θ which is well defined by Lemma
?? in the Appendix if K|D is injective and f∗ is bounded away from zero on S.
Theorem 2.1. Let Pn0 be a Gaussian measure with mean Kf∗ and covariance operator

n−1 Σ, i.e. Pn0 = GP(Kf∗ , n−1 Σ) with Σ defined in (2.4). For n fixed, if K|D is injective

and f∗ is bounded away from zero on S, then Pn0 and Pnθ are equivalent. Moreover, assume

that ∀j ≥ 0 and ∀θ ∈ Θ, ψj (θ) ∈ D∗ ∩ R(K ∗ ). Then, the Radon-Nikodym derivative is
given by

dPnθ
(rn )
(2.10)
dPn0


s
2
√
∞
∞
 1X
Y
Zj − h nK(f0θ − f∗ ), (H ∗ )−1 ψj (θ)i  { 1 kZk2 }
n−1
=
e 2 Σ
2 exp − 2
2

n−1 + ljθ
1 + nljθ
j=0
j=0

pnθ (rn ; θ) :=

where Z :=

√

n(rn − Kf∗ ), Zj := hZ, (H ∗ )−1 ψj (θ)i for all j ≥ 0, and kZkΣ := kΣ−1/2 Zk.

The quantity Σ−1/2 denotes the inverse of the positive square root of Σ and kΣ−1/2 Zk2
P
−2
2
2
∞
is defined as the limit in F of the series m
j=0 σj hZj , φj i as m → ∞ (where {σj , φj }j=0

is the eigensystem of Σ). By using (2.10), the (marginal) posterior distribution of θ takes
the form (after simplifying the terms that do not depend on θ):

µ(θ|rn ) =

=

R

pnθ (rn ; θ)µ(θ)
(2.11)
Θ pnθ (rn ; θ)µ(θ)dθ


√
2
∗ −1
Q∞ q 1
1 P∞ (Zj −h nK(f0θ −f∗ ),(H ) ψj (θ)i)
− 2 j=0
µ(θ)
2 exp
2
j=0
n−1 +ljθ
1+nljθ


√
2
∗ −1
Q∞ q 1
1 P∞ (Zj −h nK(f0θ −f∗ ),(H ) ψj (θ)i)
exp
−
µ(θ)dθ
j=0
j=0
2
n−1 +l2
1+nl2

R

Θ

jθ

jθ

and can be used to compute a point estimator of θ. We propose to use the maximum a
posterior (MAP) estimator θn defined as

12

θn := arg max µ(θ|rn )
θ∈Θ

= arg max
θ∈Θ

j=0

= arg max
θ∈Θ

∞
Y

∞
Y

j=0

s

n−1

1
2
+ ljθ

s

n−1

1
2
+ ljθ


2 
√
∞
 1X
∗
−1
Zj − h nK(f0θ − f∗ ), (H ) ψj (θ)i 
exp −
µ(θ)
2
 2

1 + nljθ
j=0


2 
√
∞
 1X
∗
−1
h n(rn − Kf0θ ), (H ) ψj (θ)i
exp −
µ(θ)
2
 2

1 + nljθ
j=0

or the posterior mean estimator E(θ|rn ) :=

(2.12)

R

Θ

θµ(θ|rn)dθ.

Remark 2.3. As discussed in Remark 2.2, when d = p and θ can be written as a linear
functional of f it is possible to use the different prior scheme described in Remark 2.2.
With this different prior scheme, the posterior of θ is recovered from the GP posterior of f
through the transformation b(f ).

2.3.2

Properties of the posterior distribution of θ

In this section we show two important results. The first one establishes that expression
(2.10) is invariant to the choice of the measure Π used to define the nuisance parameter
f and therefore the marginal posterior of θ is invariant to the choice of Π as well. More
precisely the following proposition holds.
Proposition 2.1. For a positive measure Π1 on S, let EΠ1 := L2 (S, BS , Π1 ) and z :=
dΠ
dΠ1 .

Let ϕ : E → EΠ1 be the transformation ϕ(f ) = f z and Φ be the set of measurable

transformations defined as
Φ :=



dΠ1 (x)
ϕ : E → EΠ1 ; ϕ(f ) = f z, Π1 is a positive measure and sup
<∞ .
x∈S dΠ(x)

Then, the marginal posterior distribution µ(θ|rn ) of θ is Φ-invariant.
This result shows that, once we integrate out the nuisance parameter f , the posterior
distribution of θ is not affected by the choice of the dominating measure Π which only
causes a transformation of the nuisance parameter.
Therefore, we can use two different dominating measures Π: one for the definition of
the nuisance parameter f and one for the definition of E, K and K ∗ . In particular, if
supx∈S

dF∗ (x)
dΠ(x)

< ∞ (where Π is used to define f and to construct its prior, and F∗ is used

to define E, K and K ∗ ) then, once we have specified the prior for the nuisance parameter

f = dF/dΠ, we deduce from it the prior of the transformation ϕ(f ) := f dΠ/dF∗ (see

13

the proof of Proposition 2.1). Therefore, the prior mean of ϕ(f ) is f0θ dΠ/dF∗ and the

prior covariance operator of ϕ(f ) writes in terms of an o.n.b. of L2 (S, BS , F∗ ) that we
still denote by {ϕj }j≥0 and where {ϕj }dj=1 are equal to the moment functions hj (θ, x),
P
j = 1, . . . , d, orthonormalized with respect to F∗ : Ω0θ · = j>d λj hϕj , ·iϕj . Moreover, the

operator K and its adjoint are defined by using F∗ instead of Π so that they are operators

from (resp. to) the space E = L2 (S, BS , F∗ ), and this does not change our inference on θ.

Therefore, the operator H becomes Hϕ = K|D∗ ϕ, ∀ϕ ∈ D∗ , and H −1 · = K −1 . Because

H −1 KΩ0θ K ∗ (H ∗ )−1 = Ω0θ |D∗ ∩R(K ∗ ) it follows that the eigenfunctions {ψj (θ)}j≥0 are given

by {ϕj }j≥1 ∈ D∗ ∩ L2 (S, BS , F∗ ) (which depend on θ and that are orthonormalized with

respect to F∗ ) with corresponding eigenvalues {ljθ }j≥0 = {λj 1{j > d}}j≥0 (which do not

depend on θ). Remark that if K is injective, then R(K ∗ ) = R(K ∗ ) since N(K) = {0} and

K ∗ is a closed operator, and R(K ∗ ) is dense in E. Therefore, we can simplify the expression

for µ(θ|rn ) to:

µ(θ|rn ) =

R

n
o
√
P
h n(rn −Kϕ(f0θ )),(H ∗ )−1 ψj (θ)i2
exp − 12 ∞
µ(θ)
j=0
1+nλj
.
o
n
√
1 P∞ h n(rn −Kϕ(f0θ )),(H ∗ )−1 ψj (θ)i2
µ(θ)dθ
exp − 2 j=0
1+nλj

(2.13)

Θ

This simplification and the definition of K, K ∗ and E in terms of Π = F∗ must be under-

stood in the following of the paper every time we explicitly assume supx∈S

dF∗ (x)
dΠ(x)

< ∞. It

is also important to remark that the orthonormalization of 1 and the moment functions
hj (θ, x), j = 1, . . . , d, with respect to F∗ can be implemented either by using the Cholesky

decomposition of (the empirical counterpart of) E∗ [(1, h(θ, X)T )T (1, h(θ, X)T )] or by using

the Gram-Schmidt process.
The second result we are going to show2 establishes a link between our Bayesian procedure, GEL estimators with quadratic criterion and the continuous updating GMM estimator. This relationship, given in Theorem 2.2 below, holds when the GP prior for f |θ is
allowed to become diffuse. More precisely, let us rescale the prior covariance operator of
f |θ by a positive scalar c so that the prior of f |θ may be written, for f0θ ∈ EM , as
µ(f |θ, c) ∼ GP(f0θ , cΩ0θ ),

Z

h(θ, x)f0θ (x)Π(dx) = 0,

dF∗ (x)
dΠ(x) < ∞, hj (θ, x) ∈
that E∗ [h(θ, Xi )h(θ, Xi )T ]

Theorem 2.2. Assume that supx∈S
∀j = 1, . . . , d, l > d and ∀θ ∈ Θ, and

1/2

Ω0θ (1, h(θ, ·)T )T = 0,

c ∈ R+ .

R(K ∗ ) and ϕl ∈ R(K ∗ ),

is nonsingular ∀θ ∈ Θ. Let

µ(f |θ, c) ∼ GP(f0θ , cΩ0θ ), with f0θ and Ω0θ satisfying Restrictions 1 and 2, and c ∈ R+ .
Let µ(θ|rn , c) denote the (marginal) posterior of θ obtained by integrating out f from P f
2

We thank Yuichi Kitamura for having suggested this research question.

14

with respect to µ(f |θ, c). Then,

 1
lim µ(θ|rn , c) ∝ exp −
c→∞
 2

where Vn (θ) :=

1
n

!T
n
1 X
√
h(θ, xi )
Vn (θ)−1
n
i=1

Pn

T
i=1 h(θ, xi )h(θ, xi ) .

!
n

X
1
√
h(θ, xi )
µ(θ)

n
i=1

Remarks that in the theorem the limit c → ∞ is taken after f has been marginalized

out. The result in the theorem deserves some comments. First, it shows that, as the
(conditional) prior on f becomes more and more diffuse, our marginal likelihood becomes
the quasi-likelihood function (also called LIL) that has been used often in the literature, for
instance by Chernozhukov and Hong [2003] and Kim [2002]. Therefore, the LIL naturally
arises from a nonparametric Bayesian procedure, which places a GP prior on the set of

functions in E that satisfy the EE, as the nonparametric prior becomes noninformative.

Second, Theorem 2.2 shows that, as the prior on f becomes noninformative, the MAP

objective function is the same (up to constants and up to the prior) as the GEL objective
function with quadratic criterion, see the proof of Theorem 2.1 in Newey and Smith [2004].
Moreover, as it can be deduced from Newey and Smith [2004, Theorem 2.1], the MAP
objective function becomes a monotonic transformation of the continuous updating GMM
objective function.

2.3.3

Properties of the MAP estimator of θ

By Proposition 2.1, our inference procedure is invariant to the choice of Π. Let us
dF∗ (x)
dΠ(x)

< ∞ and recall the discussion below Proposition 2.1 so that E =
R R
By Mercer’s formula, under the assumption T S |k(t, x)|2 F∗ (dx)ρ(dt) < ∞,
R
P
we can write k(t, x) = ∞
j=1 λjK φj1 (t)φj2 (x) where (λjK , φj1 , φj2 ) are such that S k(t, x)φj2 (x)F∗ (dx) =
R
λjK φj1 (t) and T k(t, x)φj1 (t)ρ(dt) = λjK φj2 (x), λjK > 0 and {φj1 }, {φj2 } are orthonorP
mal sequences of F and E, respectively. Therefore, K· = ∞
j=1 λjK φj1 (t)hφj2 , ·i and since

P∞ −1
∗
−1
(K ) ψj (θ) (t) = j=1 λjK φj1 (t)hφj2 , ϕj i under the assumptions of Theorem 2.1, we

assume supx∈S
L2 (S, BS , F∗ ).

obtain: ∀l ≥ 0

√
√
h nrn , (H ∗ )−1 ψl (θ)i = h nrn , (K ∗ )−1 ϕl i
n

∞

∞

1 XX X
1
hφj ′ 2 , ϕl i
=√
λjK hφj1 , φj ′ 1 iφj2 (xi )
λ
n
j′K
′
i=1 j=1 j =1

n

∞

1 XX
1 X
ϕl (xi ).
=√
φj2 (xi )hφj ′ 2 , ϕl i = √
n
n
i=1 j=1

15

i=1

It follows that the likelihood in (2.10) can be simplified and the MAP writes as:
θn = arg max µ(θ|rn ) = arg max (log pnθ (rn ; θ) + log µ(θ))
θ∈Θ

θ∈Θ


= arg max −
θ∈Θ

d
X
j=1

−
= arg min
θ∈Θ

d
X
j=1

+

√
h n(rn − Kϕ(f0θ )), (H ∗ )−1 ϕj i2

X√
h n(rn − Kϕ(f0θ )), (H ∗ )−1 ϕj i2
j>d

!2
n
1X
ϕj (xi ) − hϕ(f0θ ), ϕj i
n

1
1
+ 2 log µ(θ)
1 + nλj
2

i=1

X
j>d

!2
n

1
1X
ϕj (xi ) − hϕ(f0θ ), ϕj i
− log µ(θ)
n
1 + nλj

(2.14)

i=1

where we have eliminated the terms that do not depend on θ and we have used the fact
that hKϕ(f0θ ), (H ∗ )−1 ϕj i = hϕ(f0θ ), ϕj i. Equation (2.14) is quite useful and allows to

emphasize several aspects of our methodology.

I. The first term in (2.14) accounts for the EE. To see this, suppose that (1, h(θ, x)T )T
have been orthonormalized by using the Gram-Schmidt process: ϕ0 (x) = 1 and ϕj (x) =
Pj−1 ∗
cj (hj (θ, x) − l=1
E [hl (θ, X)hj (θ, X)]hl (θ, x) − E∗ [hj (θ, X)]), ∀1 ≤ j ≤ d where cj is

the normalization constant. Then, hϕ(f0θ ), ϕj i = −cj E∗ [hj (θ, X)], ∀1 ≤ j ≤ d and so
Pj−1 ∗
1 Pn
i=1 ϕj (xi ) − hϕ(f0θ ), ϕj i = cj (hj (θ, x) −
l=1 E [hl (θ, X)hj (θ, X)]hl (θ, x)) which corn
responds to the orthonormalized moment functions with respect to F∗ . Minimization of

this term corresponds to the classical GMM.

II. The second term in (2.14) accounts for the extra information that we have, namely,
the information contained in the subspace of E orthogonal to span{1, h1 (θ, ·), . . . , hd (θ, ·)}.

This information, which is in general not exploited in EE frameworks (in frequentist as well

as in Bayesian approaches), can be exploited thanks to the prior distribution and the prior
mean f0θ if the prior is not fixed but varies with n at an appropriate rate (see comment III
below). On the contrary, if the prior is fixed then, as n → ∞, the second term of (2.14)
P
converges to 0 since (1 + nλj )−1 → 0. Remark also that n−1 ni=1 ϕj (xi ) → E∗ [ϕj (X)] a.s.

and E∗ [ϕj (X)] = 0 because ϕj is orthogonal to 1 for j > d.

III. Expression (2.14) makes an explicit connection between the parametric case (infinite
number of moment restrictions) and the semiparametric case (where only the first d moment
conditions hold). The semiparametric case corresponds to the classical GEL or GMM approach while the parametric case corresponds to the maximum likelihood estimator (MLE).
Indeed, the prior distribution for f specifies a parametric model for f0θ which satisfies the d

16

moment restrictions and eventually other “extra” moment restrictions. The eigenvalues λj
of the prior covariance operator play the role of weights of the “extra” moment restrictions
and represent our “beliefs” concerning these restrictions. When we are very confident about
these “extra” conditions, or equivalently we believe that f0θ is close to f∗ , then the λj s are
close to zero or converge to 0 faster than n−1 as n → ∞. So, the prior distribution for f

is degenerate on f0θ (as n increases) when the parametric model is the true one. In that
case, the MAP estimator will essentially be equivalent to the MLE that we would obtain
if we use the prior mean function f0θ as the likelihood. When we are very uncertain about
f0θ then the λj s are very large and may tend to +∞ (uninformative prior). In this case
the MAP estimator will be close to the GMM estimator (up to a prior on θ).

2.3.4

Testing and moment selection procedures

Remark III in section 2.3.3 is important if one is interested in constructing testing
procedures or doing moment selection. We are not going to develop a formal test/selection
procedure here as this will make the object of a separated paper, but we would like to
point out that our procedure suggests an easy way to test a parametric model against a
semiparametric one characterized by a finite number of moment restrictions. We can deal
with the two following situations where we assume supx∈S dF∗ (x)/dΠ(x) < ∞:
[1. ] We know that the distribution of the data F∗ satisfies d moment restrictions and

we want to test whether it has a particular parametric form. In this case, for a given
R
probability distribution G such that h(θ, x)G(dx) = 0 for a known vector of functions
h(θ, x), the null hypothesis is H0 : F∗ = G. An example is the univariate linear regression

model: Y = Zθ + ε, where f∗ is the true joint pdf of X := (Y, Z)T and E∗ (Y |Z) = Zθ so
that h(θ, X) = Y − Zθ. We may want to test that f∗ belongs to a particular parametric

class.

[2. ] There are d EE of which we are sure and we want to test the validity of the other
D − d EE. The null hypothesis writes H0 : E∗ (hj (θ, X)) = 0 for d < j ≤ D.
To treat [1.], we have to specify f0θ = dG/dΠ. To treat [2.] we have to specify f0θ such
R
that it satisfies the extra EE: S hj (θ, x)f0θ (x)Π(dx) = 0, d < j ≤ D, and the {ϕj }D
j=d+1

as the orthonormalized hj (θ, x), d < j ≤ D, with respect to F∗ (so that the square term in

the second term of (2.14) is zero only at the true θ∗ if the extra EE are correctly specified).

Then, for both the situations, the natural approach would be to treat the λj s corresponding

to the extra EE as hyperparameters for which a prior distribution is specified. The null
hypothesis writes H0 : λj = 0 for all j > d in [1.], and H0 : λj = 0 for all d < j ≤ D in [2.].

Then, the posterior distribution of λj may be used to draw a conclusion on the test: either
by considering posterior odds ratio or by constructing encompassing tests.

17

P∞

To construct a prior for the λj s let us write: λj = cρj , j > d, where c = trΩ0 and

j=0 ρj

= 1. The {λj }dj=1 corresponding to the first d EE do not affect the trace of Ω0

since they are equal to 0. We propose two alternatives priors.

Dirichlet prior. Suppose that we want to test H0 : λj = 0 for all d < j ≤ D. Then one
may specify a Dirichlet prior for (ρd+1 , . . . , ρD−1 ):

µρ (ρd+1 , . . . , ρD−1 |ν) ∝

D−1
Y

j=d+1



ν −1 

ρj j

1−

D−1
X

j=d+1

νD −1

ρj 

D−1
Y

j=d+1



I(ρj ≥ 0)I 

D−1
X
j=1



ρj ≤ 1

where ν = (νd+1 , . . . , νD ).
Prior on c > 0. Suppose that we want to test H0 : λj = 0 for all j > d. Thus, the null

hypothesis may be written H0 : c = 0 and a prior for c may be any distribution with support
contained in the positive real semi-axis, for example an inverse gamma distribution.

3

Asymptotic Analysis
In this section we focus on the frequentist asymptotic properties of our approach for

n → ∞. For this analysis we use the true probability measure P ∗ which corresponds to

the true DGP F∗ and maintain the assumption that supx∈S

dF∗ (x)
dΠ(x)

< ∞ where Π is used to

define the nuisance parameter f and F∗ to define K, K ∗ and E (see the discussion below
Proposition 2.1). Therefore, ϕ(f ) = f∗−1 f . We show: (i) consistency of the posterior of θ

(Theorem 3.1), (ii) frequentist consistency of the MAP estimator θn (Theorem 3.2), and
(iii) convergence in Total Variation of µ(θ|rn ) towards a normal distribution (section 3.2).
In the following, for every θ̃ ∈ Θ and δ > 0 we denote by B(θ̃, δ) the closed ball centered

in θ̃ with radius δ, that is, B(θ̃, δ) = {θ ∈ Θ; kθ − θ̃k ≤ δ}, where here k · k denotes the
Euclidean norm in Rp . Moreover, denote δn = n−1/2 .

3.1

Posterior Consistency

We first state the following assumptions.
A1. The true parameter θ∗ belongs to the interior of a compact convex subset Θ of Rd and
is the unique solution of E∗ [h(θ, X)] = 0.

A2. The singular functions {ψj (θ), ρj (θ)} are continuous functions of θ.
A3. The prior mean function f0θ and the prior distribution µ(θ) are continuous in θ.
√
A4. For every j and every n, E∗ h n(rn − Kϕ(f∗ )), (H ∗ )−1 ψj (θ)i4 < ∞.

18

Assumption A1 is a standard assumption in the literature on EE. Assumptions A2 and
A3 can be easily satisfied since f0θ and the operators Ω0θ , K and Σ are chosen by the
statistician. The next theorem gives concentration of the posterior distribution around the
true value θ∗ .
Theorem 3.1. Let the assumptions of Theorem 2.1 and Assumptions A1-A4 be satisfied
and assume supx∈S

dF∗ (x)
dΠ(x)

< ∞. Then, for any prior µ(θ) that puts enough mass in a

neighborhood of θ∗ and any sequence Mn → ∞,
µ

√

in P ∗ -probability as n → ∞.


nkθ − θ∗ k > Mn |rn → 0

(3.1)

Given the result of the theorem, the next result follows.
p

Theorem 3.2. Let the assumptions of Theorem 3.1 be satisfied. Then, θn −
→ θ∗ in P ∗ probability as n → ∞.

3.2

Asymptotic Normality

In this section we establish asymptotic normality of µ(θ|rn ) for the Bayesian model
described in (2.7). This result applies to the case d > p (which is our main interest) as well
as to the case d = p. In appendix ?? we establish, under different assumptions, asymptotic
normality of µ(θ|rn ) for the case d = p described in Remarks 2.2 and 2.3 where the prior
for θ is deduced from the prior for f .
For some τ ∈ Rp , let
sn (τ ) := pn,θ∗ +δn τ (rn ; θ∗ + δn τ ).
We assume that there exist a random vector ℓ̃∗ and a nonsingular matrix I˜∗ (that depend
on the true θ∗ and f∗ ) such that the sequence ℓ̃∗ is bounded in probability, and satisfy
log

1
1
sn (τ )
= √ τ T Ie∗ ℓe∗ − τ T Ie∗ τ + op (1)
sn (0)
n
2

(3.2)

for every random sequence τ which is bounded in P ∗ -probability. Condition (3.2) is known
as the integral local asymptotic normality assumption which is used to prove asymptotic
normality of semiparametric Bayes procedures, see e.g. Bickel and Kleijn [2012]. In Appendix ?? we prove that, if supx∈S

dF∗ (x)
dΠ(x)

< ∞, then equation (3.2) holds with





 ∗

T −1 ∗ ∂h(θ∗ , X)
∗ ∂h(θ∗ , x)
e
E
E h(θ∗ , X)h(θ∗ , X)
I∗ = −E
∂θ
∂θ T

19



if E∗ h(θ∗ , X)h(θ∗ , X)T is nonsingular. For two probability measures P1 and P2 absolutely

continuous with respect to a positive measure Q, define the total variation (TV) distance
R
as ||P1 − P2 ||T V := |f1 − f2 |dQ where f1 and f2 are the Radon-Nikodym derivatives of P1

and P2 , respectively, with respect to Q. The following theorem shows that under (3.2) the
√
posterior distribution of n(θ − θ∗ ) converges in the TV distance to a Normal distribution
with mean ∆∗ := √1 ℓe∗ and variance Ie−1 .
∗

n

Theorem 3.3. Let the assumptions of Theorem 2.1 be satisfied. Assume that A1-A3,

(3.1) and (3.2) hold and that the prior µ(θ) puts enough mass in a neighborhood of θ∗ . If
√
√
µ( n(θ − θ∗ )|rn ) denotes the posterior of n(θ − θ∗ ), then:
√
kµ( n(θ − θ∗ )|rn ) − N (∆∗ , Ie∗−1 )kT V → 0

(3.3)

in P ∗ -probability as n → ∞.

4

Implementation
In this section we first explain the numerical implementation of our procedure and then

show the results of three simulations. Let M be the number of discretization points and
xeval := (x1 , . . . , xM ) (resp. teval := (t1 , . . . , tM )) be (equidistant) discretization points for
functions in E (resp. in F). The procedure consists of Markov Chain Monte Carlo (MCMC)

draws from the posterior (2.11) based on the Metropolis-Hastings (M-H) algorithm, see for
instance Robert [2002]. Beforehand, the researcher must chose: (a) the measures Π(x)
and ρ(t), (b) the function k(t, x) (for instance: k(t, x) = 1{x ≤ t}, k(x, t) = eitx , or
k(x, t) = etx ), (c) the sequence λj (for instance: λj = j −α with α > 1, or α−j with α > 1),

(d) the prior µ(θ) and (e) the auxiliary distribution, say g(θ) for the M-H. To construct Ω0θ
the series is truncated at a value J which is chosen such that λJ is small (for instance of
the order 10−3 ). In the following, functions and operators are replaced by their discretized
versions, denoted with a subindex d.
Before starting the M-H algorithm, compute the following quantities:
I. Compute rn and discretize it to get the M -vector rn,d := (rn (ti ))i=1,...,M .
II. Compute Kd as the M × M matrix with generic element (k(ti , xj )Π(xj ))i,j=1,...,M .
III. Compute Σd as the M × M matrix with generic element (σ(ti , tl )ρ(tl ))i,l=1,...,M where
P
σ(ti , tl ) = n−1 nj=1 (k(ti , xj )k(tl , xj )) − rn (ti )rn (tl ). Due to the discretization, the
matrix Σd could be ill-conditioned in some cases and then has to be regularized.

Then, implement the M-H algorithm where for every value θ̃ (j) drawn from the auxiliary
distribution g(θ), the expression of log pnθ (rn ; θ) in (2.10) is evaluated through the following

20

steps:
1. Evaluate the vector function h(θ, x) at θ̃ (j) and the M discretization points xeval so
to have a (d × M ) matrix hd := (h(θ̃ (j) , x1 ), . . . , h(θ̃ (j) , xM )). Let ιM denotes the M

dimensional row-vector of ones.

2. Orthonormalize the set of vectors given by ιM and the rows of hd (by using e.g. the
Gram-Schmidt process) to obtain the d + 1 M -vectors ϕj,d := (ϕj (x1 ), . . . , ϕj (xM ))T for
0 ≤ j ≤ d.
3. Complete the basis by first constructing J − d − 1 vectors ιM and then orthonormalizing

them sequentially by using e.g. the Gram-Schmidt process to obtain the remaining

components of the basis ϕj,d , j = d + 1, . . . , J. Denote ϕM := (ϕj (xi ))

i=1,...,M
j=0,...,J−1

ϕΠ := (ϕj (xi )Π(xi )step)

i=1,...,M ,
j=0,...,J−1

and

which have both dimension J × M and where step

denotes the discretization step.
4. Construct Ω0θ,d as Ω0θ,d = ϕTM diag(λj )j=0,...,J ϕΠ , where λj = 0 for j = 0, . . . , d.
5. Select a prior mean function f0θ and discretize it to get an M -vector f0θ,d := (f0θ (xi ))i=1,...,M .
6. Compute the singular value decomposition (ljθ , ρj (xi ), ψj (ti ))j=0,...,J−1 of the matrix
i=1,...,M
√
p
√
( Σd )−1 Kd Ω0θ,d , where for a positive definite matrix A, A denotes its positive square
root. Denote ψj,d := (ψj (t1 ), . . . , ψj (tM ))T .

7. Compute the terms in log pnθ (rn ; θ) that depend on θ by discretizing the scalar product:
2
√
J
J
 1X
(rn,d − Kd f0θ,d )T diag(ρ(ti ))i=1,...,M ( Σd )−1 ψj,d
1X
2
log 1 + nlj +
.
2
2
2
1 + nljθ
j=1

j=1

We remark that the Gram-Schmidt orthonormalization in steps 2 and 3 is extremely simple
and fast to implement even for complicated moment functions. By using the decomposition
Σ = HH ∗ with H, H ∗ as described in (2.6), the previous steps can be simplified in the
√
following way: steps III and 6 do not have to be implemented and ( Σd )−1 in step 7 has
to be replaced by (K ∗ )−1 . Moreover, ψj,d and lj2 have to be replaced by (ϕj (xi ))i=1,...,M
and λj , respectively. We now show the finite sample properties of our GP approach with

the help of three simulations. The interest in using a GP prior will be made evident in the
more complicated examples where there are overidentifying restrictions.

21

4.1

Just identification and prior on θ through µ(f )

Let the parameter θ of interest be the population mean: θ =

R

xf (x)dx and h(θ, x) =

(θ − x). This example considers the d = p case where θ is a linear functional of f and

the prior of θ is deduced from the prior of f , denoted by µ(f ). Suppose that the support
of F∗ is R, so that θ can take every value in R. Then the prior µ(f ) is a GP which is

unrestricted except for the fact that it must generate trajectories that integrate to 1 a.s.,
1/2

namely, µ(f ) ∼ GP (f0 , Ω0 ) where f0 is a pdf and Ω0 is such that Ω0 1 = 0. Therefore, the
prior distribution of θ is Gaussian with mean hf0 , ιi and variance hΩ0 ι, ιi where ι denotes

the identity functional, that is, ι(x) = x. The posterior distribution of θ is as in (??) with
−1
g = ι and A = Ω0 K ∗ n−1 Σ + KΩ0 K ∗
. We illustrate how to construct in practice the

covariance operator Ω0 . Let S = R; the Hermite polynomials {Hj }j≥0 form an orthogonal
basis of L2 (R, B, Π) for dΠ(x) = e−x

2 /2

dx and can be used to construct the eigenfunctions

of Ω0 . The first few Hermite polynomials are {1, x, x2 − 1, (x3 − 3x), . . .} and they are
√
R
2
orthogonal with respect to Π: R Hl (x)Hj (x)e−x /2 dx = 2πn!δlj , where δlj is equal to 1

if l = j and to 0 otherwise. The operator Ω0 is constructed as
Ω 0 · = σ0

∞
X
j=0

λj √

1
hHj , ·iHj
2πn!

where Hj+1 (x) = xHj (x) − jHj−1 (x), λ0 = 0 and {λj , j ≥ 1} = {aj , j ≥ 1} with a < 1.

In our simulation exercise we generate n = 1000 i.i.d. observations (x1 , . . . , xn ) from a
P
N (1, 1) distribution and construct the function rn = n−1 ni=1 etxi as the empirical Laplace
transform. Therefore, f∗ (x) =

√1 e−(1−2x)/2
2π

the operators K and K ∗ take the form
∀φ ∈ E, Kφ =

Z

R

tx

−x2 /2

e φ(x)e

dx

and θ∗ = 1. We set T = R and ρ = Π. Thus,

and ∀ψ ∈ F,

∗

K ψ=

Z

2 /2

etx ψ(t)e−t

dt.

R

The prior mean function f0 is set equal to the Π-pdf of a N (̺, 1) distribution: f0 (x) =
2
√1 e−(̺ −2̺x)/2 ,
2π

̺ = 2. We show in Figure 1 the prior and posterior distribution of θ. We

also show the prior mean (magenta asterisk), the posterior mean (blue asterisk) and the
MAP (red asterisk) of θ. The posterior mean of θ is computed by discretizing the inner
product hE(f |rn ), ιi. The pictures are obtained for a = 0.3 and σ0 = 1. The number of

discretization points, used to approximate the integrals, is equal to 1000.

4.2

Just identification and prior on θ

We consider the same framework as in the previous example where the parameter θ of
R
interest is the population mean: θ = xf (x)dx and h(θ, x) = (θ − x), but now we specify

22

Prior distributions of θ

0.4

0.7

0.3

0.6

0.25

0.5

µθ

µθ

0.35

Posterior distributions of θ

0.8
prior distribution
prior mean of θ = 2

0.2
0.15

0.3

0.1

0.2

0.05

0.1

0

posterior distribution
posterior mean of θ=0.9657
MAP=0.9656

0.4

0
-2

-1

0

1

2

3

4

5

6

-1

-0.5

θ

0

0.5

1

1.5

2

2.5

3

θ

Figure 1: Prior and posterior distributions and means of θ. The true value of θ is θ∗ = 1.
a joint proper prior distribution on (θ, f ). We specify a marginal prior µ(θ) on θ and a
conditional prior on f given θ. While µ(θ) can be arbitrarily chosen, µ(f |θ) is specified as

a GP constrained to generate functions that integrate to 1 and that have mean equal to θ

a.s., as described in section 2.1.

Compared to the approach in section 4.1, this approach allows to easily incorporate any
prior information that one may have about θ. In fact, incorporating the information on θ
through the prior distribution on f is complicated while to incorporate such an information
directly in the prior distribution of θ results to be very simple. In particular, the approach
of this section works even when θ takes values in a compact subset of Rp , while the approach
of section 4.1 does not work in this case.
Let us suppose that m = 1, S = [−1, 1] and let Π and ρ be the Lebesgue measure. Then,
the covariance operator Ω0θ can be constructed by using Legendre polynomials since the
second Legendre polynomial P1 (x) = x allows to implement the constraint on θ. Because
the moment function is separable in θ and x, the prior covariance operator does not depend
on θ, so that we denote it by Ω0 . The first few Legendre polynomials are {1, x, (3x2 −

1)/2, (5x3 − 3x)/2, . . .} and they are orthogonal with respect to the L2 inner product on
R1
[−1, 1]: −1 Pl (x)Pj (x)dx = 2/(2j + 1)δlj , where δlj is equal to 1 if l = j and to 0 otherwise.

Moreover, the Legendre polynomial obey the recurrence relation (j + 1)Pj+1 (x) = (2j +
1)xPj (x)−jPj−1 (x) which is useful for computing Ω0 in practice. The normalized Legendre
polynomials form a basis for L2 [−1, 1] so that we can construct the operator Ω0 as
Ω 0 · = σ0

∞
X
j=2

λj

2j + 1
hPj , ·iPj
2

where we have set λ0 = λ1 = 0 in order to implement the constraints. The remaining
P
λj , j ≥ 2 can be chosen in an arbitrary way provided that j≥2 λj < ∞. The constant

σ0 > 0 can be set to an arbitrary value and has the purpose of tuning the size of the prior
covariance.

23

Many orthogonal polynomials are suitable for the construction of Ω0θ and they may be
used to treat cases where S is different from [−1, 1].
We perform two simulations exercises: the first one makes use of the empirical cdf to
P
construct rn : rn (t) = Fn (t) := n−1 ni=1 1{xi ≤ t} and the second one uses the empirical
P
moment generating function rn (t) = n−1 ni=1 etxi . In both the simulations we use Legendre
polynomials and we generate n = 1000 i.i.d. observations (x1 , . . . , xn ) from a N (0, 1)

distribution truncated to the interval [−1, 1]. The prior distribution for θ is uniform over
the interval [−1, 1]. The prior mean function f0θ is taken equal to the pdf of a Beta
distribution with parameters pθ and q and with support [−1, 1]:
f0θ (x) =

(x + 1)pθ −1 (1 − x)q−1
.
B(pθ , q)2pθ +q−1

(4.1)

We use the notation pθ to stress the dependence on θ of this shape parameter. We fix q = 2
R1
and recover pθ such that −1 xf0θ (x)dx = θ. It is easy to see that for our Beta distribution:
R1
pθ −q
−1 xf0θ (x)dx = pθ +q . The covariance operator Ω0θ is constructed by using the Legendre
polynomials, λj = j −1.7 and σ0 = 5.

Since the posterior distribution µ(θ|rn ) can not be computed in a closed-form we simulate from it by using a M-H algorithm. To implement this algorithm we have to select an
auxiliary pdf g(θ). We summarize the simulation scheme.
1. Draw a n i.i.d. sample (x1 , . . . , xn ) from f∗ (where f∗ is a N (0, 1) truncated to [−1, 1]);
2. compute rn = Fn or rn = n−1

Pn

txi
i=1 e ;

3. draw θ ∼ U [−1, 1] and denote it θ̃;
4. compute pθ as pθ =

(θ̃+1)2
1−θ̃

(where we have fixed q = 2);

5. compute f0θ as in (4.1) with parameters (pθ , q = 2) and discretize it;
6. draw θ from the marginal posterior distribution of θ by using a M-H algorithm with the
following auxiliary pdf (triangular distribution):

g(ξ; θ) =

1−ξ
ξ+1
1{ξ ∈ [−1, θ)} +
1{ξ ∈ [θ, 1]}
θ+1
1−θ

as described above. We draw 10000 values and discard the first 5000. The initial value
for the algorithm is θ = 0.5.
We represent in Figure 2a the results for the simulation with rn (t) = Fn (t) and in
P
Figure 2b the results for rn = n−1 ni=1 etxi : the blue asterisk represents the posterior

24

mean estimate while the red asterisk represents the MAP estimate. These figures also show
the marginal posterior distribution of θ (dashed blue line) approximated by using a kernel
smoothing and 5000 drawings from the posterior. In both the simulations, n = 1000 and
the number of discretization points, used to approximate the integrals, is equal to 1000.
Posterior density function of θ

Posterior density function of θ
0.8

0.7

0.7

0.6

0.6

0.5

0.5
µ (θ|F

0.4

θ

0.4

θ

µ (θ|F

hat

hat

)

)

0.8

0.3

0.3
posterior distribution
posterior mean of θ = 0.0313
MAP = 0.0045

0.2
0.1
0
−1

−0.5

0
θ

0.5

0.2

posterior distribution
posterior mean of θ = −0.0092
MAP = −0.0207

0.1
0
−1

1

(a) rn (t) = Fn (t)

−0.5

(b) rn = n

0
θ

0.5

Pn
−1

i=1

1

etxi

Figure 2: Estimations of θ based on the posterior distribution: posterior mean and MAP.
The true value of θ is θ∗ = 0 and n = 1000.

4.3

Overidentified case

Let us consider the case in which x is univariate and the one-dimensional parameter
of interest θ is characterized by the moment conditions EF (h(θ, X)) = 0 with h(θ, x) =
(x − θ, 2θ 2 − x2 )T . For instance, this arises when the true DGP F∗ is an exponential distribution with parameter θ. The prior µ(θ) is specified as a U [θ∗ − 1, θ∗ + 1].

The moment conditions are incorporated in the prior µ(f |θ) for f as described in sec-

tion 2.1. We chose Π(dx) = e−x dx and rn (t) = Fn (t). We first orthonormalize the moment
functions 1, x−θ, 2θ 2 −x2 with respect to Π and then complete the bases by using the Gram-

Schmidt orthonormalization process. The inner products in E are approximating by using

the trapezoidal rule on equally spaced subintervals of the interval [min xi − 1, max xi + 1].

We use polynomially decreasing eigenvalues for Ω0θ : λj = j −1.7 . Finally, to construct Ω0θ

we truncate the series at J = 300 since after that the value of λj is of the order 10−5 and
then can be considered zero.
In our simulation, we generate n = 500 observations x1 , . . . , xn from an exponential
distribution with parameter θ∗ = 2. The measure ρ(dt), necessary to construct K ∗ , is taken

equal to the Lebesgue measure. Operators K and K ∗ are approximated by using the trapezoidal rule on equally spaced subintervals of the following intervals: [min xi − 1, max xi + 1]

for K and [min xi , max xi ] for K ∗ . The operator Σ is approximated in a similar way. Because of this discretization, the operator Σ is ill-conditioned and hence we regularize it by
adding to it the identity matrix scaled by n−1 .

25

The prior mean function f0θ is chosen by using a two-step procedure where in the first
step we compute f˜ = (0.1I + K ∗ K)−1 K ∗ rn and in the second step we project it on Λ(θ)
for a given θ. Alternatively, f0θ can be constructed a described before Remark 2.2.
To draw from the posterior distribution of θ, we use a M-H algorithm. To implement
this algorithm we use, as auxiliary distribution, a χ2⌈θ⌉ distribution. The posterior distri-

bution, its mean and its mode obtained in this simulation are plotted in Figure 3a. The
posterior density function has been obtained by kernel smoothing with a Gaussian kernel
and a bandwidth equal to 0.3.
Finally, we have repeated the same Monte Carlo simulation 100 times and have computed the average of the posterior mean estimators and MAP estimators. We report the
results in Figure 3b together with the posterior density, mean and MAP obtained in each
simulation.
Posterior density function of θ

Posterior density function of θ
0.8

0.5
0.45

0.7

0.4
0.6
0.35
0.5
µ(θ|F

µ(θ|F

hat

hat

)

)

0.3
0.25
0.2

0.4
0.3

0.15
posterior distribution
Mean posterior mean of θ = 1.9894
Mean MAP = 1.9315

0.2
posterior distribution
posterior mean of θ = 2.0186
MAP = 2.1305

0.1
0.05
0

1

1.5

2
θ

2.5

0.1
0

3

(a) Result for one simulation.

1

1.5

2
θ

2.5

3

(b) Average of posterior means and MAP
estimators over 100 Monte Carlo iterations.

Figure 3: Overidentified case d > p. Posterior distributions of θ, mean and MAP estimators. rn = Fn and the true value of θ is θ∗ = 2.

References
P. J. Bickel and B. J. K. Kleijn. The semiparametric bernstein–von mises theorem. The
Annals of Statistics, 40(1):206–237, 2012.
L. Bornn, N. Shephard, and R. Solgi. Moment conditions and bayesian nonparametrics.
Technical report, arXiv:1507.08645, 2015.
M. Broniatowski and A. Keziou. Divergences and duality for estimation and test under
moment condition models. Journal of Statistical Planning and Inference, 142(9):2554 –
2573, 2012.

26

G. Chamberlain and G. W. Imbens. Nonparametric applications of bayesian inference.
Journal of Business & Economic Statistics, 21(1):12–18, 2003.
S. Chaudhuri, D. Mondal, and T. Yin. Hamiltonian monte carlo sampling in bayesian
empirical likelihood computation. Journal of the Royal Statistical Society, Series B,
forthcoming, 2017.
S. X. Chen and I. Van Keilegom. A review on Empirical Likelihood methods for regression.
TEST, 18(3):415–447, 2009.
V. Chernozhukov and H. Hong. An mcmc approach to classical estimation. Journal of
Econometrics, 115(2):293–346, 2003.
S. Chib, M. Shin, and A. Simoni. Bayesian empirical likelihood estimation and comparison
of moment condition models. Technical report, arXiv:1606.02931, 2016.
J.-P. Florens and A. Simoni. Regularized posteriors in linear ill-posed inverse problems.
Scandinavian Journal of Statistics, 39(2):214–235, 2012.
J.-P. Florens and A. Simoni. Regularizing priors for linear inverse problems. Econometric
Theory, FirstView, 11 2014.
M. Grendar and G. Judge. Asymptotic equivalence of empirical likelihood and bayesian
map. The Annals of Statistics, 37(5A):2445–2457, 2009.
L. P. Hansen. Large sample properties of generalized method of moments estimators.
Econometrica, 50:1029–1054, 1982.
L. P. Hansen and K. J. Singleton. Generalized instrumental variables estimation of nonlinear
rational expectations models. Econometrica, 50(5):1269–1286, 1982.
L. P. Hansen, J. Heaton, and A. Yaron. Finite-sample properties of some alternative gmm
estimators. Journal of Business & Economic Statistics, 14(3):262–280, 1996.
G. W. Imbens. One-step estimators for over-identified generalized method of moments
models. The Review of Economic Studies.
J.-Y. Kim. Limited information likelihood and bayesian analysis. Journal of Econometrics,
107(1-2):175–193, 2002.
Y. Kitamura. Empirical likelihood methods with weakly dependent processes. Ann. Statist.,
25(5):2084–2102, 1997.

27

Y. Kitamura and T. Otsu. Bayesian analysis of moment condition models using nonparametric priors. Technical report, Department of Economics: Yale University, 2011.
Y. Kitamura and M. Stutzer. An information-theoretic alternative to generalized method
of moments estimation. Econometrica, 65(4):861–874, 1997.
Y. K. Kwan. Asymptotic bayesian analysis based on a limited information estimator.
Journal of Econometrics, 88(1):99–121, 1999.
N. A. Lazar. Bayesian empirical likelihood. Biometrika, 90(2):319–326, 2003.
Y. Liao and W. Jiang. Posterior consistency of nonparametric conditional moment restricted models. The Annals of Statistics, 39(6):pp. 3003–3031, 2011.
W. K. Newey and R. J. Smith. Higher order properties of gmm and generalized empirical
likelihood estimators. Econometrica, 72(1):219–255, 2004.
A. B. Owen.

Empirical likelihood ratio confidence intervals for a single functional.

Biometrika, 75(2):237–249, 1988.
A. B. Owen. Empirical Likelihood. Chapman & Hall/CRC Monographs on Statistics &
Applied Probability, 2001.
J. Qin and J. Lawless. Empirical likelihood and general estimating equations. The Annals
of Statistics, 22(1):300–325, 03 1994.
J. Rao and C. Wu. Bayesian pseudo-empirical-likelihood intervals for complex surveys.
Journal of the Royal Statistical Society Series B, 72(4):533–544, 2010.
C. Robert.

Méthodes de calcul en analyse bayésienne.

In J. F. J.J. Droesbeke and

G. Saporta, editors, Méthodes Bayésiennes en statistique, pages 149 – 181. TECHNIP,
2002.
S. M. Schennach. Bayesian exponentially tilted empirical likelihood. Biometrika, 92(1):
31–46, 2005.
S. M. Schennach. Point estimation with exponentially tilted empirical likelihood. The
Annals of Statistics, 35(2):634–672, 2007.
M. Shin. Bayesian GMM. Technical report, University of Pennsylvania, 2014.
R. J. Smith. Alternative semi-parametric likelihood approaches to generalised method of
moments estimation. The Economic Journal, 107(441):503–519, 1997.

28

