Maximum Approximated Likelihood Estimation

TL;DR

This paper extends the theoretical foundation of maximum approximated likelihood (MAL) estimators, demonstrating their consistency and asymptotic normality, and compares their performance with other approximation methods through examples and simulations.

Contribution

It generalizes existing results for maximum simulated likelihood to a broader class of MAL estimators, providing conditions for their statistical properties.

Findings

MAL estimators are consistent under general conditions.

MAL estimators are asymptotically normal.

Simulation results show MAL performs well compared to other methods.

Abstract

Empirical economic research frequently applies maximum likelihood estimation in cases where the likelihood function is analytically intractable. Most of the theoretical literature focuses on maximum simulated likelihood (MSL) estimators, while empirical and simulation analyzes often find that alternative approximation methods such as quasi-Monte Carlo simulation, Gaussian quadrature, and integration on sparse grids behave considerably better numerically. This paper generalizes the theoretical results widely known for MSL estimators to a general set of maximum approximated likelihood (MAL) estimators. We provide general conditions for both the model and the approximation approach to ensure consistency and asymptotic normality. We also show specific examples and finite-sample simulation results.

Tables

Table 1. Table 1: Total number of integrand evaluations for MALE with different integration schemes.

Setting	Error bound	Total cost $n\cdot R\left(n\right)$
Monte Carlo	$𝒪\left(r^{-1/2}\right)$	$𝒪\left(n^2\right)$
Quasi – Monte Carlo	$𝒪\left(r^{-1+\epsilon }\right)$	$𝒪\left(n^{3/2+\epsilon }\right)$
Gaussian quadrature, smoothness $k=2$	$𝒪\left(r^{-2}\right)$	$𝒪\left(n^{5/4}\right)$
Gaussian quadrature, analytic function	$𝒪\left(e^{-r}\right)$	$𝒪\left(n\log \left(n\right)\right)$
Sparse grid, bounded domain	$𝒪\left(r^{-k+\epsilon }\right)$	$𝒪\left(n^{1+\frac{1}{2k}+\epsilon }\right)$
Sparse grid, unbounded domain	$𝒪\left(r^{-k/2+\epsilon }\right)$	$𝒪\left(n^{1+\frac{1}{k}+\epsilon }\right)$

Equations

$$L_n(\mathbf{\theta })=\frac{1}{n}\sum _{i=1}^n\log f({\mathbf{z}}_i,\mathbf{\theta }).$$

$$f(\mathbf{z},\mathbf{\theta })={\int }_{\Omega }\phi (\mathbf{v},\mathbf{z},\mathbf{\theta })\omega (\mathbf{v})\mathrm{d}\mathbf{v},$$

$${\tilde{f}}_r(\mathbf{z},\mathbf{\theta }):=\sum _{j=1}^r{w}_{j,r}\phi ({\mathbf{v}}_{j,r},\mathbf{z},\mathbf{\theta })\approx f(\mathbf{z},\mathbf{\theta }).$$

$${\hat{\mathbf{\theta }}}_{MAL}^{(n)}:=\underset{\mathbf{\theta }\in \Theta }{\mathrm{argmax}}{\tilde{L}}_n(\mathbf{\theta }),\text{ where }{\tilde{L}}_n(\mathbf{\theta }):=\frac{1}{n}\sum _{i=1}^n\log {\tilde{f}}_{R(n)}({\mathbf{z}}_i,\mathbf{\theta }).$$

$${\hat{\mathbf{\theta }}}_M=\arg \underset{\mathbf{\theta }\in \Theta }{\max }{Q}_n(\mathbf{\theta }),$$

$${\hat{\mathbf{\theta }}}_{AM}=\arg \underset{\mathbf{\theta }\in \Theta }{\max }{\tilde{Q}}_n(\mathbf{\theta }).$$

$$\underset{n\to \infty }{\mathrm{plim}}{\hat{\mathbf{\theta }}}_M={\mathbf{\theta }}_0.$$

$$\underset{\mathbf{\theta }\in \Theta }{\sup }\left|{\tilde{Q}}_n(\mathbf{\theta })-Q_0(\mathbf{\theta })\right|$$

$$\le \underset{\mathbf{\theta }\in \Theta }{\sup }\left|{\tilde{Q}}_n(\mathbf{\theta })-Q_n(\mathbf{\theta })\right|+\underset{\mathbf{\theta }\in \Theta }{\sup }\left|Q_n(\mathbf{\theta })-Q_0(\mathbf{\theta })\right|.$$

$$\sqrt{n}\left({\hat{\mathbf{\theta }}}_{AM}-{\mathbf{\theta }}_0\right)\stackrel{d}{\to }N(0,H^{-1}\Sigma H^{-1}),(\text{as }n\to \infty )$$

$$\sqrt{n}{\nabla }_{\mathbf{\theta }}{\tilde{Q}}_n({\mathbf{\theta }}_0)$$

$$\displaystyle\sup_{{\boldsymbol{\theta}}\in\mathcal{N}}\big{\|}\nabla_{{\boldsymbol{\theta}}{\boldsymbol{\theta}}}\tilde{Q}_{n}({\boldsymbol{\theta}})-H({\boldsymbol{\theta}})\big{\|}$$

$$Q_n(\mathbf{\theta })=L_n(\mathbf{\theta })=\frac{1}{n}\sum _{i=1}^n\log f({\mathbf{z}}_i;\mathbf{\theta }),$$

$${\tilde{L}}_{n,r}(\mathbf{\theta }):=\frac{1}{n}\sum _{i=1}^n\log {\tilde{f}}_r({\mathbf{z}}_i;\mathbf{\theta })$$

$${\tilde{Q}}_n(\mathbf{\theta })={\tilde{L}}_n(\mathbf{\theta })=\frac{1}{n}\sum _{i=1}^n\log {\tilde{f}}_{R(n)}({\mathbf{z}}_i;\mathbf{\theta }).$$

$${\hat{\mathbf{\theta }}}_{MAL}\stackrel{p}{\to }{\mathbf{\theta }}_0.$$

$$|{\tilde{f}}_r(\mathbf{z};\mathbf{\theta })-f(\mathbf{z};\mathbf{\theta })|<\bar{\delta }-\delta .$$

$$\underset{\mathbf{\theta }\in \Theta }{\sup }\left|{\tilde{L}}_n(\mathbf{\theta })-L_n(\mathbf{\theta })\right|\stackrel{p}{\to }0.$$

$$\underset{\mathbf{\theta }\in \Theta }{\sup }\left|{\tilde{L}}_n(\mathbf{\theta })-L_n(\mathbf{\theta })\right|$$

$$\le \underset{\mathbf{\theta }\in \Theta ,\mathbf{z}\in \mathcal{Z}}{\sup }\left|\log ({\tilde{f}}_{R(n)}(\mathbf{z},\mathbf{\theta }))-\log (f(\mathbf{z},\mathbf{\theta }))\right|.$$

$$\underset{\mathbf{\theta }\in \Theta }{\sup }\left|{\tilde{L}}_n(\mathbf{\theta })-L_n(\mathbf{\theta })\right|$$

$$\mathcal{E}(r):=\underset{\mathbf{z}\in \mathcal{Z},\mathbf{\theta }\in \Theta }{\sup }\left|{\tilde{f}}_r(\mathbf{z};\mathbf{\theta })-f(\mathbf{z};\mathbf{\theta })\right|+\underset{\mathbf{z}\in \mathcal{Z},\mathbf{\theta }\in \Theta }{\sup }\Vert {\nabla }_{\mathbf{\theta }}{\tilde{f}}_r(\mathbf{z};\mathbf{\theta })-{\nabla }_{\mathbf{\theta }}f(\mathbf{z};\mathbf{\theta })\Vert ,$$

$$\sqrt{n}\left({\hat{\mathbf{\theta }}}_{MAL}-{\mathbf{\theta }}_0\right)\stackrel{d}{\to }N(0,{\mathcal{I}}^{-1}).$$

$$C_1(f):=\frac{1+{\sup }_{\mathbf{z},\mathbf{\theta }}\Vert {\nabla }_{\mathbf{\theta }}f(\mathbf{z},\mathbf{\theta })\Vert }{{\delta }^2}$$

$$C_2(f):=2\frac{1+{\sup }_{\mathbf{z},\mathbf{\theta }}\Vert {\nabla }_{\mathbf{\theta }}f(\mathbf{z},\mathbf{\theta }){\Vert }^2+{\sup }_{\mathbf{z},\mathbf{\theta }}\Vert {\nabla }_{\mathbf{\theta },\mathbf{\theta }}f(\mathbf{z},\mathbf{\theta })\Vert }{{\delta }^3}.$$

$$\sqrt{n}\underset{\mathbf{\theta }\in \Theta }{\sup }\Vert {\nabla }_{\mathbf{\theta }}{\tilde{L}}_n(\mathbf{\theta })-{\nabla }_{\mathbf{\theta }}{L}_n(\mathbf{\theta })\Vert$$

$$\le C_1(f)\sqrt{n}\mathcal{E}(R(n)),$$

$$\underset{\mathbf{\theta }\in \Theta }{\sup }\Vert {\nabla }_{\mathbf{\theta }\mathbf{\theta }}{\tilde{L}}_n(\mathbf{\theta })-{\nabla }_{\mathbf{\theta }\mathbf{\theta }}{L}_n(\mathbf{\theta })\Vert$$

$$\displaystyle\leq C_{2}(f)\Bigg{(}\mathcal{E}(R(n))+\sup_{{\boldsymbol{\theta}}\in\Theta,{\mathbf{z}}\in\mathcal{Z}}\left\|\nabla_{{\boldsymbol{\theta}}}\tilde{f}_{R(n)}({\mathbf{z}},{\boldsymbol{\theta}})-\nabla_{{\boldsymbol{\theta}}}f({\mathbf{z}},{\boldsymbol{\theta}})\right\|^{2}$$

$$\displaystyle\quad\quad+\sup_{{\boldsymbol{\theta}}\in\Theta,{\mathbf{z}}\in\mathcal{Z}}\left\|\nabla_{{\boldsymbol{\theta}}{\boldsymbol{\theta}}}\tilde{f}_{R(n)}({\mathbf{z}},{\boldsymbol{\theta}})-\nabla_{{\boldsymbol{\theta}}{\boldsymbol{\theta}}}f({\mathbf{z}},{\boldsymbol{\theta}})\right\|\Bigg{)}.$$

$$\mathcal{E}(r)\le c{r}^{-s}.$$

$$\mathcal{E}(r)\le c\exp \left(-\alpha r^{\beta }\right).$$

$$R(n)\ge \lceil c^{\frac{1}{s}}{n}^{\frac{\gamma }{s}}\rceil .$$

Full text

1 Introduction

