Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
84,543 characters · 17 sections · 54 citation commands
Maximum Approximated Likelihood Estimation
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
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
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 Butler1982, McFaddenTrain2000, Train2009, or general limited dependent variable models Hajivassiliou1994. Then, \(f\) can be approximated with an \(r\)-point quadrature rule, i.e.
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
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 (ref) 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 Hajivassiliou1994, GourierouxMonfort1996 and Train2009.
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\).\footnote{"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 Hajivassiliou1994. A recent discussion about the difference between both methodologies can also be found in kristensen2017.}
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 Butler1982 is an obvious choice. For multivariate problems, quasi-Monte-Carlo methods like Halton draws Bhat2001 or other deterministic approaches like integration on sparse grids Heiss2008 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 Butler1982. Other papers discuss the well-studied simulation-based estimators before arguing that their deterministic approaches tend to work better in practice, e.g. in Bhat2001 or in sandor2004. This might be due to the fact that the theoretical properties and prerequisites of deterministic approximation schemes are not yet very well understood. Ackerberg2009 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 (ref) by (ref) 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 (ref) we discuss conditions that imply the consistency and asymptotic normality of general extremum estimators. In Section (ref) 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 (ref), we deal with specific integration algorithms, like Gaussian quadrature, quasi--Monte Carlo and sparse grids. In Section (ref), 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 (ref).
As the most general framework, we consider extremum or M-estimators as discussed in Newey1994, NM hereafter. An M-estimator \(\hat{{\boldsymbol{\theta}}}_{M}\) maximizes (or minimizes) some objective function \(Q_n({\boldsymbol{\theta}})\), i.e.
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
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:
We abstract from any misspecifications and other problems that could violate the conditions of Lemma (ref) 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 (ref) 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 (ref) 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}})\).
Note that assumption (ref)(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\).
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}\).
Arguably the most important condition for deriving the asymptotic distribution is assumption (ref)(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\).
In the last section we discussed general extremum estimators with approximated objective functions. These results are similar to those of Ackerberg2009 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.
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
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 (ref) and is maximized by the true parameter \({\boldsymbol{\theta}}_0\).
Newey1994 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.
Next, we recall the following result on the asymptotic normality of the ML estimator.
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 (ref) and (ref) 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
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
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 (ref) and (ref) 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 (ref) is the uniform convergence in probability of the approximation \(\tilde{L}_{n}({\boldsymbol{\theta}})\) to the true, but intractable \(L_{n}({\boldsymbol{\theta}})\).
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
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.
In Theorem (ref) it is required by condition (ii) that \(\operatorname*{plim}_{n \to \infty } \sqrt{n}\,\mathcal{E}(R(n)) = 0\), i.e. the approximation error (ref) 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}\).
Now, we give general results on how to choose \(R(n)\) depending on the convergence rates of the approximation algorithm for \(\tilde{f}_r\).
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
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 Butler1982 and random coefficients models like the mixed logit model, see for example McFaddenTrain2000.
The need for approximation arises because the integral in (ref) often cannot be computed in closed form. To this end, several estimators have been proposed including the method of simulated moments (MSM; McFadden1989) and the method of maximum simulated scores (MSS; HajivassiliouMcFadden1998). 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 Bhat2001, Gaussian quadrature in Butler1982, or quadrature on sparse grids in Heiss2008.
To cover all the mentioned approaches, we consider a general approximation\footnote{Note 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.}
to the true integral of (ref). 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. rabinowitz, as well as sparse grids GerstnerGriebel: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
Their associated approximations \(\tilde{f}_r({\mathbf{z}} , {\boldsymbol{\theta}})\) stem from some cubature rule \(U_r\) of the form
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 (ref) and (ref) hold. To this end, the following notation for the partial derivative of a sufficiently smooth function $g$ will be helpful.
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 (ref) and (ref), we resort to
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 2 p \bar{\mathcal{E}}_1(r)\) and hence it is sufficient (and convenient) to work with the quantity defined in (ref).
Now, \(\lim_{r \to \infty} \bar{\mathcal{E}}_k(r) = 0\) basically implies that the quadrature rules \((U_r)_{r \in \mathbb{N}}\) in (ref) converge for \(\varphi(\cdot, {\mathbf{z}} , {\boldsymbol{\theta}})\) and all partial derivatives with respect to \({\boldsymbol{\theta}}\) up to order \(k\).
First, we deal with consistency, which follows as a simple application of Theorem (ref) to the specific setting where \(f\) has the form (ref) and \(\tilde{f}_r\) has the form (ref).
Condition (ref)(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 (ref) and \(\tilde{f}_r\) has the form (ref). The critical part of the properties of \(\hat{\boldsymbol{\theta}}_{MAL}\) is the worst-case approximation error, i.e. (ref) for \(k=1\), which not only has to go to zero, but also has to vanish at a certain speed.
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 (ref) 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 (ref). Several concrete examples are provided in Table (ref).
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 (ref) 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. \[ \sup_{{\mathbf{z}} , {\boldsymbol{\theta}}} \int_\Omega \left( D^{({\boldsymbol{\theta}})}_{\boldsymbol{\alpha}} \varphi({\mathbf{v}} , {\mathbf{z}} , {\boldsymbol{\theta}}) \right)^2 \, \omega({\mathbf{v}} ) \mathrm{d} {\mathbf{v}} < \infty \quad \text{ for all } |{\boldsymbol{\alpha}} |_1 \leq 2. \] Then, with high probability and also in expectation, \(\bar{\mathcal{E}}_1 \leq c r^{-1/2}\) and \(\bar{\mathcal{E}}_1\) converges to zero at a rate of \(\mathcal{O}(r^{-1/2})\). It follows by Corollary (ref) and Theorem (ref) 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\).
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. SloanJoe:1994 or Sobol points, cf. Sobol:1967. A rather recent development are higher-order QMC points, cf. DiPi10,Hinrichs.Markhasin.Oettershagen.ea: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 DiPi10,Niederreiter:1992,owen2005multidimensional. 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.
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 (ref) 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 (ref) and Theorem (ref) 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.
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
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. rabinowitz. 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.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. Patterson1967, Clenshaw-Curtis quadrature, cf. Imhof1963 or Leja points, cf. jantsch2016lebesgue,Griebel.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 (ref) 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\).
\paragraph{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. rabinowitz. However, if the integrand is only \(k < \infty\) times differentiable, the convergence rate deteriorates to the order \(\mathcal{O}(r^{-k})\).
\paragraph{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. Smith2018 or Mastroianni94. To be more precise, if there exists a constant \(c > 0\) such that
then it holds
where the constant \( C(k,g) > 0\) may depend on \(k\) and the integrand \(g\), but not on \(r\). Basically, condition (ref) 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 boyd1984. Finally, for certain classes of entire functions, even exponential convergence \(\mathcal{O}(e^{-\alpha r})\) is possible, cf. Kuo2011 or Irrgeher2015.
For multivariate integration problems it is possible to use sparse grids, cf. GerstnerGriebel:1998,GerstnerGriebel:2003,Griebel.Oettershagen:2016,Novak1996 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. GerstnerGriebel:1998 and Novak1996. In the case of sparse grid Gauss-Hermite quadrature on the unbounded domain \(\mathbb{R}\) it holds, cf. Zhang2013,
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
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.
\paragraph{Example Ia} In the logit model, one needs to compute integrals of the form
Since the parameter vector \({\boldsymbol{\theta}} = (\theta_1,\theta_2) = (\mu, \sigma)\) shall only contribute to \(\varphi\) in (ref) we reformulate the integral (ref) using the change of variable \(v = (x-\mu)/\sigma\) to obtain
Hence, in the notation of (ref) we have with \({\boldsymbol{\theta}} = (\theta_1,\theta_2) := (\mu, \sigma)\)
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 (ref) from the Appendix, which ensures that (ref) is satisfied for all \(k \in \mathbb{N}\), i.e.
Therefore, we can expect arbitrary large algebraic convergence of order \(\mathcal{O}(r^{-k})\).
\paragraph{Example Ib} In the multivariate random coefficients logit model one needs to compute integrals of the form
We use the variable transformation \[ {\mathbf{u}} = {\boldsymbol{C}}^{-1}({\mathbf{v}} - {\boldsymbol{\mu}} ), \] 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
Hence, in the notation of (ref) we have with \({\boldsymbol{\theta}} = (\theta_1,\theta_2, \ldots, \theta_p) := ({\boldsymbol{\mu}} , {\boldsymbol{C}})\), where \(p = d + d(d+1)/2\), that
Again, we can use Lemma (ref) from the Appendix to ensure that (ref) is fulfilled for all \(k \in \mathbb{N}\), i.e.
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})\).
\paragraph{Example 2} In the Butler-Moffitt model one has to solve univariate integrals of the form \[ f({\mathbf{z}} _i, {\boldsymbol{\theta}}) = \int_\mathbb{R} \prod_{t=1}^T \Phi( z_{i,t} \beta + \sigma v) \frac{e^{-\frac{v^2}{2}}}{\sqrt{2 \pi}} \, \mathrm{d} v , \] where \(\Phi: \mathbb{R} \to (0,1)\) denotes the cumulative distribution functions of the standard normal distribution.
Hence, in the notation of (ref) we have with \({\boldsymbol{\theta}} = (\theta_1,\theta_2) := (\sigma, \beta)\)
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, (ref) is fulfilled for all \(k \in \mathbb{N}\), i.e.
Thus, we can expect arbitrary large algebraic convergence of order \(\mathcal{O}(r^{-k})\).
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.
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
\[ y_{i}=x_{i}\beta_{i}+\varepsilon_{i}\quad\text{with \ensuremath{\varepsilon_{i}\sim N(0,1)\quad \beta_{i}\sim N(\bar{\beta},1)}}. \]
The likelihood contribution is
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
\[ \tilde{f}\left(z_i, \theta \right) =\sum_{j=1}^{r}w_{j,r}g(y-xv_{j,r}), \] 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 quadrature\footnote{We also tested quasi Monte-Carlo draws following Sobol:1967 and, representative for variance reduction techniques, we use Modified Latin Hypercube Sampling (MLHS) \citep*[see][]{Hess2006} 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 (ref) 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 (ref). The RMSE is a more common measure of convergence used here for the comparison of the approximation methods.
As shown in Figure (ref) 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 (ref)). 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 (ref) 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 ((ref))). This model has a discontinuity at $z=x$ and is therefore not smooth in terms of section (ref). We consider
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
\[ \tilde{f}\left(z\right) =\sum_{j=1}^{r}w_{j,r}1(v_{j,r} \leq z) ). \]
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 (ref).
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 (ref) the convergence rate of an approximation method translates into the required link-function $R(n)$ to achieve asymptotic normality. Figure (ref) shows the results for the previous model $f(\cdot)$ from ((ref)) were the scaled approximation error $\mathcal{E}(r)$ (see equation ((ref))) 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.
The overall estimation error ($\hat{\theta}_{MAL} - \theta_0$) in MAL has two components (compare also equation ((ref))). 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 ((ref)) 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 (ref) 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 (ref). 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 kristensen2017. 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 (ref) 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.
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.