Consider classical maximum-likelihood estimation, i.e. the estimated parameter vector $\hat{{\boldsymbol{\theta}}}_{ML}\in\Theta\subset\mathbb{R}^{p}$ is obtained by maximizing

[TABLE]

Here, $f({\mathbf{z}}_{i},{\boldsymbol{\theta}})$ denotes the individual likelihood contribution of the sample ${\mathbf{z}}_{i}\in\mathcal{Z}\subset\mathbb{R}^{q}$. As $n$ tends to infinity, the ML-estimator converges to the true parameter ${\boldsymbol{\theta}}_{0}$, which defines the distribution the sample ${\mathbf{z}}_{1},\ldots,{\mathbf{z}}_{n}$ was drawn from, i.e. $\operatorname*{plim}_{n\to\infty}\hat{{\boldsymbol{\theta}}}_{ML}={\boldsymbol{\theta}}_{0}$.

Often, the function $f$ stems from an integral representation

[TABLE]

where $\varphi:\Omega\times\mathcal{Z}\times\Theta\to\mathbb{R}$ is some (usually non-negative) function that is integrated with respect to the variable ${\mathbf{v}}$ and a weight function $\omega:\Omega\to\mathbb{R}_{+}$ that is defined on a domain $\Omega\subset\mathbb{R}^{d}$. In many relevant cases, these integrals cannot be computed in closed form, e.g. for models of discrete choice (Butler and Moffitt, 1982; McFadden and Train, 2000; Train, 2009), or general limited dependent variable models (Hajivassiliou and Ruud, 1994). Then, $f$ can be approximated with an $r$-point quadrature rule, i.e.

[TABLE]

Here, the quadrature points and weights $({\mathbf{v}}_{j,r},w_{j,r})_{j=1}^{r}$ are independed of ${\mathbf{z}}$ and need to be chosen properly to guarantee a certain accuracy, provided that specific assumptions on $\varphi$ are valid. Moreover, increasing the parameter $r$ allows to increase the accuracy on the one hand, but leads to additional cost on the other hand. Therefore, it is necessary to determine the required accuracy for the approximation of $f$ by $\tilde{f}_{r}$ such that the resulting estimator maintains consistency and asymptotical normality.

To this end, we introduce a link function $R:\mathbb{N}\to\mathbb{N}$ that couples the number of integration points $r$ to the sample size $n$, i.e. $r=R(n)$. Now, $R(n)$ can be chosen such that $\tilde{f}_{R(n)}({\mathbf{z}}_{i},{\boldsymbol{\theta}})$ approximates $f({\mathbf{z}}_{i},{\boldsymbol{\theta}})$ well enough to ignore this additional approximation error but not better, to keep the overall cost at a tractable level. Altogether, this leads to the maximum approximated likelihood estimator (MALE) given by

[TABLE]

Here, the cost for one evaluation of $\tilde{L}_{n}$ in terms of evaluations of $\varphi$ is $n\cdot R(n)$.

A classical and flexible tool to deal with the numerical integration problems in (1.3) is Monte Carlo simulation (MC). Here, all weights are chosen uniformly as $w_{j,r}=1/r,j=1,\ldots,r$ and the points ${\mathbf{v}}_{j,r}$ are drawn identically and independently from the probability distribution induced by the weight function $\omega$. Since MC only requires weak assumptions on $\varphi$, simulation-based estimation has become part of the standard econometrics toolkit. It is implemented in software packages and taught at graduate schools. Comprehensive surveys can be found in Hajivassiliou and Ruud (1994), Gouriéroux and Monfort (1996) and Train (2009).

Maximum simulated likelihood (MSL) is consistent under the usual assumptions if the number of simulation draws $r$ increases with the sample size $n$. In order to achieve asymptotic efficiency for identical sampling, $r$ has to grow at least linearly in $n$, i.e. $r=R(n)\geq n$.111”Identical sampling” denotes the use of one draw-set for all likelihood contributions, i.e. the same quadrature rule is used for all samples ${\mathbf{z}}_{i}$. In contrast, ”independent sampling” denotes the use of one draw-set per likelihood contribution, i.e. for each likelihood contribution an individual quadrature rule is applied. In case of independent sampling, $r$ only has to increase faster than $\sqrt{n}$, see for example Hajivassiliou and Ruud (1994). A recent discussion about the difference between both methodologies can also be found in Kristensen and Salanié (2017).

While asymptotically the simulation error disappears under the appropriate conditions, the computational burden imposed by a sufficient number of simulation draws to achieve well-behaved estimators can be cumbersome or prohibitive in practice. Many studies have found that approximation algorithms other than Monte-Carlo simulation can achieve much higher accuracy and better behaved estimators with a reduced level of computational costs. For one-dimensional problems that are sufficiently smooth, Gaussian quadrature (Butler and Moffitt, 1982) is an obvious choice. For multivariate problems, quasi-Monte-Carlo methods like Halton draws (Bhat, 2001) or other deterministic approaches like integration on sparse grids (Heiss and Winschel, 2008) have been applied successfully.

These deterministic approximation algorithms work well enough in practice to warrant their routine use in software packages like the mixed logit implementation of Stata. But while there is plenty of literature on properties of simulation-based approaches, little is known for estimators based on the other methods. Some papers that use deterministic numerical integration methods other than pure simulation ignore the approximation error in the discussion of the estimator. This applies mainly to examples with one-dimensional integration problems which are tackled with Gaussian quadrature such as Butler and Moffitt (1982). Other papers discuss the well-studied simulation-based estimators before arguing that their deterministic approaches tend to work better in practice, e.g. in Bhat (2001) or in Sándor and Train (2004). This might be due to the fact that the theoretical properties and prerequisites of deterministic approximation schemes are not yet very well understood. Ackerberg et al. (2009) took a first step by giving a set of conditions for the approximated log likelihood contributions that imply consistency and asymptotic normality. However their results remain on the rather abstract level of log-likelihood approximation and provide no guidance on how to check these conditions for specific applications and approximation algorithms. In particular, it remains unclear how the integration error for the approximation of (1.2) by (1.3) influences the statistical properties of the approximated likelihood estimator.

This paper aims at closing this gap and provides a general and unified discussion of the asymptotic properties of a large class of estimators that is based on a broad range of integration methods. The well-known results of simulation-based approaches emerge as special cases. We provide specific conditions under which these estimators are consistent and under which additional conditions the approximation error is irrelevant for their asymptotic distribution. For example, we can derive that a logarithmic growth of the number of samples $r$ is sufficient, if the convergence in $r$ is fast enough. Therefore, in the setting of Gaussian quadrature, it suffices to have $r=R(n)=\log(n)$. This is a huge reduction compared to $R(n)\approx n$ or $R(n)\approx\sqrt{n}$. As an application and demonstration of our framework to a specific model, we deal with mixed logit models and the Butler-Moffitt model, which both are estimated using an approximation that is based on Gauss-Hermite quadrature or Gauss-Hermite sparse grids.

The remainder of this article is organized as follows: In Section 2 we discuss conditions that imply the consistency and asymptotic normality of general extremum estimators. In Section 3 we specialize on maximum approximated likelihood estimators and the assumptions that have to be made such that the conditions from the previous section are fulfilled. Then, in Section 4, we deal with specific integration algorithms, like Gaussian quadrature, quasi–Monte Carlo and sparse grids. In Section 5, we put our theory into practice, by analyzing specific econometric models in the light of our results. This analysis is supplemented by numerical results in Section 6.

2 Asymptotic theory with approximated objective functions

As the most general framework, we consider extremum or M-estimators as discussed in Newey and McFadden (1994), NM hereafter. An M-estimator $\hat{{\boldsymbol{\theta}}}_{M}$ maximizes (or minimizes) some objective function $Q_{n}({\boldsymbol{\theta}})$, i.e.

[TABLE]

where $n$ refers to the number of samples ${\mathbf{z}}_{1},\ldots,{\mathbf{z}}_{n}$ that contribute to $Q_{n}({\boldsymbol{\theta}})$. Examples include least squares, maximum likelihood, GMM, and minimum distance. Throughout this paper, we will maintain the assumption that $\hat{{\boldsymbol{\theta}}}_{M}$ would have all desired properties if we were able to compute the respective objective function $Q_{n}({\boldsymbol{\theta}})$ and therefore the estimator itself. NM comprehensively discuss the conditions to ensure these properties.

We are interested in the problems that arise if $Q_{n}({\boldsymbol{\theta}})$ cannot be evaluated analytically and therefore needs to be approximated. Let $\tilde{Q}_{n}({\boldsymbol{\theta}})$ denote the approximated objective function. Then, our approximated M-estimator is simply

[TABLE]

We now give general conditions to ensure consistency and the asymptotic distribution of $\hat{{\boldsymbol{\theta}}}_{AM}$. We start by recalling Theorem 2.1 of NM:

Lemma 1**.**

Assume that there is a function $Q_{0}({\boldsymbol{\theta}})$ such that (i) $Q_{0}({\boldsymbol{\theta}})$ is uniquely maximized at ${\boldsymbol{\theta}}_{0}$; (ii) ${\boldsymbol{\theta}}\in\Theta$ and $\Theta$ is compact; (iii) $Q_{0}({\boldsymbol{\theta}})$ is continuous; (iv) the function $Q_{n}({\boldsymbol{\theta}})$ converges uniformly in probability to $Q_{0}({\boldsymbol{\theta}})$. Then, $\hat{{\boldsymbol{\theta}}}_{M}=\arg\max_{{\boldsymbol{\theta}}\in\Theta}Q_{n}({\boldsymbol{\theta}})$ is a consistent estimator of ${\boldsymbol{\theta}}_{0}$, i.e. $\operatorname*{plim}_{n\rightarrow\infty}\hat{{\boldsymbol{\theta}}}={\boldsymbol{\theta}}_{0}.$

We abstract from any misspecifications and other problems that could violate the conditions of Lemma 1 in order to focus on the inaccuracies introduced by the approximation of the objective function. Indeed, We can apply the same arguments as in Lemma 1 to the approximated M-estimator $\hat{\boldsymbol{\theta}}_{AM}$ if we can ensure that assumption (iv) also holds for the approximated objective function. The following Theorem 2 states that this is the case as long as $\tilde{Q}_{n}({\boldsymbol{\theta}})$ converges uniformly in probability to the unavailable exact objective function $Q_{n}({\boldsymbol{\theta}})$.

Theorem 2**.**

Assume that (i) the assumptions of Lemma 1 hold (ii) $\tilde{Q}_{n}({\boldsymbol{\theta}})$ converges uniformly in probability to $Q_{n}({\boldsymbol{\theta}})$, i.e. $\operatorname*{plim}_{n\rightarrow\infty}\sup_{{\boldsymbol{\theta}}\in\Theta}\left|\tilde{Q}_{n}({\boldsymbol{\theta}})-Q_{n}({\boldsymbol{\theta}})\right|=0.$ Then, $\hat{{\boldsymbol{\theta}}}_{AM}$ is a consistent estimator of ${\boldsymbol{\theta}}_{0}$, i.e.

[TABLE]

Proof.

We assume in 2(i) that the M-estimator with the exact objective function is consistent. We can use the same arguments as NM to show Lemma 1 if we can establish that $\tilde{Q}_{n}({\boldsymbol{\theta}})$ converges uniformly in probability to $Q_{0}({\boldsymbol{\theta}})$ (Assumption 1(iv)). To see why this holds true, note that, by the triangle inequality for norms there holds

[TABLE]

Both terms converge to zero in probability: the first by assumption 2(ii), the second by assumption 2(i). ∎

Note that assumption 2(ii) requires the approximation accuracy to increase with $n$. In the well-known example of approximation by Monte Carlo simulation, we can increase the number of simulation draws as $n\rightarrow\infty$. We will come back to this more explicitly when we discuss specific approximations approaches.

accuracy parameter $r$ to the number of observations $n$ we introduce a function $R:\mathbb{N}\to\mathbb{N}$. We will assume that $R(n)$ is monotonically increasing, i.e. $R(n+1)\geq R(n)$.

In order to derive the asymptotic distribution of $\hat{\boldsymbol{\theta}}_{AM}$, we again recall the relevant Theorem 3.1 of NM for extremum estimators $\hat{{\boldsymbol{\theta}}}=\operatorname*{arg\,max}_{\boldsymbol{\theta}}{Q}_{n}$.

Lemma 3**.**

Suppose that the assumptions of Lemma 1 hold and that (i) ${\boldsymbol{\theta}}_{0}\in\text{interior}(\Theta)$; (ii) $Q_{n}({\boldsymbol{\theta}})$ is twice continuously differentiable in a neighborhood $\mathcal{N}$ of ${\boldsymbol{\theta}}_{0}$; (iii) $\sqrt{n}\nabla_{{\boldsymbol{\theta}}}Q_{n}({\boldsymbol{\theta}}_{0})\overset{d}{\rightarrow}N(0,\Sigma)$; (iv) there is $H({\boldsymbol{\theta}})$ that is continuous at ${\boldsymbol{\theta}}_{0}$ and $\sup_{{\boldsymbol{\theta}}\in\mathcal{N}}\left\|\nabla_{{\boldsymbol{\theta}}{\boldsymbol{\theta}}}Q_{n}({\boldsymbol{\theta}})-H({\boldsymbol{\theta}})\right\|\overset{p}{\rightarrow}0$; (v) $H=H({\boldsymbol{\theta}}_{0})$ is nonsingular. Then the M-estimator is asymptotically normal distributed, i.e $\sqrt{n}(\hat{{\boldsymbol{\theta}}}-{\boldsymbol{\theta}}_{0})\overset{d}{\rightarrow}N(0,H^{-1}\Sigma H^{-1})$.

To analyze our approximate M-estimator, we again assume that the conditions of this lemma hold and then give additional assumptions such that the same results apply to the approximate M-estimator $\hat{{\boldsymbol{\theta}}}_{AM}$.

Theorem 4**.**

Assume that (i) the assumptions of Lemma 3 hold; (ii) With probability one, $\tilde{Q}_{n}({\boldsymbol{\theta}})$ is twice continuously differentiable in a neighborhood $\mathcal{N}$ of ${\boldsymbol{\theta}}_{0}$; (iii) $\operatorname*{plim}_{n\rightarrow\infty}\sup_{{\boldsymbol{\theta}}\in\Theta}\left|\tilde{Q}_{n}({\boldsymbol{\theta}})-Q_{n}({\boldsymbol{\theta}})\right|=0$; (iv) $\operatorname*{plim}_{n\rightarrow\infty}\sqrt{n}\,\sup_{{\boldsymbol{\theta}}\in\Theta}\left\|\nabla_{{\boldsymbol{\theta}}}\tilde{Q}_{n}({\boldsymbol{\theta}})-\nabla_{{\boldsymbol{\theta}}}Q_{n}({\boldsymbol{\theta}})\right\|=0$; (v) $\operatorname*{plim}_{n\rightarrow\infty}\sup_{{\boldsymbol{\theta}}\in\Theta}\left\|\nabla_{{\boldsymbol{\theta}}{\boldsymbol{\theta}}}\tilde{Q}_{n}({\boldsymbol{\theta}})-\nabla_{{\boldsymbol{\theta}}{\boldsymbol{\theta}}}Q_{n}({\boldsymbol{\theta}})\right\|=0$. Then

[TABLE]

so $\hat{{\boldsymbol{\theta}}}_{AM}$ has the same limiting distribution as the infeasible $\hat{{\boldsymbol{\theta}}}_{M}$.

Proof.

We can apply the same arguments used by NM to show Lemma 3, applied to our approximated objective function $\tilde{Q}_{n}({\boldsymbol{\theta}})$. By our assumptions 4(i) and 4(ii), assumptions 3(i), 3(ii), and 3(v) are directly implied. To check 3(iii), we write

[TABLE]

The first term converges in distribution to $N(0,\Sigma)$ by assumption 4(i) and 3(iii). The second term converges in probability to zero by 4(iv). Finally, we confirm 3(iv) by noting that

[TABLE]

is implied by the triangle inequality. Both terms converge to zero in probability: the first by assumption 4(v), the second by assumption 4(i). ∎

Arguably the most important condition for deriving the asymptotic distribution is assumption 4(iv). It not only requires the approximated gradient to converge to the exact value, but the rate of convergence also needs to be faster than $1/\sqrt{n}$.

3 Maximum approximated likelihood

In the last section we discussed general extremum estimators with approximated objective functions. These results are similar to those of Ackerberg et al. (2009) who study general approximation algorithms for maximum likelihood estimation. For the remainder of this paper, we focus on the case of maximum likelihood. This covers a large share of the applications of approximate M-estimation and allows us to be more specific. We will derive general conditions for the likelihood contributions and the approximation algorithm to ensure favorable properties of the resulting MAL estimator.

3.1 Asymptotic analysis with respect to likelihood contributions

We consider an i.i.d. random sample $[{\mathbf{z}}_{i}\in\mathcal{Z};\ i=1,\dots,n]$ from a population distribution characterized by the family of probability mass or density functions $f({\mathbf{z}};{\boldsymbol{\theta}}_{0})$ and the sample space $\mathcal{Z}\subset\mathbb{R}^{d}$. Here, ${\mathbf{z}}$ includes all variables. In most econometric models, ${\mathbf{z}}$ includes some “endogenous” variables ${\mathbf{y}}$ and some “exogenous” variables ${\mathbf{x}}$. In these cases, $f({\mathbf{z}};{\boldsymbol{\theta}}_{0})$ is actually conditional on ${\mathbf{x}}$. For notational convenience and consistency with the literature, we will not explicitly make this distinction. The exact log likelihood function is

[TABLE]

where the individual likelihood contributions of sample ${\mathbf{z}}_{i}$ are given by $f({\mathbf{z}}_{i};{\boldsymbol{\theta}})$. The maximum likelihood estimator $\hat{{\boldsymbol{\theta}}}_{ML}$ maximizes $L_{n}({\boldsymbol{\theta}})$, i.e. $\hat{{\boldsymbol{\theta}}}_{ML}=\operatorname*{arg\,max}_{\boldsymbol{\theta}}L_{n}({\boldsymbol{\theta}})$. Moreover, in the context of the preceding section, we have $\operatorname*{plim}_{n\to\infty}Q_{n}({\boldsymbol{\theta}})=Q_{0}({\boldsymbol{\theta}})$, where $Q_{0}$ was introduced in Lemma 1 and is maximized by the true parameter ${\boldsymbol{\theta}}_{0}$.

Newey and McFadden (1994, Theorems 2.5, 3.3) provide conditions for the ML estimator to have favorable properties like consistency and asymptotic normality. We recall them in the following, treating the consistency of the ML estimator first.

Lemma 5**.**

Assume that (i) For all ${\boldsymbol{\theta}}\neq{\boldsymbol{\theta}}_{0}$, we have $f({\mathbf{y}};{\boldsymbol{\theta}})\neq f({\mathbf{y}};{\boldsymbol{\theta}}_{0})$; (ii) ${\boldsymbol{\theta}}_{0}\in\Theta$ and $\Theta$ is compact; (iii) $\log f({\mathbf{y}};{\boldsymbol{\theta}})$ is continuous at each ${\boldsymbol{\theta}}\in\Theta$ with probability one; (iv) $\mathbb{E}_{\mathbf{y}}[\sup_{{\boldsymbol{\theta}}\in\Theta}|\log f({\mathbf{y}};{\boldsymbol{\theta}})|]<\infty$. Then, $\operatorname*{plim}_{n\rightarrow\infty}\hat{{\boldsymbol{\theta}}}_{ML}={\boldsymbol{\theta}}_{0}$ (consistency).

Next, we recall the following result on the asymptotic normality of the ML estimator.

Lemma 6**.**

Assume that (i) the assumptions of Lemma 5 hold and that (ii) ${\boldsymbol{\theta}}_{0}\in\text{interior}(\Theta)$; (iii) $f({\mathbf{y}};{\boldsymbol{\theta}})$ is twice continuously differentiable and $f({\mathbf{y}};{\boldsymbol{\theta}})>0$ in a neighborhood $\mathcal{N}$ of ${\boldsymbol{\theta}}_{0}$; (iv) $\int\sup_{\boldsymbol{\theta\in\mathcal{N}}}\left\|\nabla_{{\boldsymbol{\theta}}}f({\mathbf{y}};{\boldsymbol{\theta}})\right\|\,d{\mathbf{y}}<\infty$ and $\int\sup_{\boldsymbol{\theta\in\mathcal{N}}}\left\|\nabla_{{\boldsymbol{\theta}}{\boldsymbol{\theta}}}f({\mathbf{y}};{\boldsymbol{\theta}})\right\|\,d{\mathbf{y}}<\infty$; (v) $\mathcal{I}=\mathbb{E}_{\mathbf{y}}\left[\nabla_{{\boldsymbol{\theta}}}f({\mathbf{y}};{\boldsymbol{\theta}}_{0})\left(\nabla_{{\boldsymbol{\theta}}}f({\mathbf{y}};{\boldsymbol{\theta}}_{0})\right)^{T}\right]$ exists and is nonsingular, (vi) $\mathbb{E}_{\mathbf{y}}\left[\sup_{\boldsymbol{\theta\in\mathcal{N}}}\left\|\nabla_{{\boldsymbol{\theta}}}\log f({\mathbf{y}};{\boldsymbol{\theta}}_{0})\right\|\right]<\infty$. Then $\sqrt{n}\left(\hat{{\boldsymbol{\theta}}}_{ML}-{\boldsymbol{\theta}}_{0}\right)\overset{d}{\rightarrow}N\left(0,\mathcal{I}^{-1}\right)$ (asymptotic normality).

If the likelihood contributions $f({\mathbf{z}}_{i};{\boldsymbol{\theta}})$ cannot be computed analytically, we need to approximate them using some algorithm. The approximated likelihood contributions are denoted by $\tilde{f}_{r}({\mathbf{z}}_{i};{\boldsymbol{\theta}}),i=1,\dots,n$ and depend on an accuracy parameter $r$, but not on the sample ${\mathbf{z}}_{i}$, i.e. the same approximation approach to $f$ is used for all likelihood contributions. For maximum simulated likelihood, $r$ might be the number of simulation draws. In general, it determines the amount of approximation error as well as the computational costs.

Theorems 2 and 4 require the approximated objective function $\tilde{Q}_{n}({\boldsymbol{\theta}})$ to converge to the exact function $Q_{n}({\boldsymbol{\theta}})$ as $n\rightarrow\infty$. In general, this requires $r$ to increase with $n$. Therefore, it is not sufficient to keep $r$ fixed for all $n$, but to consider a sequence of approximation functions $(\tilde{f}_{r}({\mathbf{z}}_{i};{\boldsymbol{\theta}}))_{r=1}^{\infty}$ which, in a sense that will be clarified later, converge to the exact $f({\mathbf{z}}_{i};{\boldsymbol{\theta}})$ as $r\rightarrow\infty$. Then, the approximated log-likelihood

[TABLE]

fulfills $\lim_{n,r\to\infty}\tilde{L}_{n,r}({\boldsymbol{\theta}})=Q_{0}({\boldsymbol{\theta}})$. In practice for any finite $n\in\mathbb{N}$, we need to choose a finite accuracy level $r$. Choosing $r$ too large will result in unneccesary cost, while using a too small $n$ will result in an additional error contribution. This trade-off will be the subject of the remainder of this section.

In order to link the accuracy parameter $r$ to the number of observations $n$ we introduce a function $R:\mathbb{N}\to\mathbb{N}$. We will assume that $R(n)$ is monotonically increasing, i.e. $R(n+1)\geq R(n)$ and $\lim_{n\to\infty}R(n)=\infty$. The speed with which $R(n)$ needs to increase with $n$ will depend on the specific approximation algorithms as will be discussed below.

Then, the maximum approximated likelihood estimator $\hat{{\boldsymbol{\theta}}}_{MAL}$ maximizes

[TABLE]

We now discuss general properties of $\tilde{f}_{R(n)}({\mathbf{z}}_{i};{\boldsymbol{\theta}})$ that imply the consistency and the asymptotic distribution of $\hat{\boldsymbol{\theta}}_{MAL}$. In order to focus on the effects of approximation, we assume that the infeasible ML estimator would have all desired properties if it were available, i.e. that the assumptions of Lemma 5 and 6 hold.

In addition to the consistency of the ML estimator, the only additional assumption we need for consistency of $\hat{{\boldsymbol{\theta}}}_{MAL}$ according to Theorem 2 is the uniform convergence in probability of the approximation $\tilde{L}_{n}({\boldsymbol{\theta}})$ to the true, but intractable $L_{n}({\boldsymbol{\theta}})$.

Theorem 7**.**

Assume that (i) the conditions of Lemma 5 hold and that (ii) there exists $\bar{\delta}>0$ such that for all ${\mathbf{z}}\in\mathcal{Z}$ and all ${\boldsymbol{\theta}}\in\Theta$ it holds that $f({\mathbf{z}},{\boldsymbol{\theta}})\geq\bar{\delta}$; (iii) $\lim_{r\to\infty}\sup_{{\mathbf{z}}\in\mathcal{Z},{\boldsymbol{\theta}}\in\Theta}\left|\tilde{f}_{r}({\mathbf{z}},{\boldsymbol{\theta}})-f({\mathbf{z}},{\boldsymbol{\theta}})\right|=0$; (iv) $R(n)$ is monotonically increasing in $n$ and $\lim_{n\to\infty}R(n)=\infty$. Then, $\hat{{\boldsymbol{\theta}}}_{MAL}$ is a consistent estimator of ${\boldsymbol{\theta}}_{0}$, i.e.

[TABLE]

Proof.

We use our general results from Theorem 2. Given the (infeasible) MLE is consistent by assumption (i), we only have to verify Assumption 2(ii).

Note at this point that there exists $r_{0}\in\mathbb{N}$ and a real number $\delta\in(0,\bar{\delta})$ such that for all $r\geq r_{0}$, all ${\mathbf{z}}\in\mathcal{Z}$ and all ${\boldsymbol{\theta}}\in\Theta$ it holds that both, $\tilde{f}_{r}({\mathbf{z}},{\boldsymbol{\theta}})\geq\delta$ and $f({\mathbf{z}},{\boldsymbol{\theta}})\geq\bar{\delta}\geq\delta$. To see why this is true, we use that there exists $\bar{\delta}>0$ such that for all ${\boldsymbol{\theta}}\in\Theta$ and ${\mathbf{z}}\in\mathcal{Z}$ it holds $f({\boldsymbol{\theta}},{\mathbf{z}})>\bar{\delta}$. Then, we choose some $\delta$ with $0<\delta<\bar{\delta}$. Because of 7(iii), we have $\lim_{r\to\infty}\tilde{f}_{r}({\mathbf{z}};{\boldsymbol{\theta}})=f({\mathbf{z}};{\boldsymbol{\theta}})$ for all ${\mathbf{z}}\in\mathcal{Z}$ and ${\boldsymbol{\theta}}\in\Theta$. Hence, there exists $r_{0}(\delta)\in\mathbb{N}$ such that for all $r\geq r_{0}(\delta)$ it holds for all ${\mathbf{z}}\in\mathcal{Z}$ and all ${\boldsymbol{\theta}}\in\Theta$ that

[TABLE]

This implies $\tilde{f}_{r}({\mathbf{z}};{\boldsymbol{\theta}})\geq\delta$ for all $r\geq r_{0}(\delta)$.

Now we are in the position to show

[TABLE]

To this end, we write

[TABLE]

For sufficiently large $n$ (i.e. $R(n)\geq r_{0}$), we can use Corollary 19 (i) in Appendix A. By assumptions (i) and (ii) with $h({\boldsymbol{\theta}})=f({\mathbf{z}}_{i},{\boldsymbol{\theta}})$ and $g({\boldsymbol{\theta}})=\tilde{f}_{R(n)}({\mathbf{z}}_{i},{\boldsymbol{\theta}})$ it follows that

[TABLE]

Moreover, by assumption (iii) and (iv), this expression converges to zero in probability. ∎

In order to study the asymptotic distribution of the MAL estimator, we assume that the (intractable) ML estimator $\hat{\boldsymbol{\theta}}_{ML}$ is asymptotically normally distributed and efficient and that the MAL estimator $\hat{\boldsymbol{\theta}}_{MAL}$ is consistent. Then we make additional assumptions to derive the asymptotic distribution of $\hat{\boldsymbol{\theta}}_{MAL}$. To this end, we define the quantity

[TABLE]

which measures the worst-case approximation error of both, the function $f$ and its gradient $\nabla f$ by $\tilde{f}_{r}$ and $\nabla\tilde{f}_{r}$, respectively.

Theorem 8**.**

Assume that the assumptions of Lemma 6 and Theorem 7 hold. Also assume that (i) $\tilde{f}_{r}({\mathbf{z}};{\boldsymbol{\theta}})$ is twice continuously differentiable in ${\boldsymbol{\theta}}\in\Theta$ and bounded for all ${\mathbf{z}}\in\mathcal{Z}$ and $r\in\mathbb{N}$; (ii) $\operatorname*{plim}_{n\to\infty}\sqrt{n}\,\mathcal{E}(R(n))=0$; (iii) $\sup_{{\mathbf{z}}\in\mathcal{Z},{\boldsymbol{\theta}}\in\Theta}\left\|\nabla_{{\boldsymbol{\theta}}{\boldsymbol{\theta}}}\tilde{f}_{R(n)}({\mathbf{z}};{\boldsymbol{\theta}})-\nabla_{{\boldsymbol{\theta}}{\boldsymbol{\theta}}}{f}({\mathbf{z}};{\boldsymbol{\theta}})\right\|\overset{p}{\rightarrow}0$. Then

[TABLE]

Proof.

We proceed by verifying the assumptions of Theorem 4. Assumption 4(i) directly follows from Assumption 8(i).

The critical conditions we have to check are 4(iv) and 4(v). The arguments are similar to that in the proof of Theorem 7: Again, there exists a $\delta\in(0,\bar{\delta})$ such that for sufficiently large $r$ it holds $\tilde{f}_{r}>\delta$. In order to show 4(iv), we define

[TABLE]

and

[TABLE]

Then, we use Corollary 19 (ii) and the definition of $\mathcal{E}_{r}$ to derive

[TABLE]

which converges to zero in probability by Assumption 8(iv).

Using similar arguments, we employ Corollary 19(iii) to verify 4(v). We compute

[TABLE]

Clearly, as $n\to\infty$, by the given assumptions, all three summands within the bracket tend to zero in probability. This implies that the whole expression tends to zero, because $C_{2}(f)$ is independent of $n$.

∎

3.2 The level of approximation accuracy

In Theorem 8 it is required by condition (ii) that $\operatorname*{plim}_{n\to\infty}\sqrt{n}\,\mathcal{E}(R(n))=0$, i.e. the approximation error (3.7) of the gradient $\mathcal{E}(R(n))$ needs to decay faster than $n^{-1/2}$. Different approximation algorithms behave differently in terms of how the error bounds $\mathcal{E}(r)$ change with $r$. We discuss the two most common forms of convergence and their implication on the choice of $R:\mathbb{N}\to\mathbb{N}$.

• (a)

Algebraic convergence: For constants $c>0$ and $s>0$ it holds

[TABLE]

• (b)

Exponential convergence: For constants $c>0$ and $\alpha,\beta>0$ it holds

[TABLE]

Now, we give general results on how to choose $R(n)$ depending on the convergence rates of the approximation algorithm for $\tilde{f}_{r}$.

Theorem 9**.**

1. (a)

Assume that an algebraic convergence rate $\mathcal{E}(r)\leq cr^{-s}$ holds for sufficiently large $r\in\mathbb{N}$ and $s>0$. Then, for all $\gamma>\frac{1}{2}$ it holds that $\sqrt{n}\mathcal{E}(R(n))\overset{p}{\rightarrow}0$ if

[TABLE] 2. (b)

Assume that an exponential convergence rate \mathcal{E}(r)\leq c\exp\big{(}-\alpha r^{\beta}\big{)} holds for sufficiently large $r\in\mathbb{N}$ and $\alpha,\beta>0$. Then, for all $\gamma>\frac{1}{2}$ it holds that $\sqrt{n}\mathcal{E}(R(n))\overset{p}{\rightarrow}0$ if

[TABLE]

Proof.

For (a), we insert (3.10) into (3.8). This yields

[TABLE]

which, for $\gamma>1/2$ tends to zero as $n\to\infty$.

For (b), we insert (3.11) into (3.9) and obtain

[TABLE]

which also tends to zero as $n\to\infty$, if $\gamma>1/2$. ∎

4 Specific approximation algorithms

Up to this point, we did not specify any particular approximation. In the following, we examine a very common setting, i.e. a likelihood contribution of the form

[TABLE]

which is a (possibly multivariate) integral of a function $\varphi:\Omega\times\mathcal{Z}\times\Theta\to\mathbb{R}$ over the domain of integration $\Omega\subset\mathbb{R}^{d}$ with respect to a given weight function $\omega:\Omega\to\mathbb{R}_{+}$. The integral typically represents the expectation over a nonlinear function $\varphi$ of a set of random variables with density function $\omega$. This class of models includes nonlinear random effects models like Butler and Moffitt (1982) and random coefficients models like the mixed logit model, see for example McFadden and Train (2000).

The need for approximation arises because the integral in (4.1) often cannot be computed in closed form. To this end, several estimators have been proposed including the method of simulated moments (MSM; McFadden (1989)) and the method of maximum simulated scores (MSS; Hajivassiliou and McFadden (1998)). Possibly the most widely used approach is the maximum simulated likelihood (MSL) estimator. It approximates $f({\mathbf{z}}_{i};{\boldsymbol{\theta}})$ with a Monte Carlo estimate, i.e. $\tilde{f}_{r}({\mathbf{z}}_{i};{\boldsymbol{\theta}})=\frac{1}{r}\sum_{j=1}^{r}\varphi({\mathbf{v}}_{j},{\mathbf{z}}_{i},{\boldsymbol{\theta}})$, where $({\mathbf{v}}_{j})_{j=1}^{r}$ is a set of random draws distributed according to the weight function $\omega$.

Many authors have found that integration rules other than pure Monte-Carlo simulation often perform much better in practice. Examples include Quasi-Monte Carlo rules in Bhat (2001), Gaussian quadrature in Butler and Moffitt (1982), or quadrature on sparse grids in Heiss and Winschel (2008).

To cover all the mentioned approaches, we consider a general approximation222Note that both, the weights $w_{j,r}$ and the points ${\mathbf{v}}_{j,r}$ do not depend on the sample ${\mathbf{z}}_{i}$. This implies that for each likelihood contribution the same approximation algorithm is employed.

[TABLE]

to the true integral of (4.1). This formulation includes Monte-Carlo, where all weights ($\omega_{j,r}$) are equal to $r^{-1}$ and the nodes (${\mathbf{v}}_{j,r}$) are random draws, but also Quasi-Monte Carlo (QMC; like Halton or Sobol sequences), cf. Halton (1964) or Sobol (1967), where the weights are also equal to $r^{-1}$ but the nodes are deterministic. Moreover, it includes classical (Gaussian) quadrature rules, cf. Davis and Rabinowitz (2007), as well as sparse grids (Gerstner and Griebel, 1998). We will come back to specific algorithms below.

Remember that for continuous $\varphi(\cdot,{\mathbf{z}},{\boldsymbol{\theta}})$ our likelihood contributions $f({\mathbf{z}},{\boldsymbol{\theta}})$ and their derivatives now have the following general form

[TABLE]

Their associated approximations $\tilde{f}_{r}({\mathbf{z}},{\boldsymbol{\theta}})$ stem from some cubature rule $U_{r}$ of the form

[TABLE]

The goal in this section will be to impose conditions on $\varphi$ as well as on the integration rule $U_{r}$ such that the assumptions of Theorems 7 and 8 hold. To this end, the following notation for the partial derivative of a sufficiently smooth function $g$ will be helpful.

[TABLE]

where ${\boldsymbol{\alpha}}=(\alpha_{1},\ldots,\alpha_{p})$ is a multi-index that contains the order of the partial derivatives for different coordinate directions $\theta_{1},\ldots,\theta_{p}$.

In order to bound the magnitude of $\sup_{{\mathbf{z}}\in\mathcal{Z},{\boldsymbol{\theta}}\in\Theta}\left|\tilde{f}_{r}({\mathbf{z}},{\boldsymbol{\theta}})-f({\mathbf{z}},{\boldsymbol{\theta}})\right|$, as well as $\sup_{{\mathbf{z}}\in\mathcal{Z},{\boldsymbol{\theta}}\in\Theta}\|\nabla_{{\boldsymbol{\theta}}}\tilde{f}_{r}({\mathbf{z}};{\boldsymbol{\theta}})-\nabla_{{\boldsymbol{\theta}}}f({\mathbf{z}};{\boldsymbol{\theta}})\|$ and $\sup_{{\mathbf{z}}\in\mathcal{Z},{\boldsymbol{\theta}}\in\Theta}\|\nabla_{{\boldsymbol{\theta}}{\boldsymbol{\theta}}}\tilde{f}_{r}({\mathbf{z}};{\boldsymbol{\theta}})-\nabla_{{\boldsymbol{\theta}}{\boldsymbol{\theta}}}f({\mathbf{z}};{\boldsymbol{\theta}})\|$ defined by (4.3) and (4.4), we resort to

[TABLE]

for $k=0,1,2$. Since all matrix- and vector-norms can be bounded by the $\ell_{\infty}$ norm it holds $\mathcal{E}(r)\leq 2p\bar{\mathcal{E}}_{1}(r)$ and hence it is sufficient (and convenient) to work with the quantity defined in (4.6).

Now, $\lim_{r\to\infty}\bar{\mathcal{E}}_{k}(r)=0$ basically implies that the quadrature rules $(U_{r})_{r\in\mathbb{N}}$ in (4.4) converge for $\varphi(\cdot,{\mathbf{z}},{\boldsymbol{\theta}})$ and all partial derivatives with respect to ${\boldsymbol{\theta}}$ up to order $k$.

4.1 Consistency and asymptotic normality

First, we deal with consistency, which follows as a simple application of Theorem 4 to the specific setting where $f$ has the form (4.3) and $\tilde{f}_{r}$ has the form (4.4).

Corollary 10**.**

Assume that the conditions of Lemma 1 hold and that (i) there exists $\bar{\delta}>0$ such that for all values of ${\mathbf{z}}\in\mathcal{Z}$ and ${\boldsymbol{\theta}}\in\Theta$ it holds $f({\mathbf{z}};{\boldsymbol{\theta}})>\bar{\delta}>0$; (ii) the function $\varphi({\mathbf{v}},{\mathbf{z}},{\boldsymbol{\theta}})$ is continuous in ${\boldsymbol{\theta}}\in\Theta$; (iii) $\operatorname*{plim}_{r\to\infty}\bar{\mathcal{E}}_{0}(r)=0$; (iv) $R(n)$ is monotonically increasing in $n$. Then, $\hat{{\boldsymbol{\theta}}}_{MAL}$ is a consistent estimator of ${\boldsymbol{\theta}}_{0}$, i.e.

[TABLE]

Condition 10(iii) states that the cubature rule $U_{r}$ converges for $\varphi$, which is a rather mild assumption since no requirements are made regarding the speed of convergence.

Next, we consider the asymptotic distribution of $\hat{{\boldsymbol{\theta}}}_{MAL}$. Again, $f$ has the form (4.3) and $\tilde{f}_{r}$ has the form (4.4). The critical part of the properties of $\hat{\boldsymbol{\theta}}_{MAL}$ is the worst-case approximation error, i.e. (4.6) for $k=1$, which not only has to go to zero, but also has to vanish at a certain speed.

Corollary 11**.**

Assume that the assumptions of Lemma 6 and Corollary 10 hold. Also assume that (i) for all ${\mathbf{v}}\in\Omega$ and all ${\mathbf{z}}\in\mathcal{Z}$ it holds that $\varphi({\mathbf{v}},{\mathbf{z}},\cdot)$ is twice continuously differentiable; (ii) for all ${\mathbf{z}}\in\mathcal{Z}$ and all ${\boldsymbol{\theta}}\in\Theta$ all terms in (4.3) remain bounded. (iv) $\operatorname*{plim}_{r\to\infty}\bar{\mathcal{E}}_{2}(r)=0$; (v) $R(n)$ rises fast enough with $n$ to ensure \sqrt{n}\bar{\mathcal{E}}_{1}\big{(}R(n)\big{)}\overset{p}{\rightarrow}0. Then

[TABLE]

Remark 12**.**

Note that at no point the positivity of the integration weights is assumed. Only (sufficiently fast) convergence of $U_{r}$ for all partial derivatives of $\varphi$ up to order $2$ is required.

Remark 13**.**

Note that the results obtained in the last three Sections are of asymptotical nature, i.e. they hold for $n\to\infty$. In practice, $n$ and also $r=R(n)$ need to be finite numbers, where $R(n)$ has to be chosen such that $R(n)\geq r_{0}(\delta)$ because otherwise $\tilde{f}_{R(n)}>\delta$ might not hold true. Especially for $\tilde{f}_{R(n)}\leq 0$ the whole approach does not work anymore because of the logarithm in the log-likelihood contributions.

4.2 Examples for integration rules

In this section we will turn from the previous abstract framework to more explicit approximation algorithms and likelihood constructions. We will discuss several choices for the approximation of the integrals (4.3) and their resulting complexities. To this end, remember that the likelihood function consists of $n$ terms where each one requires the numerical approximation of an integral with $r=R(n)$ evaluations of $\varphi$. Therefore, the total complexity for one evaluation of the approximated likelihood function is $n\cdot R(n)$, with the requirements on $R(n)$ given in Theorem 9. Several concrete examples are provided in Table 1.

4.2.1 Monte Carlo simulation

In order verify the consistency of the Monte Carlo simulation (MC) approach it only is required that $\varphi(\cdot,{\mathbf{z}},{\boldsymbol{\theta}})$ is uniformly bounded in $L^{2}(\Omega,\omega)$ for all data ${\mathbf{z}}\in\mathcal{Z}$ and all parameters ${\boldsymbol{\theta}}\in\Theta$. This means that $\sup_{{\mathbf{z}},{\boldsymbol{\theta}}}\int_{\Omega}\varphi({\mathbf{v}},{\mathbf{z}},{\boldsymbol{\theta}})^{2}\,\omega({\mathbf{v}})\mathrm{d}{\mathbf{v}}<\infty$. Then, the choice $w_{j,r}=1/r$ and independent and identical draws ${\mathbf{v}}_{j,r}$ from the probability distribution induced by $\omega$ yield a quadrature rule that converges both, in expectation and with high probability at a rate of $\mathcal{O}(r^{-1/2})$. By Corollary 10 the resulting estimator is consistent.

For asymptotical normality, also the partial derivatives of $\varphi$ with respect to ${\boldsymbol{\theta}}$ must be square integrable, i.e.

[TABLE]

Then, with high probability and also in expectation, $\bar{\mathcal{E}}_{1}\leq cr^{-1/2}$ and $\bar{\mathcal{E}}_{1}$ converges to zero at a rate of $\mathcal{O}(r^{-1/2})$. It follows by Corollary 11 and Theorem 9 that the Monte Carlo based MAL estimator is asymptotically normal if the number of integration points $r$ increases at least linearly in the number of data samples $n$, i.e. $r=R(n)=n^{\beta}$ with $\beta>1$.

4.2.2 Quasi – Monte Carlo integration

Another approach to multivariate integration on the $d$-dimensional unit cube $\Omega=(0,1)^{d}$ is the quasi–Monte Carlo method. Similar to Monte Carlo it has all weights $w_{j}=r^{-1}$, but the points are not drawn randomly, but determined deterministically by certain number-theoretic considerations. Quasi–Monte Carlo methods exist in differenct variants, depending on the specific selection of points. Examples are Halton points, cf. Halton (1964), lattice rules, cf. Sloan and Joe (1994) or Sobol points, cf. Sobol (1967). A rather recent development are higher-order QMC points, cf. Dick and Pillichshammer (2010); Hinrichs et al. (2016).

For integrands that possess bounded variation in the sense of Hardy and Krause, cf. Niederreiter (1992), most QMC constructions achieve a convergence rate of order $\mathcal{O}(r^{-1}\log(r)^{d-1})$. Asymptotically, this bound behaves like $\mathcal{O}(r^{-1+\varepsilon})$, where the $\varepsilon>0$ asymptotically suppresses the $d$-dependent power of $\log r$.

Details on the definition of the Hardy-Krause variation can be found in Dick and Pillichshammer (2010); Niederreiter (1992); Owen (2005). We just note at this point that bounded variation is a much stronger assumption than bounded variance (as it is required for Monte Carlo integration) since for bounded Hardy-Krause variation also the mixed derivatives need to exist and to be bounded, i.e.

[TABLE]

If this property is fulfilled for $g({\mathbf{v}})=\varphi({\mathbf{v}},{\mathbf{z}},{\boldsymbol{\theta}})$, the resulting estimator will be consistent. If in addition also (4.8) is fulfilled for all partial derivatives of $\varphi$ up to order $2$ with respect to ${\boldsymbol{\theta}}$, i.e. $g({\mathbf{v}})=D^{({\boldsymbol{\theta}})}_{\boldsymbol{\alpha}}\varphi({\mathbf{v}},{\mathbf{z}},{\boldsymbol{\theta}}),|{\boldsymbol{\alpha}}|_{1}\leq 2$ then both, $\bar{\mathcal{E}}_{1}(r)$ and $\bar{\mathcal{E}}_{2}(r)$ decay at a rate of $\mathcal{O}(r^{-1+\varepsilon})$ and Corollary 11 and Theorem 9 yield that the number of QMC integration points needs to increase as $r=R(n)=n^{\beta}$ with any $\beta>1/2$.

We remark that in some cases the condition of bounded Hardy-Krause variation can be relaxed to deal with integrands that possess mild singularities, cf. Owen (2006).

4.2.3 Gaussian quadrature and related methods

The classical approach to numerical integration aims at the so-called degree of exactness. This means that the quadrature rule is constructed such that it integrates polynomials up to a certain degree $D$ exactly, i.e. the points $v_{1,r},\ldots,v_{r,r}$ and weights $w_{1,r},\ldots,w_{r,r}$ are determined such that

[TABLE]

The quadrature rule with the best possible degree of exactness $2r-1$ is Gaussian quadrature. Depending on the smoothness of the integrand, Gaussian quadrature yields algebraic convergence rates (for integrands with finite smoothness) or (sub-)exponential convergence rates for infinitely differentiable or analytic integrands, cf. Davis and Rabinowitz (2007). We will discuss two examples for Gaussian quadrature in more detail below. Note that it is also possible to extend the Gaussian approach to non-polynomial basis functions, e.g. to deal with certain boundary singularities that occur with the GHK sampling approach, cf. Griebel and Oettershagen (2014). Moreover, there are further approaches that yield quadrature rules with a polynomial degree of exactness that involve nested point sets. Examples are Gauss-Patterson quadrature rules, cf. Patterson (1968), Clenshaw-Curtis quadrature, cf. Imhof (1963) or Leja points, cf. Jantsch et al. (2016); Griebel and Oettershagen (2016).

Besides favorable convergence rates, Gaussian quadrature rules have the property that the Stone-Weierstrass Theorem and its variations ensure their convergence for any continuous integrand. Therefore, by Theorem 10 maximum approximated likelihood estimators based on Gaussian quadrature rules and other stable polynomial-based quadrature rules usually are consistent if $\varphi(\cdot,{\mathbf{z}},{\boldsymbol{\theta}})$ is continuous for all $({\mathbf{z}},{\boldsymbol{\theta}})\in\mathcal{Z}\times\Theta$.

Gauss-Legendre quadrature

If the weight function is constant, i.e. $\omega(x)=1$ and $\Omega$ is bounded, e.g. $\Omega=[0,1]$, Gauss-Legendre quadrature achieves the best possible degree of polynomial exactness. As a consequence, it achieves exponential convergence rates of type $\mathcal{O}(e^{-\alpha r}),\alpha>0$ for integrands that are analytic in certain ellipses or circles that enclose the domain of integration $\Omega$, cf. Davis and Rabinowitz (2007). However, if the integrand is only $k<\infty$ times differentiable, the convergence rate deteriorates to the order $\mathcal{O}(r^{-k})$.

Gauss-Hermite quadrature

If the weight function corresponds to a standard normal density, e.g. $\omega(x)=\frac{1}{\sqrt{2\pi}}e^{-x^{2}/2}$ and $\Omega=\mathbb{R}$, then Gauss-Hermite quadrature is optimal with respect to polynomial exactness. The analysis of Gaussian quadrature on unbounded domains $\Omega$ is more complicated than in the case of bounded $\Omega$. Yet, there are a number of useful results available: For integrands with $k$ continuous and integrable derivatives, the error of Gauss-Hermite quadrature can be bounded by $\mathcal{O}(r^{-k/2})$, cf. Smith et al. (1983) or Mastroianni and Monegato (1994). To be more precise, if there exists a constant $c>0$ such that

[TABLE]

then it holds

[TABLE]

where the constant $C(k,g)>0$ may depend on $k$ and the integrand $g$, but not on $r$. Basically, condition (4.10) bounds the growth of the $k$-th derivative along the real axis. Moreover, for integrands that are analytic in an infinite complex strip that contains the real axis $\mathbb{R}$, sub-exponential convergence rates of type $\mathcal{O}(e^{-\alpha n^{\beta}})$ are shown in Boyd (1984). Finally, for certain classes of entire functions, even exponential convergence $\mathcal{O}(e^{-\alpha r})$ is possible, cf. Kuo and Woźniakowski (2011) or Irrgeher et al. (2015).

4.2.4 Sparse grid cubature

For multivariate integration problems it is possible to use sparse grids, cf. Gerstner and Griebel (1998, 2003); Griebel and Oettershagen (2016); Novak and Ritter (1996) to turn a univariate quadrature rule to a multivariate integration method.

For integrands with bounded mixed derivatives up to order $k$ on bounded domains a classical result is the convergence rate of order $\mathcal{O}(r^{-k}(\log r)^{(d-1)(r+1)})$, cf. Gerstner and Griebel (1998) and Novak and Ritter (1996). In the case of sparse grid Gauss-Hermite quadrature on the unbounded domain $\mathbb{R}$ it holds, cf. Zhang et al. (2013),

[TABLE]

which asymptotically behaves like $\mathcal{O}(r^{-k/2+\varepsilon})$. Here, as before $C(k,g)$ depends on $k$ and $g$ and the integrand $g$ must have continuous partial derivatives up to order $k$ which fulfill

[TABLE]

i.e. the mixed derivative of order $k$ is continuous and does not grow too fast as ${\mathbf{v}}$ approaches infinity. For integrands that are infinitely often differentiable, $k\in\mathbb{N}$ can be chosen as large as desired, but then also the constant $C(k,g)$ might become arbitrary large.

5 Examples

Example Ia

In the logit model, one needs to compute integrals of the form

[TABLE]

Since the parameter vector ${\boldsymbol{\theta}}=(\theta_{1},\theta_{2})=(\mu,\sigma)$ shall only contribute to $\varphi$ in (4.1) we reformulate the integral (5.1) using the change of variable $v=(x-\mu)/\sigma$ to obtain

[TABLE]

Hence, in the notation of (4.1) we have with ${\boldsymbol{\theta}}=(\theta_{1},\theta_{2}):=(\mu,\sigma)$

[TABLE]

In order to determine the convergence rate of Gauss-Hermite quadrature, we have to analyze the function $v\mapsto\varphi(v,z_{i},{\boldsymbol{\theta}})$, as well as its derivatives with respect to ${\boldsymbol{\theta}}$, i.e. the functions $v\mapsto D_{\boldsymbol{\alpha}}^{({\boldsymbol{\theta}})}\varphi(v,z_{i},{\boldsymbol{\theta}})$ for $|{\boldsymbol{\alpha}}|_{1}\leq 2$. To this end, we use Lemma 20 from the Appendix, which ensures that (4.10) is satisfied for all $k\in\mathbb{N}$, i.e.

[TABLE]

Therefore, we can expect arbitrary large algebraic convergence of order $\mathcal{O}(r^{-k})$.

Example Ib

In the multivariate random coefficients logit model one needs to compute integrals of the form

[TABLE]

We use the variable transformation

[TABLE]

where ${\boldsymbol{C}}$ denotes the Cholesky factorization of $\Sigma$, i.e. $\Sigma={\boldsymbol{C}}{\boldsymbol{C}}^{t}$. Then, for any $f:\mathbb{R}^{d}\to\mathbb{R}$ we obtain

[TABLE]

Hence, in the notation of (4.3) we have with ${\boldsymbol{\theta}}=(\theta_{1},\theta_{2},\ldots,\theta_{p}):=({\boldsymbol{\mu}},{\boldsymbol{C}})$, where $p=d+d(d+1)/2$, that

[TABLE]

Again, we can use Lemma 20 from the Appendix to ensure that (4.13) is fulfilled for all $k\in\mathbb{N}$, i.e.

[TABLE]

Therefore, for any $k\in\mathbb{N}$ we can (asymptotically) expect algebraic convergence of order $\mathcal{O}(r^{-k/2}(\log n)^{(d-1)(k/2+1)})$, which asymptotically behaves like $\mathcal{O}(r^{-k/2+\varepsilon})$.

Example 2

In the Butler-Moffitt model one has to solve univariate integrals of the form

[TABLE]

where $\Phi:\mathbb{R}\to(0,1)$ denotes the cumulative distribution functions of the standard normal distribution.

Hence, in the notation of (4.1) we have with ${\boldsymbol{\theta}}=(\theta_{1},\theta_{2}):=(\sigma,\beta)$

[TABLE]

In order to determine the convergence rate of Gauss-Hermite quadrature, we have to analyze the function $v\mapsto\varphi(v,{\mathbf{z}}_{i},{\boldsymbol{\theta}})$, as well as its derivatives with respect to ${\boldsymbol{\theta}}$, i.e. the functions $v\mapsto D_{\boldsymbol{\alpha}}^{({\boldsymbol{\theta}})}\varphi(v,{\mathbf{z}}_{i},{\boldsymbol{\theta}})$ for $|{\boldsymbol{\alpha}}|_{1}\leq 2$. To this end, we note that $\Phi$ is bounded and all of its derivatives are also bounded. Therefore, by the general Leibnitz / product rule, (4.10) is fulfilled for all $k\in\mathbb{N}$, i.e.

[TABLE]

Thus, we can expect arbitrary large algebraic convergence of order $\mathcal{O}(r^{-k})$.

Remark 14**.**

Note that an arbitrary large algebraic rate of convergence usually amounts to an (sub-)exponential convergence of type $\mathcal{O}(\exp(-\alpha r^{\beta}))$. However, to determine the exact rate of decay, i.e. $\alpha$ and $\beta$ or even to check the required assumptions, usually requires advanced techniques from complex analysis, cf. e.g. Davis and Rabinowitz (2007) or Boyd (1984), which are beyond the scope of this article. Therefore, at this point, we will stay with algebraic convergence rates for any $k\in\mathbb{N}$ and their much simpler to check prerequisites.

6 Pre-Asymptotics

So far we have discussed the asymptotic behavior of the maximum approximated likelihood estimator where $n$ and $r$ approach infinity. In the following simulation study we focus on the illustration of the theoretical findings as well as practical relevant insights into the pre-asymptotic behavior. First we demonstrate the relative performance of common approximation methods and link the result to our theory. Second we consider the consequences caused by the approximation error for the estimation accuracy.

Relative Performance of Approximation Methods

To ensure consistency in the context of Corollary 10 an approximation method is required for which the approximation error disappears if $r\rightarrow\infty$. For all reasonable approaches, including Monte-Carlo methods, quasi Monte-Carlo and numerical quadrature this assumption is fulfilled. However, for a finite $r$ the approximation accuracy can differ dramatically.

For our illustration we analyze a simplified random coefficient regression model of the form

[TABLE]

The likelihood contribution is

[TABLE]

where $g(\cdot)$ denotes the standard normal density function with $g(x)=\frac{1}{\sqrt{2\pi}}e^{-x^{2}/2}$. The likelihood contribution is evaluated at $\theta=0$ using

[TABLE]

where $v_{j,r}$ denotes the draws or nodes regarding the standard normal distribution and $w_{j,r}$ the corresponding weights. For every run in the experiment we randomly generate $x_{i}\sim N(0,1)$, $\beta_{i}\sim N(0,1)$ and $\varepsilon_{i}\sim N(0,1)$. Based on the generated data we compute the approximated likelihood contribution.

We apply three approximation methods, i.e. (1) ordinary Monte-Carlo sampling using pseudo-random draws, (2) quasi Monte-Carlo sampling proposed by Halton (1964) and (3) Gauss-Hermite quadrature333We also tested quasi Monte-Carlo draws following Sobol (1967) and, representative for variance reduction techniques, we use Modified Latin Hypercube Sampling (MLHS) (see Hess, Train and Polak, 2006) as well as Gauss-Legendre quadrature and the midpoint-rule. The general result is that Monte-Carlo performs worse, quasi Monte-Carlo and MLHS perform roughly the same and Gaussian-Quadrature performs best.. As reference $f_{\text{ref}}\simeq f$ we compute the approximation $\tilde{f}(z_{i},\theta=0)$ with a very fine Gauss-Hermite quadrature using 100 grid-points.

For an increasing sequence of $r$ we repeat the experiment $m=5.000$ times and compute for each repetition the approximation error ($\tilde{f}-f_{\text{ref}}$). We aggregate the results in Figure 6.1 using the "maximum absolute error" (left panel) and the RMSE$=\sqrt{\frac{1}{m}\sum_{i=1}^{m}(\tilde{f}-f_{\text{ref}})^{2}}$ (right panel). The maximum absolute error complies with the worst case error in assumption (iii) of Theorem 7. The RMSE is a more common measure of convergence used here for the comparison of the approximation methods.

As shown in Figure 6.1 the three methods are converging in the worst case error and in the average error, respectively. But the methods differ in their convergence rates. The function $f(\cdot)$ used here is smooth and therefore we observe the corresponding theoretical error bounds of the different employed methods (see Table 1). As expected, Gauss-quadrature outperforms the other methods and obtains even an exponential convergence rate. Thus, to reach at least the same accuracy as Monte-Carlo with 16.384 draws, just 128 Halton-draws and even only 16 nodes for the Gauss-Hermite quadrature are needed.

As stated in section 4.2 the convergence behavior and the relative performance is related to the smoothness of the approximated function. In general there is no absolute dominant approximation method but, based on the properties of the likelihood function, one method could perform better than another, provided that more smoothness is present, compare Table 1.

To illustrate this result we also analyze the Accept-Reject-Sampling (ARS) solution for the standard-normal-cdf (see equation (6.2)). This model has a discontinuity at $z=x$ and is therefore not smooth in terms of section 4.2. We consider

[TABLE]

where $1(\cdot)$ denotes the indicator function, which is one if the expression in brackets is true and zero otherwise. This function is approximated using

[TABLE]

For our simulation study we randomly draw 5.000 z-values from a standard-normal-distribution for each accuracy level $r$. We apply Monte-Carlo, Halton draws and Gauss-Legendre quadrature. As reference solution we take the error-function representation $f_{\text{ref}}=\frac{1}{2}\left[1+\operatorname{erf}\left(\frac{z}{\sqrt{2}}\right)\right]$. The results are summarized in Figure 6.2.

The result for the non-smooth ARS confirms that the relative performance depends strongly on the smoothness of the approximated function. Here Monte-Carlo achieves the same convergence rate as before. The Halton draws and Gauss-Legendre achieve a higher rate, but Gauss-quadrature does not outperform the other two methods any more.

As stated in section 3.2 the convergence rate of an approximation method translates into the required link-function $R(n)$ to achieve asymptotic normality. Figure 6.3 shows the results for the previous model $f(\cdot)$ from (6.1) were the scaled approximation error $\mathcal{E}(r)$ (see equation (3.7)) is plotted against the sample size in terms of Theorem 8 (ii). For Theorem 8 we require a link function which ensures $\operatorname*{plim}_{n\rightarrow\infty}\sqrt{n}\mathcal{E}(R(n))=0$. For the constant link-function (upper left panel) this condition is not met by any approximation method. For the logarithmic link-function (upper right panel) and for the $\sqrt{n}$ link-function (lower left panel) only the more efficient Gauss-Hermite or the Gauss-Hermite and Halton draws, respectively, meet the condition. And for the computationally most costly link function $n$ (lower right panel) all methods meet the condition, even the Monte-Carlo method (blue line), albeit with a very flat slope.

Consequences for the practical application of MALE

The overall estimation error ($\hat{\theta}_{MAL}-\theta_{0}$) in MAL has two components (compare also equation (2.4)). First, the sampling error and, second, the approximation error. In most practical applications, the sample size is typically fixed and there is only a small variation allowed in the approximation accuracy. Thus, we now run a simulation were we generate $m=2000$ data sets for our random coefficient model with sample size fixed at $n=50$ and $n=5000$, respectively. We estimate the parameter $\bar{\beta}$ using the likelihood function (6.1) from above. To illustrate the convergence behavior we compute the $RMSE=\sqrt{\frac{1}{m}\sum_{i=1}^{m}(\hat{\theta}_{MALE}-\theta_{0})^{2}}$ based on the estimated parameter $\theta=\bar{\beta}$. Figure 6.4 shows the result.

If there would be no approximation error only the sample error would cause variation in the result and would not be affected by $r$. This sample-error expressed in the $RMSE$ is given by the dashed horizontal line in Figure 6.4. For increasing sample size this error decreases (compare left and right panel). Indeed, the estimator’s variance needed for standard hypothesis testing is just of this type.

For finite $r$ there is also an approximation error whose size depends on the approximation method and $r$. This error component increases the variance for the estimator and therefore potentially affects the interpretation of hypothesis tests. To correct for the approximation error in the inference one could apply corrections as proposed in Kristensen and Salanié (2017). So, also for practical applications, it is recommended to choose the most efficient method to minimize this error-component for a given computational cost. Moreover, Figure 6.4 shows an increasing relevance of the approximation error due a rising number of observations. This leads to the conclusion that, for a large dataset, an appropriate approximation is of higher interest than for a small dataset.

7 Summary and conclusions

This paper discusses maximum approximated likelihood (MAL) estimators that generalize maximum simulated likelihood (MSL) estimators. A major advantage of MSL is that the underlying Monte Carlo simulation techniques provide favorable asymptotic properties under very general conditions. This not only makes it a versatile tool for the practitioner, it also simplifies the theoretical analysis. However, it has been frequently found that the computational costs required to achieve a sufficient approximation quality can be burdensome or infeasible. So it has become common practice to use more accurate numerical approximation algorithms such as quasi–Monte Carlo simulation, Gaussian quadrature, and integration on sparse grids. This paper contributes to the theoretical underpinning of these approaches. It establishes sets of conditions on the model, the approximation approach and their interactions that ensure the consistency and the asymptotic efficiency of general MAL estimators. We also provide discussions of specific algorithms and models and show their behavior both asymptotically and in a simulation analysis.

Overall, numerical approximation methods have stronger requirements on properties of the models such as the smoothness of the likelihood contributions than MSL. Given that these conditions are met, not only their finite sample properties, but also their asymptotic behavior can be superior. This is manifested mainly in the fact that the speed with which the computational burden has to be asymptotically increased (as $n\rightarrow\infty$) can be dramatically decreased relative to MSL.

The verification of the general conditions provided in this paper for specific models and approximation methods can be a nontrivial task. We provide examples to demonstrate the approach. But we expect that more work discussing the specifics for different classes of models and algorithms would be useful for guiding the practitioner.

Appendix A Appendix: Technical results

In Section 3, we have worked with distances of the log likelihood function from its quadrature-approximation as well as their derivatives. In order to deal with the logarithm in the log likelihood contribution, we now present a general result that relates the distance of a function or its derivatives up to order $2$ to the respective distance between the logarithmized functions and its derivatives.

To this end, we will again use the notation $D_{{\boldsymbol{\alpha}}}^{({\boldsymbol{\theta}})}g({\boldsymbol{\theta}})$ to denote the partial derivatives with respect to ${\boldsymbol{\theta}}$, cf. (4.5). Moreover, we use $\|\nabla_{\boldsymbol{\theta}}g\|$ to denote the $\ell_{2}$-(vector-)norm of the gradient of $g$ and $\|\nabla_{{\boldsymbol{\theta}},{\boldsymbol{\theta}}}g\|$ to denote the $\ell_{2}$-(matrix-)norm of the Hessian of $g$.

Moreover, we will use the $L^{\infty}(X)$ space which contains all functions that are bounded on some domain $X$, i.e.

[TABLE]

Theorem 15**.**

Consider real-valued functions $g$ and $h$ defined on a compact set $\Theta\subset\mathbb{R}^{p}$. Assume that (i) the image of $g$ and $h$ is bounded away from zero, i.e. for all ${\boldsymbol{\theta}}\in\Theta$ it holds $g({\boldsymbol{\theta}}),h({\boldsymbol{\theta}})\in J=[\delta,D]$ with $0<\delta<D<\infty$; (ii) the partial derivatives of $g$ and $h$ up to order $k\in\mathbb{N}_{0}$ exist and are bounded in $L^{\infty}(\Theta)$. Moreover, consider a function $\psi:J\to\mathbb{R}$, whose first $k+1$ derivatives exist and are bounded in $L^{\infty}(J)$. Then, it holds for

(i) $k=0$, i.e. $g,h\in L^{\infty}(\Theta)$ and $\psi,\psi^{\prime}\in L^{\infty}(J)$, that

[TABLE]

(ii) $k=1$, i.e. $g,h,\|\nabla_{\boldsymbol{\theta}}g\|,\|\nabla_{\boldsymbol{\theta}}h\|\in L^{\infty}(\Theta)$ and $\psi,\psi^{\prime},\psi^{\prime\prime}\in L^{\infty}(J)$, that

[TABLE]

(iii) $k=2$, i.e. $g,h,\|\nabla_{\boldsymbol{\theta}}g\|,\|\nabla_{\boldsymbol{\theta}}h\|,\|\nabla_{{\boldsymbol{\theta}},{\boldsymbol{\theta}}}g\|,\|\nabla_{{\boldsymbol{\theta}},{\boldsymbol{\theta}}}h\|\in L^{\infty}(\Theta)$ and $\psi,\psi^{\prime},\psi^{\prime\prime},\psi^{\prime\prime\prime}\in L^{\infty}(J)$, that

[TABLE]

The theorem is proven by the following three Lemmas.

Lemma 16**.**

For a differentiable function $\psi:J\to\mathbb{R}$, with $[a,b]\subseteq J\subset\mathbb{R}$, it holds that

[TABLE]

Proof.

The claim follows immediately from the mean-value theorem, i.e. there exists $\xi\in[a,b]$ such that it holds

[TABLE]

∎

Lemma 16 immediately implies

[TABLE]

which implies Theorem 15 for the case $k=0$, i.e. (A.1).

Next, we deal with the gradient.

Lemma 17**.**

For a twice differentiable function $\psi:J\to\mathbb{R}$ with $\psi,\psi^{\prime},\psi^{\prime\prime}\in L^{\infty}(J)$ and $J$-valued functions $g,h:\Theta\to J$ with bounded derivatives of first order it holds that

[TABLE]

Proof.

Using the chain rule for the gradient $\nabla_{\boldsymbol{\theta}}\psi\circ g({\boldsymbol{\theta}})=\psi^{\prime}(g({\boldsymbol{\theta}}))\cdot\nabla g({\boldsymbol{\theta}})$ we compute

[TABLE]

By $|\psi^{\prime}\circ g({\boldsymbol{\theta}})|\leq\|\psi^{\prime}\|_{L_{\infty}(J)}$ and by applying Lemma 16 we obtain

[TABLE]

Taking the $\sup$ over all ${\boldsymbol{\theta}}\in\Theta$ we arrive at the desired result, which is (A.2), i.e. (ii) in Theorem 15.

∎

Next, we deal with the Hessian matrix $\nabla_{{\boldsymbol{\theta}}{\boldsymbol{\theta}}}f({\boldsymbol{\theta}})=(\partial_{\theta_{k}}\partial_{\theta_{j}}f({\boldsymbol{\theta}}))_{k,j=1}^{p}$ .

Lemma 18**.**

For a $3$-times differentiable function $\psi:J\to\mathbb{R}$ with $\psi,\psi^{\prime},\psi^{\prime\prime},\psi^{\prime\prime\prime}\in L^{\infty}(J)$ and $J$-valued functions $g,h:\Theta\to J$ with bounded derivatives up to second order it holds that

[TABLE]

Proof.

First we note that it holds by the multivariate chain-rule that

[TABLE]

where we abbreviated $p_{\boldsymbol{\theta}}(g):=(\nabla_{\boldsymbol{\theta}}g({\boldsymbol{\theta}}))\cdot(\nabla_{\boldsymbol{\theta}}g({\boldsymbol{\theta}}))^{t}$. Here, $p_{\boldsymbol{\theta}}(g)\in\mathbb{R}^{p\times p}$ is a rank-1 matrix that consists of the outer product of the gradient of $g$ at ${\boldsymbol{\theta}}$ with itself. For its norm it holds that $\|p_{\boldsymbol{\theta}}(g)\|=\|\nabla_{\boldsymbol{\theta}}g({\boldsymbol{\theta}})\|^{2}$. Moreover, $\nabla_{{\boldsymbol{\theta}}{\boldsymbol{\theta}}}g({\boldsymbol{\theta}})$ denotes the Hessian matrix of $g$ at ${\boldsymbol{\theta}}$.

Now, we proceed analogously to the proof of the preceeding Lemma and compute

[TABLE]

Next, we will take care of the last summand in (A.6) and derive

[TABLE]

Inserting this back into (A.6) we obtain

[TABLE]

It remains to bound the term \big{\|}p_{\boldsymbol{\theta}}(g)-p_{\boldsymbol{\theta}}(h)\big{\|} in (A.6). To this end, we derive for vectors ${\mathbf{v}},{\mathbf{w}}\in\mathbb{R}^{p}$

[TABLE]

Now, using (A.8) with ${\mathbf{v}}=\nabla_{\boldsymbol{\theta}}g({\boldsymbol{\theta}})$ and ${\mathbf{w}}=\nabla_{\boldsymbol{\theta}}h({\boldsymbol{\theta}})$ we obtain

[TABLE]

Hence, inserting (A.9) into (A.7) we arrive at

[TABLE]

which concludes the proof.

∎

For the special case $\psi(x)=\log(x)$ we obtain the following result.

Corollary 19**.**

Consider real-valued functions $g,h$ on some compact domain $\Theta\subset\mathbb{R}^{p}$. Assume that the image of $g$ and $h$ is bounded away from zero, i.e. for all ${\boldsymbol{\theta}}\in\Theta$ it holds $g({\boldsymbol{\theta}}),h({\boldsymbol{\theta}})\in J=[\delta,D]$ with $0<\delta<1\leq D<\infty$. Then

(i) for $g,h\in L^{\infty}(\Theta)$ it holds

[TABLE]

(ii) for differentiable $g,h$ with bounded first derivatives it holds

[TABLE]

where

[TABLE]

(iii) for $2$-times differentiable $g,h$ with bounded derivatives up to second order it holds

[TABLE]

where

[TABLE]

Proof.

First note that

[TABLE]

Then, for (i) and (ii) we use that $\|\psi^{(j)}\|_{L_{\infty}(J)}=\frac{(j-1)!}{\delta^{j}}$ with $j\in\mathbb{N}$ and $J=[\delta,\infty)$. For (iii) we additionally use the inequality $\max(x,x^{2})\leq 1+x^{2}$. ∎

We finally prove the following result that was employed in the analysis of the logit model in Section 5.

Lemma 20**.**

Let ${\boldsymbol{\theta}},{\mathbf{v}}\in\mathbb{R}^{n}$ and

[TABLE]

Then it holds for all ${\boldsymbol{\alpha}},{\boldsymbol{\beta}}\in\mathbb{N}_{0}^{n}$ that

1. (i)

there exists a constant $0<c({\boldsymbol{\alpha}},{\boldsymbol{\beta}})<\infty$ such that

[TABLE] 2. (ii)

there exists a constant $0<\tilde{c}({\boldsymbol{\alpha}},{\boldsymbol{\beta}})<\infty$ such that

[TABLE]

Proof.

First, we note that (ii) follows from (i), because the function $t\mapsto e^{t^{2}/2}/\sqrt{1+t^{2}}$ grows faster than any polynomial.

In order to prove (i), we write $f({\boldsymbol{\theta}},{\mathbf{v}}):=g({\boldsymbol{\theta}}\cdot{\mathbf{v}})$, where $g(t)=(1+e^{-t})^{-1}$. Note that $g$ and hence $g^{p}$ are bounded on $\mathbb{R}$ for all $p\in\mathbb{N}$. Next, we use that

[TABLE]

and likewise also $\frac{\partial}{\partial v_{i}}g({\boldsymbol{\theta}}\cdot{\mathbf{v}})=\theta_{i}(g({\boldsymbol{\theta}}\cdot{\mathbf{v}})-g^{2}({\boldsymbol{\theta}}\cdot{\mathbf{v}}))$. This proves (i) for $|{\boldsymbol{\alpha}}|=1$ and $|{\boldsymbol{\beta}}|=1$. The case of higher order derivatives follows by induction by showing that $D^{({\boldsymbol{\theta}})}_{\boldsymbol{\alpha}}D^{({\mathbf{v}})}_{\boldsymbol{\beta}}g({\boldsymbol{\theta}}\cdot{\mathbf{v}})$ always has a representation as $\sum_{j}c_{j}{\mathbf{v}}^{{\mathbf{a}}_{j}}{\boldsymbol{\theta}}^{{\mathbf{b}}_{j}}g^{p_{j}}({\boldsymbol{\theta}}\cdot{\mathbf{v}})$ with ${\mathbf{a}}_{j}\leq{\boldsymbol{\alpha}}$, ${\mathbf{b}}_{j}\leq{\boldsymbol{\alpha}}$ and $p_{j}\leq|{\boldsymbol{\alpha}}|+|{\boldsymbol{\beta}}|+1$. ∎