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.
87,552 characters · 17 sections · 83 citation commands
The ABC of Simulation Estimation with Auxiliary Statistics
JEL Classification: C22, C23.\\
Keywords: Indirect Inference, Simulated Method of Moments, Efficient Method of Moments, Laplace Type Estimator.
\baselineskip=18.0pt \thispagestyle{empty} \setcounter{page}{0}
As knowledge accumulates, scientists and social scientists incorporate more and more features into their models to have a better representation of the data. The increased model complexity comes at a cost; the conventional approach of estimating a model by writing down its likelihood function is often not possible. Different disciplines have developed different ways of handling models with an intractable likelihood. An approach popular amongst evolutionary biologists, geneticists, ecologists, psychologists and statisticians is Approximate Bayesian Computation (ABC). This work is largely unknown to economists who mostly estimate complex models using frequentist methods that we generically refer to as the method of Simulated Minimum Distance (SMD), and which include such estimators as Simulated Method of Moments, Indirect Inference, or Efficient Methods of Moments.\footnote{ Indirect Inference is due to gmr, the Simulated Method of moments is due to duffie-singleton, and the Efficient Method of Moments is due to gallant-tauchen-emm.}
The ABC and SMD share the same goal of estimating parameters $\theta$ using auxiliary statistics $\hat\psi$ that are informative about the data. An SMD estimator minimizes the $\mathsf L_2$ distance between $\hat\psi$ and an average of the auxiliary statistics simulated under $\theta$, and this distance can be made as close to zero as machine precision permits. An ABC estimator evaluates the distance between $\hat\psi$ and the auxiliary statistics simulated for each $\theta$ drawn from a proposal distribution. The posterior mean is then a weighted average of the draws that satisfy a distance threshold of $\delta>0$. There are many ABC algorithms, each differing according to the choice of the distance metric, the weights, and sampling scheme. But the algorithms can only approximate the desired posterior distribution because $\delta$ cannot be zero, or even too close to zero, in practice.
While both SMD and ABC use simulations to match $\psi(\theta)$ to $\hat\psi$ (hence likelihood-free), the relation between them is not well understood beyond the fact that they are asymptotically equivalent under some high level conditions. To make progress, we focus on the MCMC-ABC algorithm due to mmpt-03. The algorithm applies uniform weights to those $\theta$ satisfying $\|\hat\psi-\psi(\theta)\|\le \delta$ and zero otherwise. Our main insight is that this $\delta$ can be made very close to zero if we combine optimization with Bayesian computations. In particular, the desired ABC posterior distribution can be targeted using a `Reverse Sampler' (or RS for short) that applies importance weights to a sequence of SMD solutions. Hence, seen from the perspective of the RS, the ideal MCMC-ABC estimate with $\delta=0$ is a weighted average of SMD modes. This offers a useful contrast with the SMD estimate, which is the mode of the average deviations between the model and the data. We then use stochastic expansions to study sources of variations in the two estimators in the case of exact identification. The differences are illustrated using simple analytical examples as well as simulations of the dynamic panel model.
Optimization of models with a non-smooth objective function is challenging, even when the model is not complex. The Quasi-Bayes (LT) approach due to chernozhukov-hong use Bayesian computations to approximate the mode of a likelihood-free objective function. Its validity rests on the Laplace (asymptotic normal) approximation of the posterior distribution with the goal of valid asymptotic frequentist inference. The simulation analog of the LT (which we call SLT) further uses simulations to approximate the intractable relation between the model and the data. We show that both the LT and SLT can also be represented as a weighted average of modes with appropriately defined importance weights.
A central theme of our analysis is that the mean computed from many likelihood-free posterior distributions can be seen as a weighted average of solutions to frequentist objective functions. Optimization permits us to turn the focus from computational to analytical aspects of the posterior mean, and to provide a bridge between the seemingly related approaches. Although our optimization-based samplers are not intended to compete with the many ABC algorithms that are available, they can be useful in situations when numerical optimization of the auxiliary model is fast. This aspect is studied in our companion paper jjng-15 in which implementation of the RS in the overidentified case is also considered. The RS is independently proposed in meeds-welling with emphasis on efficient and parallel implementations. Our focus on the analytical properties complements their analysis.
The paper proceeds as follows. After laying out the preliminaries in Section 2, Section 3 presents the general idea behind ABC and introduces an optimization view of the ideal MCMC-ABC. Section 4 considers Quasi-Bayes estimators and interprets them from an optimization perspective. Section 5 uses stochastic expansions to study the properties of the estimators. Section 6 uses analytical examples and simulations to illustrate their differences. Throughout, we focus the discussion on features that distinguish the SMD from ABC which are lesser known to economists.\footnote{ The class of SMD estimators considered are well known in the macro and finance literature and with apologies, many references are omitted. We also do not consider discrete choice models; though the idea is conceptually similar, the implementation requires different analytical tools. smith-palgrave provides a concise overview of these methods. The finite sample properties of the estimators are studied in michaelides-ng. Readers are referred to the original paper concerning the assumptions used.}
As a matter of notation, we use $L(\cdot)$ to denote the likelihood, $p(\cdot)$ to denote posterior densities, $q(\cdot)$ for proposal densities, and $\pi(\cdot) $ to denote prior densities. A `hat' denotes estimators that correspond to the mode and a `bar' is used for estimators that correspond to the posterior mean. We use $(s,S)$ and $(b,B)$ to denote the (specific, total number of) draws in frequentist and Bayesian type analyses respectively. A superscript $s$ denotes a specific draw and a subscript $S$ denotes the average over $S$ draws. For a function $f(\theta)$, we use $f_\theta(\theta_0)$ to denote $\frac{\partial} {\partial \theta}f(\theta)$ evaluated at $\theta_0$, $f_{\theta\theta_j}(\theta_0)$ to denote $\frac{\partial }{\partial \theta_j} f_\theta (\theta)$ evaluated at $\theta_0$ and $f_{\theta,\theta_j,\theta_k}(\theta_0)$ to denote $\frac{\partial^2 }{\partial \theta_j \theta_k} f_\theta (\theta)$ evaluated at $\theta_0$.
Throughout, we assume that the data $\mathbf y=(y_1,\ldots,y_T)^\prime $ are covariance stationary and can be represented by a parametric model with probability measure $\mathcal P_\theta$ where $\theta \in \Theta\subset \mathbbm R^K$. The true value of $\theta$ is denoted by $\theta_0$. Unless otherwise stated, we write $\mathbb E[\cdot]$ for expectations taken under $P_{\theta_0}$ instead of $\mathbb E_{\mathcal P_{\theta_0}}[\cdot]$. If the likelihood $L(\theta)=L(\theta|\mathbf y)$ is tractable, maximizing the log-likelihood $\ell (\theta)=\log L(\theta) $ with respect to $\theta$ gives \[ \hat\theta_{ML}=\argmax_\theta \ell(\theta).\]
Bayesian estimation combines the likelihood with a prior $\pi(\theta)$ to yield the posterior density
For any prior $\pi(\theta)$, it is known that $\hat\theta_ {ML}$ solves $\argmax_\theta \ell(\theta)=\lim_{\lambda\rightarrow \infty} \frac { \int_\Theta \theta\exp(\lambda \ell(\theta)) \pi(\theta )d\theta}{\int_\Theta \exp(\lambda \ell(\theta))\pi(\theta)d\theta}$. That is, the maximum likelihood estimator is a limit of the Bayes estimator using $\lambda\rightarrow\infty$ replications of the data $\mathbf y$.\footnote{See robert-casella, jjp-07.} The parameter $\lambda$ is the cooling temperature in simulated annealing, a stochastic optimizer due to kirkpatrick-gellatt-vecchi for handling problems with multiple modes.
In the case of conjugate problems, the posterior distribution has a parametric form which makes it easy to compute the posterior mean and other quantities of interest. For non-conjugate problems, the method of Monte-Carlo Markov Chain (MCMC) allows sampling from a Markov Chain whose ergodic distribution is the target posterior distribution $p(\theta|\mathbf y)$, and without the need to compute the normalizing constant. We use the Metropolis-Hastings (MH) algorithm in subsequent discussion. In classical Bayesian estimation with proposal density $q(\cdot)$, the acceptance ratio is \[\rho_{BC}(\theta^b,\theta^{b+1}) = \min \Big( \frac{L(\theta^{b+1})\pi(\theta^{b+1})q(\theta^b|\theta^{b+1}) }{L (\theta^b)\pi(\theta^b)q(\theta^{b+1}|\theta^b)},1 \Big).\] When the posterior mode $\hat\theta_{BC}=\argmax_\theta p(\theta|y)$ is difficult to obtain, the posterior mean
is often the reported estimate, where $\theta^b$ are draws from the Markov Chain upon convergence. Under quadratic loss, the posterior mean minimizes the posterior risk $ Q(a)= \int_\Theta |\theta-a|^2 p(\theta|\mathbf y) d\theta$.
The method of generalized method of moments (GMM) is a likelihood-free frequentist estimator developed in hansen-82,hansen-singleton:82. For example, it allows for the estimation of $K$ parameters in a dynamic model without explicitly solving the full model. It is based on a vector of $L\ge K$ moment conditions $ g_t (\theta) $ whose expected value is zero at $\theta=\theta_0$, i.e. $\mathbb E [g_t(\theta_0)]=0$. Let $ \bar g(\theta)=\frac{1}{T}\sum_{t=1}^T g_t(\theta)$ be the sample analog of $\mathbb E[g_t(\theta)]$. The estimator is
where $W$ is a $L\times L$ positive-definite weighting matrix. Most estimators can be put in the GMM framework with suitable choice of $g_t$. For example, when $g_t$ is the score of the likelihood, the maximum likelihood estimator is obtained.
Let $\hat\psi\equiv \hat\psi(\mathbf y(\theta_0))$ be $L$ auxiliary statistics with the property that $\sqrt{T}(\hat\psi -\psi(\theta_0)) \dconv \mathcal N(0,\Sigma)$. It is assumed that the mapping $\psi(\theta)=\lim_{T\rightarrow\infty} \mathbb E[\hat\psi(\theta)]$ is continuously differentiable in $\theta$ and locally injective at $\theta_0$. gmr refer to $\psi(\theta)$ as the {\em binding function} while jiang-turnbull:94 use the term {\em bridge function}. The minimum distance estimator is a GMM estimator which specifies \[ \bar g(\theta)=\hat \psi- \psi(\theta),\] with efficient weighting matrix $W=\Sigma^{-1}$. Classical MD estimation assumes that the binding function $\psi(\theta)$ has a closed form expression so that in the exactly identified case, one can solve for $\theta$ by inverting $\bar g(\theta)$.
Simulation estimation is useful when the asymptotic binding function $\psi(\theta_0)$ is not analytically tractable but can be easily evaluated on simulated data. The first use of this approach in economics appears to be due to smith-93. The simulated analog of MD, which we will call SMD, minimizes the weighted difference between the auxiliary statistics evaluated at the observed and simulated data:
where \[ \bar g_S(\theta) = \hat\psi -\frac{1}{S}\sum_{s=1}^S \hat \psi^s(\mathbf y^s(\theta)), \] $\mathbf y^s(\theta) \equiv \mathbf y^s(\varepsilon^s,\theta)$ are data simulated under $\theta$ with errors $\varepsilon^s$ drawn from an assumed distribution $F_\varepsilon$, and $\hat\psi^s(\theta)\equiv \hat\psi^s (\mathbf y^s (\varepsilon^s,\theta))$ are the auxiliary statistics computed using $\mathbf y^s(\theta)$. Of course, $\bar g_S(\theta)$ is also the average over $S$ deviations between $\hat\psi$ and $\hat\psi^s(\mathbf y^s(\theta))$. To simplify notation, we will write $\mathbf y^s$ and $\hat\psi^s(\theta)$ when the context is clear. As in MD estimation, the auxiliary statistics $\psi(\theta)$ should `smoothly embed' the properties of the data in the terminology of gallant-tauchen-emm. But SMD estimators replace the asymptotic binding function $\psi(\theta_0)=\lim_{T\rightarrow\infty} \mathbb E [\hat\psi (\theta_0)]$ by a finite sample analog using Monte-Carlo simulations. While the SMD is motivated with the estimation of complex models in mind, grt-99 show that simulation estimation has a bias reduction effect like the bootstrap. Hence in the econometrics literature, SMD estimators are used even when the likelihood is tractable, as in gpy.
The steps for implementing the SMD are as follows:
The SMD is the $\theta$ that makes $J_S(\theta)$ smaller than the tolerance specified for the numerical optimizer. In the exactly identified case, the tolerance can be made as small as machine precision permits. When $\hat\psi$ is a vector of unconditional moments, the SMM estimator of duffie-singleton is obtained. When $\hat\psi$ are parameters of an auxiliary model, we have the `indirect inference' estimator of gmr. These are Wald-test based SMD estimators in the terminology of smith-palgrave. When $\hat\psi$ is the score function associated with the likelihood of the auxiliary model, we have the EMM estimator of gallant-tauchen-emm, which can also be thought of as an LM-test based SMD. If $\hat\psi$ is the likelihood of the auxiliary model, $J_S(\theta)$ can be interpreted as a likelihood ratio and we have a LR-test based SMD. g-monfort-simulation provide a framework that unifies these three approaches to SMD estimation. nickl-potscher show that an SMD based on non-parametrically estimated auxiliary statistics can have asymptotic variance equal to the Cramer-Rao bound if the tuning parameters are optimally chosen.
The Wald, LM, and LR based SMD estimators minimize a weighted $\mathsf L_2$ distance between the data and the model as summarized by auxiliary statistics. creel-kristensen-il consider a class of estimators that minimize the Kullback-Leibler distance between the model and the data.\footnote{ In the sequel, we take the more conventional $\mathsf L_2$ definition of SMD as given above.} Within this class, their MIL estimator maximizes an `indirect likelihood', defined as the likelihood of the auxiliary statistics. Their BIL estimator uses Bayesian computations to approximate the mode of the indirect likelihood. In practice, the indirect likelihood is unknown. Estimating it by kernel smoothing of the simulated statistics, the SBIL estimator combines Bayesian computations with non-parametric estimation. gao-hong show that using local linear regressions instead of kernel estimation can reduce the variance and the bias. Using non-parametric estimation in ABC has previously been considered in beaumont-zhang-balding. cghk:16 show that not only can such an ABC implementation bypass MCMC altogether, it can provide asymptotically valid frequentist inference. Bounds for the number of simulations that achieve the parametric rate of convergence and asymptotic normality are derived.
The ABC literature often credits Donald Rubin to be the first to consider the possibility of estimating the posterior distribution when the likelihood is intractable. diggle-gratton-84 propose to approximate the likelihood by simulating the model at each point on a parameter grid and appear to be the first implementation of simulation estimation for models with intractable likelihoods. Subsequent developments adapted the idea to conduct posterior inference, giving the prior an explicit role. The first ABC algorithm was implemented by tbfd and pspf:99 to study population genetics. Their Accept/Reject algorithm is as follows: (i) draw $\theta^b$ from the prior distribution $\pi(\theta)$, (ii) simulate data using the model under $\theta^b$ (iii) accept $\theta^b$ if the auxiliary statistics computed using the simulated data are close to $\hat\psi$. As in the SMD literature, the auxiliary statistics can be parameters of a regression or unconditional sample moments. heggland-frigessi, drovandi-pettitt-faddy,drovandi-15 use simulated auxiliary statistics.
Since simulating from a non-informative prior distribution is inefficient, subsequent work suggests to replace the rejection sampler by one that takes into account the features of the posterior distribution. The likelihood of the full dataset $L(y|\theta)$ is intractable, as is the likelihood of the finite dimensional statistic $L(\hat \psi|\theta)$. However, the latter can be consistently estimated using simulations. The general idea is to set as a target the intractable posterior density \[p^*_{ABC}(\theta|\hat \psi) \propto \pi(\theta)L (\hat \psi|\theta)\] and approximate it using Monte-Carlo methods. Some algorithms are motivated from the perspective of non-parametric density estimation, while others aim to improve properties of the Markov chain.\footnote { Recent surveys on ABC can be found in mprr-12, blum-nunes-prangle-sisson among others. See drovandi-15,drovandi-pettitt-faddy for differences amongst ABC estimators.} The main idea is, however, using data augmentation to consider the joint density $p_{ABC}(\theta,x|\hat\psi)\propto L(\hat\psi|x,\theta)L(x|\theta)\pi(\theta)$, putting more weight on the draws with $x$ close to $\hat\psi$. When $x=\hat\psi$, $L(\hat\psi|\hat\psi,\theta)$ is a constant, $p_{ABC} (\theta,\hat\psi|\hat\psi)\ \propto L(\hat\psi|\theta)\pi(\theta)$, and the target posterior is recovered. If $\hat\psi$ are sufficient statistics, one recovers the posterior distribution associated with the intractable likelihood $L(\theta|y)$, not just an approximation.
To better understand the ABC idea and its implementation, we will write $\mathbf y^{b}$ instead of $\mathbf y^{b} (\varepsilon^{b},\theta^{b})$ and $\hat \psi^{b}$ instead of $\hat\psi^{b}(\mathbf y^{b} (\varepsilon^{b},\theta^{b}))$ to simplify notation. Let $\mathbb K_\delta(\hat \psi^b,\hat \psi|\theta)\geq 0$ be a kernel function that weighs deviations between $\hat\psi$ and $\hat\psi^b$ over a window of width $\delta$. Suppose we keep only the draws that satisfy $\hat \psi^b=\hat \psi$ and hence $\delta=0$. Note that $\mathbb K_0(\hat\psi^b,\hat\psi|\theta)=1$ if $\hat\psi=\hat\psi^b$ for any choice of the kernel function. Once the likelihood of interest \[ L(\hat \psi|\theta) = \int L(x|\theta)\mathbb K_0(x,\hat \psi|\theta)dx \] is available, moments and quantiles can be computed. In particular, for any measurable function $\varphi$ whose expectation exists, we have:
Since $\hat \psi^b|\theta^b \sim L(\cdot|\theta^b)$, the expectation can be approximated by averaging over draws from $L(\cdot|\hat\theta^b)$. More generally, draws can be taken from an importance density $q(\cdot)$. In particular, \[ \hat{\mathbb E}\left[\varphi (\theta)|\hat\psi=\hat\psi^b\right]=\frac{\sum_{b=1}^B \varphi(\theta^b) \mathbb K_0(\hat\psi^b,\hat \psi|\theta^b)\frac{\pi(\theta^b)}{q (\theta^b)}}{\sum_{b=1}^B \mathbb K_0(\hat\psi^b,\hat \psi|\theta^b)\frac{\pi(\theta^b)}{q(\theta^b)}}.\] The importance weights are then \[w_0^b \propto \mathbb K_0(\hat\psi^b,\hat \psi|\theta^b)\frac{\pi(\theta^b)}{q (\theta^b)}.\] By a law of large numbers, $\hat{\mathbb{E}} \left[ \varphi(\theta)|\hat \psi\right]\rightarrow \mathbb {E} \left[ \varphi(\theta)|\hat \psi\right]$ as $B\rightarrow\infty$.
There is, however, a caveat. When $\hat \psi$ has continuous support, $\hat \psi^b = \hat \psi$ is an event of measure zero. Replacing $\mathbb K_0$ with $\mathbb K_\delta$ where $\delta$ is close to zero yields the approximation:
Since $\mathbb K_\delta(\cdot)$ is a kernel function, consistency of the non-parametric estimator for the conditional expectation of $\varphi(\theta)$ follows from, for example, pagan-ullah. This is the approach considered in beaumont-zhang-balding, creel-kristensen-il and gao-hong. The case of a rectangular kernel $\mathbb K_\delta (\hat \psi,\hat \psi^b) = \mathbbm I_{\|\hat \psi-\hat\psi^b\|\leq \delta}$ corresponds to the ABC algorithm proposed in mmpt-03. This is the first ABC algorithm that exploits MCMC sampling. Hence we refer to it as MCMC-ABC. Our analysis to follow is based on this algorithm. Accordingly, we now explore it in more detail.
\paragraph{Algorithm MCMC-ABC} Let $q (\cdot)$ be the proposal distribution. For $b=1,\ldots, B$ with $\theta^0$ given,
As with all ABC algorithms, the success of the MCMC-ABC lies in augmenting the posterior with simulated data $\hat \psi^b$, i.e. $p^*_{ABC}(\theta^b,\hat\psi^b|\hat\psi)\propto L (\hat\psi|\theta^b,\hat\psi^b)L (\hat\psi^b|\theta^b)\pi(\theta^b)$. The joint posterior distribution that the MCMC-ABC would like to target is \[ p^0_{\text{ABC}}\left( \theta^b, \hat \psi^b | \hat \psi \right) \propto \pi(\theta^b)L(\hat \psi^b|\theta^b)\mathbbm I_{\|\hat \psi^b-\hat \psi\| = 0} \] since integrating out $\varepsilon^b$ would yield $p^*_{ABC}(\theta|\hat\psi)$. But it would not be possible to generate draws such that $\|\hat\psi^b-\hat\psi\|$ equals zero exactly. Hence as a compromise, the MCMC-ABC algorithm allows $\delta>0$ and targets \[ p^\delta_{\text{ABC}}\left( \theta^b, \hat \psi^b | \hat \psi \right) \propto \pi(\theta^b)L(\hat \psi^b|\theta^b)\mathbbm I_{\|\hat \psi^b-\hat \psi\| \leq \delta}. \] The adequacy of $p_{ABC}^\delta$ as an approximation of $p^0_{ABC}$ is a function of the tuning parameter $\delta$.
To understand why this algorithm works, we follow the argument in sisson-fan. If the initial draw $\theta^1$ satisfies $\|\hat \psi - \hat \psi^1 \| \leq \delta$, then all subsequent $b>1$ draws are such that $ \mathbbm I_{\| \hat \psi^b - \hat \psi \| \leq \delta } =1$ by construction. Furthermore, since we draw $\theta^{b+1}$ and then independently simulate data $\hat \psi^{b+1}$, the proposal distribution becomes $ q(\theta^{b+1},\hat\psi^{b+1}|\theta^b) = q(\theta^{b+1}|\theta^b)L(\hat\psi^ {b+1}|\theta^{b+1}). $ The two observations together imply that
The last equality shows that the acceptance ratio is in fact the ratio of two ABC posteriors times the ratio of the proposal distribution. Hence the MCMC-ABC effectively targets the joint posterior distribution $p_{ABC}^\delta$.
Thus far, we have seen that the SMD estimator is the $\theta$ that makes $\|\hat\psi-\frac{1}{S}\sum_{s=1}^S\hat\psi^s(\theta)\|$ no larger than the tolerance of the numerical optimizer. We have also seen that the feasible MCMC-ABC accepts draws $\theta^b$ satisfying $\|\hat\psi-\hat\psi^b(\theta^b)\|\leq \delta $ with $\delta>0$. To view the MCMC-ABC from a different perspective, suppose that setting $\delta=0$ was possible. Then each accepted draw $\theta^b$ would satisfy: \[ \hat \psi^b(\theta^b)=\hat \psi. \] For fixed $\varepsilon^b$ and assuming that the mapping $\hat \psi^b : \theta \rightarrow \hat \psi^b(\theta)$ is continuously differentiable and one-to-one, the above statement is equivalent to: \[ \theta^b = \text{argmin}_\theta \left( \hat \psi^b(\theta)-\hat \psi\right)^\prime \left( \hat \psi^b(\theta)-\hat \psi\right). \] Hence each accepted $\theta^b$ is the solution to a SMD problem with $S=1$. Next, suppose that instead of drawing $\theta^b$ from a proposal distribution, we draw $\varepsilon^b$ and solve for $\theta^b$ as above. Since the mapping $\hat \psi^b$ is invertible by assumption, a change of variable yields the relation between the distribution of $\hat\psi^b$ and $\theta^b$. In particular, the joint density, say $h(\theta^b,\varepsilon^b)$, is related to the joint density $L(\hat\psi^b(\theta^b),\varepsilon^b)$ via the determinant of the Jacobian $|\hat \psi^b_\theta (\theta^b)|$ as follows: \[ h(\theta^b,\varepsilon^b|\hat \psi) = |\hat \psi^b_\theta (\theta^b)|L(\hat \psi^b(\theta^b), \varepsilon^b|\hat \psi). \] Multiplying the quantity on the right-hand-side by $w^b(\theta^b)=\pi(\theta^b)|\hat \psi^b_\theta (\theta^b)|^{-1}$ yields $\pi(\theta^b)L(\hat\psi,\varepsilon^b|\theta^b)$ since $\hat\psi^b(\theta^b)=\hat\psi$ and the mapping from $\theta^b$ to $\psi^b(\theta^b)$ is one-to-one. This suggests that if we solve the SMD problem $B$ times each with $S=1$, re-weighting each of the $B$ solutions by $w^b(\theta^b)$ would give the target the joint posterior $p_{ABC}^*(\theta|\hat\psi)$ after integrating out $\varepsilon^b$.
\paragraph{Algorithm RS }
The RS has the optimization aspect of SMD as well as the sampling aspect of the MCMC-ABC. We call the RS the reverse sampler for two reasons. First, typical Bayesian estimation starts with an evaluation of the prior probabilities. The RS terminates with the evaluation of the prior. Furthermore, we use the SMD estimates to reverse engineer the posterior distribution.
Consistency of each RS solution (i.e. $\theta^b$) is built on the fact that the SMD is consistent even with $S=1$. The RS estimate is thus an average of a sequence of SMD modes. In contrast, the SMD is the mode of an objective function defined from a weighted average of the simulated auxiliary statistics. Optimization effectively allows $\delta$ to be as close to zero as machine precision permits. This puts the joint posterior distribution as close to the infeasible target as possible, but has the consequence of shifting the distribution from $(\mathbf y^b,\hat\psi^b)$ to $(\mathbf y^b,\theta^b)$. Hence a change of variable is required. The importance weight depends on the Jacobian matrix, making the RS an optimization based importance sampler.
The proof is given in jjng-15. By convergence, we mean that for any measurable function $\varphi(\theta)$ such that the expectation exists, a law of large numbers implies that \newline $ \sum_{b=1}^B \bar w^b(\theta^b) \varphi(\theta^b) \asconv \mathbb E_{p^*(\theta|\hat\psi)}(\varphi(\theta))$. In general, $\bar w^b(\theta^{b})\ne \frac1B$. The RS draws and moments can be interpreted as if they were taken from $p^*_{\text {ABC}}$, the posterior distribution had the likelihood $p(\hat\psi|\theta)$ been available.
That the draws of the MCMC-ABC at $\delta=0$ can be seen from an optimization perspective allows us to subsequently use the RS as a conceptual framework to understand the differences between the ideal MCMC-ABC and SMD. It should be noted that the RS is not the same as the MCMC-ABC or any ABC estimator implemented with $\delta>0$ as they necessarily have an acceptance rate strictly less than one. Indeed, a challenge of many ABC implementations is the low acceptance rate. The RS draws are always accepted and can be useful in situations when numerical optimization of the auxiliary model is easy. Properties of the RS are further analyzed in jjng-15. meeds-welling independently propose an ABC sampling algorithm similar to the RS. Their focus is on ways to implement it efficiently using embarrassingly parallel methods.
The GMM objective function $J(\theta)$ defined in ((ref)) is not a proper density. Noting that $\exp(-J(\theta))$ is the kernel of the Gaussian density, jiang-turnbull:94 define an {\em indirect likelihood} (distinct from the one defined in creel-kristensen-il) as \[ L_{IND}(\theta|\hat\psi) \equiv \frac{1}{\sqrt{2\pi}} |\Sigma|^{-1} \exp( - J(\theta)).\] Associated with the indirect likelihood is the indirect score, indirect Hessian, and a generalized information matrix equality, just like a conventional likelihood. Though the indirect likelihood is not a proper density, its maximizer has properties analogous to the maximum likelihood estimator provided by $\mathbb E[g_t (\theta_0)]=0$.
In chernozhukov-hong, the authors observe that extremum estimators can be difficult to compute if the objective function is highly non-convex, especially when the dimension of the parameter space is large. These difficulties can be alleviated by using Bayesian computational tools, but this is not possible when the objective function is not a likelihood. chernozhukov-hong take an exponential of $-J(\theta)$, as in jiang-turnbull:94, but then combine $\exp(-J(\theta))$ with a prior density $\pi(\theta)$ to produce a quasi-posterior density. Chernozhukov and Hong initially termed their estimator `Quasi-Bayes' because $\exp(-J(\theta))$ is not a standard likelihood. They settled on the term `Laplace-type estimator' (LT), so-called because Laplace suggested to approximate a smooth pdf with a well defined peak by a normal density, see tierney-kadane:86. If $\pi(\theta)$ is strictly positive and continuous over a compact parameter space $\Theta$, the `quasi-posterior' LT distribution
is proper. The LT posterior mean is thus well-defined even when the prior may not be proper. As discussed in chernozhukov-hong, one can think of the LT under a flat prior as using simulated annealing to maximize $\exp(-J(\theta))$ and setting the cooling parameter $\tau$ to 1. Frequentist inference is asymptotically valid because as the sample size increases, the prior is dominated by the pseudo likelihood which, by the Laplace approximation, is asymptotically normal.\footnote{ For loss function $d(\cdot)$, the LT estimator is $ \hat\theta_{LT}(\vartheta)=\argmin_\theta \int_\Theta d(\theta-\vartheta) p_{LT}(\theta|\mathbf y)d\theta. $ If $d(\cdot)$ is quadratic, the posterior mean minimizes quasi-posterior risk.}
In practice, the LT posterior distribution is targeted using MCMC methods. Upon replacing the likelihood $L(\theta)$ by $\exp(-J(\theta))$, the MH acceptance probability is \[\rho_{LT}(\theta^b,\vartheta) = \min \Big( \frac{\exp(-J(\vartheta))\pi(\vartheta)q(\theta^b|\vartheta)} {\exp( -J(\theta^ b))\pi(\theta^b)q(\vartheta|\theta^b)},1 \Big).\] The quasi-posterior mean is $ \bar \theta_{LT}= \frac{1}{B} \sum_{b=1}^B \theta^b$ where each $\theta^b$ is a draw from $p_{LT}(\theta|\mathbf y)$. Chernozhukov and Hong suggest to exploit the fact that the quasi-posterior mean is much easier to compute than the mode and that, under regularity conditions, the two are first order equivalent. In practice, the weighting matrix can be based on some preliminary estimate of $\theta$, or estimated simultaneously with $\theta$. In exactly identified models, it is well known that the MD estimates do not depend on the choice of $W$. This continues to be the case for the LT posterior mode $\hat\theta_{LT}$. However, the posterior mean is affected by the choice of the weighting matrix even in the just-identified case.\footnote{kormiltsina-nekipelov:14 suggests to scale the objective function to improve coverage of the confidence intervals.}
The LT estimator is built on the validity of the asymptotic normal approximation in the second order expansion of the objective function. nekipolov-kormilitsina:15 show that in small samples, this approximation can be poor so that the LT posterior mean may differ significantly from the extremum estimate that it is meant to approximate. To see the problem in a different light, we again take an optimization view. Specifically, the asymptotic distribution $\sqrt{T}(\hat\psi(\theta_0)-\psi(\theta_0))\dconv \mathcal N (0,\Sigma (\theta_0))\equiv \mathbb A_\infty(\theta_0)$ suggests to use \[\hat\psi^b(\theta) \approx \psi(\theta) +\frac{\mathbb A^b_\infty(\theta_0)}{\sqrt{T}}\] where $\mathbb A^b_\infty(\theta_0) \sim \mathcal N (0,\hat\Sigma(\theta))$. Given a draw of $\mathbb A^b_\infty$, there will exist a $ \theta^b$ such that $ (\hat\psi^b(\theta)-\hat\psi)^\prime W (\hat\psi^b(\theta)-\hat\psi )$ is minimized. In the exactly identified case, this discrepancy can be driven to zero up to machine precision. Hence we can define \[ \theta^b= \argmin_\theta \|\hat\psi^b(\theta)-\hat\psi\|.\] Arguments analogous to the RS suggest the following will produce draws of $\theta$ from $p_{LT}(\theta|\mathbf y)$.
Seen from an optimization perspective, the LT is a weighted average of MD modes with the determinant of the Jacobian matrix as importance weight, similar to the RS. It differs from the RS in that the Jacobian here is computed from the asymptotic binding function $\psi(\theta)$, and the draws are based on the asymptotic normality of $\hat\psi$. As such, simulation of the structural model is not required.
When $\psi(\theta )$ is not analytically tractable, a natural modification is to approximate it by simulations as in the SMD. This is the approach taken in lise-meghir-robin. We refer to this estimator as the Simulated Laplace-type estimator, or SLT. The steps are as follows:
The SLT algorithm has two loops, one using $S$ simulations for each $b$ to approximate the asymptotic binding function, and one using $B$ draws to approximate the `quasi-posterior' SLT distribution
The above SLT algorithm has features of SMD, ABC, and LT, it also requires simulations of the full model. As a referee pointed out, though the SLT resembles the ABC algorithm when used with a Gaussian kernel, $\exp (-J_S(\theta))$ is not a proper density, and $p_{SLT}(\theta|\mathbf y,\varepsilon^1,\ldots,\varepsilon^S)$ is not a conventional likelihood-based posterior distribution. While the SLT targets the pseudo likelihood, ABC algorithms target the proper but intractable likelihood. Furthermore, the asymptotic distribution of $\hat\psi$ is known from a frequentist perspective. In ABC estimation, lack of knowledge of the likelihood of $\hat\psi$ motivates the Bayesian computation.
The optimization implementation of SLT presents a clear contrast with the ABC.
While the SLT is a weighted average of SMD modes, the draws of $\hat\psi^b(\theta)$ are taken from the (frequentist) asymptotic distribution of $\hat\psi$ instead of solving the model at each $b$. gao-hong use a similar idea to make draws of what we refer to as $\bar g(\theta)$ in their extension of the BIL estimator of creel-kristensen-il to non-separable models.
The SMD, RS, ABC, and SLT all require specification and simulation of the full model. At a practical level, the innovations $\varepsilon^1,\dots, \varepsilon^s$ used in SMD and SLT are only drawn from $F_\varepsilon$ once and held fixed across iterations. Equivalently, the seed of the random number generator is fixed so that the only difference in successive iterations is due to change in the parameters to be estimated. In contrast, ABC draws new innovations from $ F_\varepsilon$ each time a $\theta^{b+1}$ is proposed. We need to simulate $B$ sets of innovations of length $T$, not counting those used in draws that are rejected, and $B$ is generally much bigger than $S$. The SLT takes $B$ draws from an asymptotic distribution of $\hat\psi$. Hence even though some aspects of the algorithms considered seem similar, there are subtle differences.
This section studies the finite sample properties of the various estimators. Our goal is to compare the SMD with the RS, and by implication, the infeasible MCMC-ABC. Note that our RS is different from the original kernel based ABC methods. To do so in a tractable way, we only consider the expansion up to order $\frac{1}{T}$. As a point of reference, we first note that under assumptions in rsu-96,bao-ullah:07, $\hat\theta_{ML}$ admits a second order expansion \[ \hat\theta_{ML}=\theta_0+\frac{A_{ML}(\theta_0)}{\sqrt{T}}+\frac{C_{ML}(\theta_0)}{T}+o_p(\frac1T). \] where $A_{ML}(\theta_0)$ is a mean-zero asymptotically normal random vector and $C_{ML}(\theta_0)$ depends on the curvature of the likelihood. These terms are defined as
where the normalized score $\frac{1}{\sqrt{T}} \ell_\theta(\theta_0)$ and centered Hessian $\frac{1}{\sqrt{T}} ( \ell_{\theta\theta}(\theta_0)-\mathbb E[\ell_{\theta\theta}(\theta_0)])$ converge in distribution to the normal vectors $Z_{S}$ and $Z_H$ respectively. The order $\frac1T$ bias is large when Fisher information is low.
Classical Bayesian estimators are likelihood based. Hence the posterior mode $\hat\theta_{ BC}$ exhibits a bias similar to that of $\hat\theta_{ML}$. However, the prior $\pi(\theta)$ can be thought of as a constraint, or penalty since the posterior mode maximizes $\log p(\theta|\mathbf y)= \log L(\theta|\mathbf y)+\log \pi(\theta)$. Furthermore, kass-tierney-kadane show that the posterior mean deviates from the posterior mode by a term that depends on the second derivatives of the log-likelihood. Accordingly, there are three sources of bias in the posterior mean $\bar\theta_{BC}$: a likelihood component, a prior component, and a component from approximating the mode by the mean. Hence
Note that the prior component is under the control of the researcher.
In what follows, we will show that posterior means based on auxiliary statistics $\hat\psi$ generically have the above representation, but the composition of the terms differ.
Minimum distance estimators depend on auxiliary statistics $\hat\psi$. Its properties have been analyzed in newey-smith-04 within an empirical-likelihood framework. To facilitate subsequent analysis, we follow g-monfort-simulation and directly expand $\hat \psi$ around $\psi(\theta_0)$, under the assumption that it admits a second-order expansion. In particular, since $\hat\psi$ is $\sqrt{T}$ consistent for $\psi(\theta_0)$, $\hat\psi$ has expansion
It is then straightforward to show that the minimum distance estimator $\hat\theta_{MD}$ has expansion
The bias in $\hat\theta_{MD}$ depends on the curvature of the binding function and the bias in the auxiliary statistic $\hat\psi$, $\mathbbm C(\theta_0)$. Then following grt-99, we can analyze the SMD as follows. In view of ((ref)), we have, for each $s$:
The estimator $\hat\theta_{SMD}$ satisfies $\hat\psi= \frac{1}{S}\sum_{s=1}^S \hat \psi^s(\hat\theta_ {SMD})$ and has expansion $\hat\theta_{SMD}= \theta_0+\frac{A_{SMD}(\theta_0)}{\sqrt {T}}+\frac{C_{SMD}(\theta_0)}{T}+o_p(\frac1T)$. Plugging it in the Edgeworth expansions gives:
Expanding $\psi(\hat\theta_{SMD})$ and $\mathbbm A^s(\hat\theta_{SMD})$ around $\theta_0 $ and equating terms in the expansion of $\hat\theta_{SMD}$,
The first order term can be written as $A_{SMD}=A_{MD}+\frac{1}{B}[\psi_\theta(\theta_0)]^{-1}\sum_{b=1}^B \mathbb A^b(\theta_0)$, the last term has variance of order $1/B$ which accounts for simulation noise. Note also that $\mathbb{E}\left( \frac{1}{S}\sum_{s=1}^S \mathbbm C^s(\theta_0) \right) = \mathbb E[\mathbbm C (\theta_0)]$. Hence, unlike the MD, $\mathbb E[C_{SMD}(\theta_0)]$ does not depend on the bias $\mathbb{C}(\theta_0)$ in the auxiliary statistic. In the special case when $\hat\psi$ is a consistent estimator of $\theta_0$, $\psi_ {\theta}(\theta_0)$ is the identity map and the term involving $\psi_{\theta\theta_j} (\theta_0)$ drops out. Consequently, the SMD has no bias of order $\frac{1}{T}$ when $S\rightarrow\infty$ and $\psi(\theta)=\theta$. In general, the bias of $\hat\theta_ {SMD}$ depends on the curvature of the binding function as
This is an improvement over $\hat\theta_{MD}$ because as seen from ((ref)),
The bias in $\hat\theta_{MD}$ has an additional term in $\mathbb C(\theta_0)$.
The convergence properties of the ABC algorithms have been well analyzed but the theoretical properties of the estimates are less understood. dsjp establish consistency of the ABC in the case of hidden Markov models. The analysis considers a scheme so that maximum likelihood estimation based on the ABC algorithm is equivalent to exact inference under the perturbed hidden Markov scheme. The authors find that the asymptotic bias depends on the ABC tolerance $\delta$. calvet-czellar:14 provide an upper bound for the mean-squared error of their ABC filter and study how the choice of the bandwidth affects properties of the filter. Under high level conditions and adopting the empirical likelihood framework of newey-smith-04, creel-kristensen-il show that the infeasible BIL is second order equivalent to the MIL after bias adjustments, while MIL is in turn first order equivalent to the continuously updated GMM. The feasible SBIL (which is also an ABC estimator) has additional errors compared to the BIL due to simulation noise and kernel smoothing, but these errors vanish as $S\rightarrow\infty$ for an appropriately chosen bandwidth. gao-hong show that local-regressions have better variance properties compared to kernel estimations of the indirect likelihood. cghk:16 show that the number of simulations can affect the parametric convergence rate and asymptotic normality of the estimator, which is important for frequentist inference.
ABC algorithms are traditionally implemented using kernel smoothing, the first implementation being beaumont-zhang-balding. The bias due to kernel smoothing is rigorously studied in cghk:16 under the assumption that the draws are taken directly from the prior. Our RS is an importance sampler that does not use kernel smoothing. Instead it uses optimization to set $\delta$ equal to zero. This offers different insight as we look at the bias in the ideal case where $\delta$ is exactly zero.
As shown above, $\bar\theta_ {RS }$ is the weighted average of a sequence of SMD modes. Analysis of the weights $w^b(\theta^b)$ requires an expansion of $\hat\psi^b_\theta (\theta^b)$ around $\psi_\theta(\theta_0)$. From such an analysis, shown in the Appendix, we find that
where
The SMD and RS are first order equivalent, but $\bar\theta_{RS }$ has an order $\frac {1}{T}$ bias. The bias, given by $C_{RS}(\theta_0)$, has three components. The $C^M_ {RS}(\theta_0)$ term (defined in Appendix A) can be traced directly to the weights, or to the interaction of the weights with the prior, and is a function of $A_{RS}(\theta_0)$. Some but not all the terms vanish as $B\rightarrow \infty$. The second term will be zero if a uniform prior is chosen since $\pi_\theta=0$. A similar result is obtained in creel-kristensen-il. The first term is
{
} The term $ \mathbbm C(\theta_0)-\frac{1}{B}\sum_{b=1}^B\mathbbm C^b(\theta_0)$ is exactly the same as in $C_{SMD}(\theta_0)$. The middle term involves $\psi_{\theta\theta_j}(\theta_0)$ and is zero if $\psi(\theta)=\theta$. But because the summation is over $\theta^b$ instead of $\hat \psi^s$,
As a consequence $ \mathbb E[C_{RS}(\theta_0)]\ne 0$ even when $\psi(\theta)=\theta$. In contrast, $\mathbb E[ C_{SMD}(\theta_0)]=0$ when $\psi(\theta)=\theta$ as seen from ((ref)). The reason is that the comparable term in $C_{SMD}(\theta_0)$ is
The difference boils down to the fact that the SMD is the mode of the average over simulated auxiliary statistics, while the RS is a weighted average over the modes. As will be seen below, this difference is also present in the LT and SLT and comes from averaging over $\theta^b$. The result is based on fixing $\delta$ at zero and holds for any $B$. Proposition (ref) implies that the ideal MCMC-ABC with $\delta=0$ also has a non-negligible second-order bias. Note that Proposition (ref) is stated for the exactly identified case. When $dim(\hat \psi)>dim(\theta)$, the analysis is more complicated. Essentially, when the model is overidentified, weighting is needed since all moments cannot be made equal to zero simultaneously in general. This introduces additional biases. A result analogous to Proposition (ref) is given in jjng-15 for the overidentified case.
In theory, the order $\frac{1}{T}$ bias can be removed if $\pi(\theta)$ can be found to put the right hand side of $C^{RS}(\theta_0)$ defined in ((ref)) to zero. Then $\bar\theta_{RS }$ will be second order equivalent to SMD when $\psi (\theta)=\theta$ and may have a smaller bias than SMD when $\psi(\theta)\ne \theta$ since SMD has a non-removable second order bias in that case. That the choice of prior will have bias implications for likelihood-free estimation echoes the findings in the parametric likelihood setting. ArellanoBonhomme show in the context of non-linear panel data models that the first-order bias in Bayesian estimators can be eliminated with a particular prior on the individual effects. bester-hansen-06 also show that in the estimation of parametric likelihood models, the order $\frac{1}{T}$ bias in the posterior mode and mean can be removed using objective Bayesian priors. They suggest to replace the population quantities in a differential equation with sample estimates. Finding the bias-reducing prior for the RS involves solving the differential equation: \[0= \mathbb{E}[C_{RS}^b(\theta_0)] + \frac {\pi_\theta(\theta_0)} {\pi(\theta_0)}\mathbb{E}[(A_{RS}^b(\theta_0)-\bar A_{RS}(\theta_0))A_{RS}^b(\theta_0)] +\mathbb{E}[C^M_{RS}(\theta_0),\pi(\theta_0)] \] which has the additional dependence on $\pi$ in $C^M_{RS}(\theta_0,\pi (\theta_0))$ that is not present in bester-hansen-06. A closed-form solution is available only for simple examples as we will see Section 6.1 below. For realistic problems, how to find and implement the bias-reducing prior is not a trivial problem. A natural starting point is the plug-in procedure of bester-hansen-06 but little is known about its finite sample properties even in the likelihood setting for which it was developed.
This section has studied the RS, which is the best that the MCMC-ABC can achieve in terms of $\delta$. This enables us to make a comparison with the SMD holding the same $\mathsf L_2$ distance between $\hat\psi$ and $\psi(\theta)$ at zero by machine precision. However, the MCMC-ABC algorithm with $\delta>0$ will not produce draws with the same distribution as the RS. To see the problem, suppose that the RS draws are obtained by stopping the optimizer before $\|\hat\psi-\psi(\theta^b)\|$ reaches the tolerance guided by machine precision. This is analogous to equating $\psi(\theta^b)$ to the pseudo estimate $\hat\psi+\delta$. Inverting the binding function will yield an estimate of $\theta$ that depends on the random $\delta$ in an intractable way. The RS estimate will thus have an additional bias from $\delta\ne 0$. By implication, the MCMC-ABC with $\delta>0$ will be second order equivalent to the SMD only after a bias adjustment even when $\psi(\theta)=\theta$.
The mode of $\exp(-J(\theta))\pi(\theta)$ will inherit the properties of a MD estimator. However, the quasi-posterior mean has two additional sources of bias, one arising from the prior, and another one from approximating the mode by the mean. The optimization view of $\bar\theta_{LT}$ facilitates an understanding of these effects. As shown in Appendix B, each draw $\theta^b_{LT}$ has expansion terms
Even though the LT has the same objective function as MD, simulation noise enters both $A^b_{LT}(\theta_0)$ and $C^b_{LT}(\theta_0)$. Compared to the extremum estimate $\hat\theta_{MD}$, we see that $ A_{LT}= \frac{1}{B}\sum_{b=1}^B A_{LT}^b(\theta_0)\ne A_{MD}(\theta_0)$ and $ C_{LT}(\theta_0)\ne C_{MD}(\theta_0)$. Although $C_{LT}(\theta_0)$ has the same terms as $C_{RS}(\theta_0)$, they are different because the LT uses the asymptotic binding function, and hence $A^b_ {LT}(\theta_0)\ne A^b_{RS}(\theta_0)$.
A similar stochastic expansion of each $\theta^b_{SLT}$ gives:
Following the same argument as in the RS, an optimally chosen prior can reduce bias, at least in theory, but finding this prior will not be a trivial task. Overall, the SLT has features of the RS (bias does not depend on $\mathbbm {C}(\theta_0)$) and the LT (dependence on $\mathbbm{A}^b_\infty$) but is different from both. Because the SLT uses simulations to approximate the binding function $\psi(\theta)$, $\mathbb E[\mathbb C(\theta_0)-\frac{1}{S}\sum_{s=1}^S \mathbb C^s(\theta_0)]=0$. The improvement over the LT is analogous to the improvement of SMD over MD. However, the $A^b_{SLT}(\theta_0)$ is affected by estimation of the binding function (the term with superscript $s$) and of the quasi-posterior density (the terms with superscript $b$). This results in simulation noise with variance of order $1/S$ plus another of order $1/B$. Note also that the SLT bias has an additional term \[ \frac{1}{B}\sum_{b=1}^B \left(\frac{1}{S}\sum_{s=1}^S (\mathbbm{A}_{\theta}^s(\theta_0)+\mathbbm {A}_{\infty,\theta}^b (\theta_0)) A^b_{SLT}(\theta_0)\right) \overset{S \to \infty}{\to} \frac{1}{B}\sum_{b=1}^B \mathbbm {A}_{\infty,\theta}^b (\theta_0) A^b_{LT}(\theta_0).\] The main difference with the RS is that $\mathbbm{A}^b$ is replaced with $\mathbbm{A}^b_\infty$. For $S=\infty$ this term matches that of the LT.
We started this section by noting that the Bayesian posterior mean has two components in its bias, one arising from the prior which acts like a penalty on the objective function, and another due to approximating the mean with the mode. We are now in a position to use the results in the foregoing subsections to show that for $d$=(MD, SMD, RS, LT) and SLT and $D=$ (RS,LT,SLT) these estimators can be represented as
where with $ A_d^b(\theta_0) = [\psi_\theta(\theta_0)]^{-1} \Big( \mathbbm A (\theta_0) - \mathbbm A_d^b(\theta_0) \Big)$,
The term $C^P_d(\theta_0)$ is a bias directly due to the prior. The term $C^M_d(\theta_0)$, defined in the Appendix, depends on $A_d(\theta_0)$, the curvature of the binding function, and their interaction with the prior. Hence at a general level, the estimators can be distinguished by whether or not Bayesian computation tools are used, as the indicator function is null only for the two frequentist estimators (MD and SMD). More fundamentally, the estimators differ because of $A_d (\theta_0)$ and $C_d(\theta_0)$, which in turn depend on $\mathbb A^b_d(\theta_0)$ and $\mathbb C_d(\theta_0)$. We compactly summarize the differences as follows:
The MD is the only estimator that is optimization based and does not involve simulations. Hence it does not depend on $b$ or $s$ and has no simulation noise. The SMD does not depend on $b$ because the optimization problem is solved only once. The LT simulates from the asymptotic binding function. Hence its errors are associated with parameters of the asymptotic distribution.
The MD and LT have a bias due to asymptotic approximation of the binding function. In such cases, cabrera-fernholz suggest to adjust an initial estimate $\tilde\theta$ such that if the new estimate $\hat\theta$ were the true value of $\theta$, the mean of the original estimator equals the observed value $\tilde\theta$. Their {\em target estimator} is the $\theta$ such that $\mathbb E_{\mathcal P_{\theta}}[\hat \theta]=\tilde \theta$. While the bootstrap directly estimates the bias, a target estimator corrects for the bias implicitly. cabrera-hu show that the bootstrap estimator corresponds to the first step of a target estimator. The latter improves upon the bootstrap estimator by providing more iterations.
An auxiliary statistic based target estimator is the $\theta$ that solves $\mathbb E_{\mathcal P _\theta}[ \hat\psi (\mathbf y( \theta))] =\hat\psi(\mathbf y(\theta_0))$. It replaces the asymptotic binding function $\lim_{T\rightarrow\infty} \mathbb E[\hat\psi(\mathbf y(\theta_0))]$ by $\mathbb E_{\mathcal P _\theta}[ \hat\psi (\mathbf y( \theta))]$ and approximates the expectation under $\mathcal P_\theta$ by stochastic expansions. The SMD and SLT can be seen as target estimators that approximate the expectation by simulations. Thus, they improve upon the MD estimator even when the binding function is tractable and is especially appealing when it is not. However, the improvement in the SLT is partially offset by having to approximate the mode by the mean.
The preceding section can be summarized as follows. A posterior mean computed through auxiliary statistics generically has a component due to the prior, and a component due to the approximation of the mode by the mean. The binding function is better approximated by simulations than asymptotic analysis. It is possible for simulation estimation to perform better than $\hat\psi_{MD}$ even if $\psi(\theta)$ were analytically and computationally tractable.
In this section, we first illustrate the above findings using a simple analytical example. We then evaluate the properties of the estimators using the dynamic panel model with fixed effects.
We consider the simple DGP $y_i\sim N( m,\sigma^2)$. The parameters of the model are $\theta=(m,\sigma^2)^\prime$. We focus on $\sigma^2$ since the estimators have more interesting properties.
The MLE of $\theta$ is \[\hat m=\frac{1}{T}\sum_{t=1}^T y_t, \quad\quad \hat\sigma^2= \frac{1}{T}\sum_{t=1}^T (y_t-\bar y)^2.\]
While the posterior distribution is dominated by the likelihood in large samples, the effect of the prior is not negligible in small samples. We therefore begin with a analysis of the effect of the prior on the posterior mean and mode in Bayesian analysis. Details of the calculations are provided in Appendix D.1.
We consider the prior $\pi(m,\sigma^2)= (\sigma^2)^{-\alpha} \mathbbm I_{\sigma^2>0} $, $\alpha>0$ so that the log posterior distribution is \[ \log p(\theta|y)=\log p(\theta|\hat m,\hat\sigma^2 )\propto \frac{-T}{2}\bigg [\log (2\pi \sigma^2) -\alpha \log \sigma^2 - \frac{1}{2\sigma^2}\sum_{t=1}^T (y_t-m)^2\bigg]\mathbbm I_{\sigma^2>0}.\] The posterior mode and mean of $\sigma^2$ are $ \sigma^2_{mode}= \frac{T \hat\sigma^2}{T+2\alpha} $ and $ \sigma^2_{mean} = \frac{T\hat\sigma^2}{T+2\alpha-5}.$ respectively. Using the fact that $E[\hat\sigma^2]= \frac{(T-1)}{T}\sigma^2$, we can evaluate $\sigma^2_{mode}$, $\sigma^2_{mean}$ and their expected values for different $\alpha$.
Two features are of note. For a given prior (here indexed by $\alpha$), the mean does not coincide with the mode. Second, the statistic (be it mean or mode) varies with $\alpha$. The Jeffrey's prior corresponds to $\alpha=1$, but the bias-reducing prior is $\alpha=2$. In the Appendix, we show that the bias reducing prior for this model is $\pi^R(\theta)\propto \frac{1}{\sigma^4}$.
Next, we consider estimators based on auxiliary statistics: \[ \hat\psi(\mathbf y)^\prime =
. \] As these are sufficient statistics, we can also consider (exact) likelihood-based Bayesian inference. For SMD estimation, we let $(\hat m_S, \hat \sigma^2_S)=(\frac{1}{S}\sum_{s=1}^S \hat m^ s, \frac{1}{S}\sum_ {s=1}^S \hat \sigma^{2,s})$. The LT quasi-likelihood using the variance of preliminary estimates of $m$ and $\sigma ^2$ as weights is: \[\exp(- J(m,\sigma^2)) = \exp\bigg(-\frac{T}{2} \bigg[ \frac{(\hat m-m)^2}{\hat\sigma^2} + \frac{ (\hat\sigma^2-\sigma^2)^2}{2\hat\sigma^4}\bigg]\bigg). \] The LT posterior distribution is $p(m,\sigma^2|\hat m,\hat \sigma^2)\propto \pi(m,\sigma^2)\exp(-J(m,\sigma^2))$. Integrating out $m$ gives $p(\sigma^2|\hat m,\hat\sigma ^2 )$. We consider a flat prior $\pi^U(\theta) \propto \mathbbm I_{\sigma^2 \geq 0}$ and the bias-reducing prior $\pi^R(\theta) \propto 1/\sigma^4\mathbbm I_{\sigma^2 \geq 0}$. The RS is the same as the SMD under a bias-reducing prior. Thus,
For completeness, the parametric Bootstrap bias corrected estimator $\hat\sigma^2_{\text{Bootstrap}}=2\hat\sigma^2 - \mathbb{E}_{\text{Bootstrap}} (\hat\sigma^2)$ is also considered:
$\mathbb{E}_{\text{Bootstrap}}(\hat\sigma^2)$ computes the expected value of the estimator replacing the true value $\sigma^2$ with $\hat \sigma^2$, the plug-in estimate. In this example the bias can be computed analytically since $\mathbb{E}(\hat \sigma^2(1+\frac{1}{T}))=\sigma^2(1-\frac{1}{T})(1+\frac{1}{T})=\sigma^2(1-\frac{1}{T^2}).$ While the bootstrap does not involve inverting the binding function, this computational simplicity comes at the cost of adding a higher order bias term (in $1/T^2$).
A main finding of this paper is that the reverse sampler can replicate draws from $p^*_{ABC}(\theta_0)$, which in turn equals the Bayesian posterior distribution if $\hat\psi$ are sufficient statistics. The weight for each SMD estimate is the prior times the Jacobian. To illustrate the importance of the Jacobian transformation, the top panel of Figure (ref) plots the Bayesian/ABC posterior distribution and the one obtained from the reverse sampler. They are indistinguishable. The bottom panel shows an incorrectly constructed reverse sampler that does not apply the Jacobian transformation. Notably, the two distributions are not the same.
The properties of the estimators are summarized in Table (ref). It should be reminded that increasing $S$ improves the approximation of the binding function in SMD estimation while increasing $B$ improves the approximation to the target distribution in Bayesian type estimation. For fixed $T$, only the Bayesian estimator with the bias reducing prior is unbiased. The SMD and RS (with bias reducing prior) have the same bias and mean-squared error in agreement with the analysis in the previous section. These two estimators have smaller errors than the RS estimator with a uniform prior. The SLT posterior mean differs from that of the SMD by $\kappa_ {SLT}$ that is not mean-zero. This term, which is a function of the Mills-ratio, arises as a consequence of the fact that the $\sigma^2$ in SLT are drawn from the normal distribution and then truncated to ensure positivity.
The dynamic panel model $y_{it}= \alpha_i + \rho y_{it-1} + \sigma e_{it}$ is known to be severely biased when $T$ is small because the unobserved heterogeneity $\alpha_i$ is imprecisely estimated. Various approaches have been suggested to improve the precision of the least squares dummy variable (LSDV) estimator $\hat\beta$.\footnote{See hsiao-book for a detailed account of this incidental parameter problem.} An interesting approach, due to gpy, is to exploit the bias reduction properties of the indirect inference estimator. Using the dynamic panel model as auxiliary equation, i.e. $\psi (\theta)=\theta$, the authors reported estimates of $\beta$ that are sharply more accurate than the LSDV, even when an exogenous regressor and a linear trend is added to the model. Their simulation experiments hold $\sigma^2$ fixed. We reconsider their exercise but also estimate $\sigma^2$.
With $\theta=(\rho,\beta,\sigma^2)^\prime$, we simulate data from the model:
Let $A=I_T-1_T1_T^{\prime}/T$ $\underline{A}=A \otimes I_T$, $\underline{y}=\underline A \; vec (y), \underline{y}_{-1}=\underline A\; vec(y_{-1}), \underline{x}=\underline A \;vec (x)$, where $y_{-1}$ are the lagged $y$. For this model, Bayesian inference is possible since the likelihood in de-meaned data is \[ L( \underline {\mathbf y},\underline {\mathbf x}|\theta)= \frac{1}{\sqrt{2\pi|\sigma^2\Omega|}^N}\exp \left( -\frac{1}{2\sigma^2} \sum_{i=2}^N (\underline y_i - \rho \underline y_{i,-1} - \beta \underline x_i)^\prime \Omega^{-1} (\underline y_i - \rho \underline y_{i,-1} - \beta \underline x_i)\right)\] where $\Omega = I_{T-1}-1_{T-1}1_{T-1}^{\prime}/T $. We use the following moment conditions for MD estimation:
with $\bar g(\hat\rho,\hat\beta,\hat\sigma^2)=0$. The simulated quantity $\bar g_S(\theta)$ for SMD and $\bar g^b(\theta)$ for ABC are defined analogously. The MD estimator in this case is also the LSDV. The auxiliary estimates for the ABC, RS, SLT and SMD are the LSDV estimates. Recall that while the weighting matrix $W$ is irrelevant to finding the mode in exactly identified models, $W$ affects computation of the posterior mean. We use $W = (\frac{1}{NT} \sum_ {i,t} g_{it}^\prime g_{it} - \bar{g}^\prime \bar{g})^{-1}$ for LT, MCMC-ABC, and SMD. The prior is $\pi(\theta)= \mathbbm I_{\sigma^2 \geq 0, \rho \in [-1,1], \beta \in \mathbb{R}}$. Since the demeaned data are used in LSDV estimation, the estimates are invariant to the specification of the fixed effects. Accordingly, we set them to zero both in the assumed DGP and the auxiliary model. The innovations $\varepsilon^s$ used to simulate the auxiliary model and to construct $\hat\psi^s$ are drawn from the standard normal distribution once and held fixed.
Table (ref) reports results from 5000 replications for $T=6$ time periods and $N=100$ cross-section units, as in gpy. Both $\hat\rho$ and $\hat\sigma^2$ are significantly biased. The LT is the same as the MD except that it is computed using Bayesian tools. Hence its properties are similar to the MD. The simulation estimators have much improved properties. The properties of $\bar\theta_{RS }$ are similar to those of the SMD. Figure (ref) illustrates for one simulated dataset how the posteriors for RS /SLT are shifted towards the true value compared to the one based on the direct likelihood.
The MCMC-ABC results in Table (ref) are for $\delta=0.10$ which has an acceptance rate of 0.58. These estimates are clearly more precise than MLE but more biased than SMD or RS. The dependence of MCMC-ABC on $\delta$ is investigated in further detail in jjng-15. In brief, when we set $\delta=0.25$, we achieve an acceptance ratio of 0.72 but the estimates are severely biased, as shown in Figure (ref). Bias similar to SMD and RS can be obtained if we set $\delta$ to 0.025. But the corresponding acceptance rate is 0.28, meaning that the MCMC-ABC needs at least three times more draws than the RS for a comparable level of bias. The choice of $\delta$ is more important for the properties of MCMC-ABC than the RS which associates $\delta$ with the tolerance of optimization.
Different disciplines have developed different estimators to overcome the limitations posed by an intractable likelihood. These estimators share many similarities: they rely on auxiliary statistics and use simulations to approximate quantities that have no closed form expression. We suggest an optimization framework that helps understand the estimators from the perspective of classical minimum distance estimation. All estimators are first-order equivalent as $S\rightarrow\infty$ and $T\rightarrow\infty$ for any choice of $\pi(\theta)$. Nonetheless, up to order $1/T$, the estimators are distinguished by biases due to the prior and approximation of the mode by the mean, the very two features that distinguish Bayesian and frequentist estimation.
We have only considered regular problems when $\theta_0$ is in the interior of $\Theta$ and the objective function is differentiable. When these conditions fail, the posterior is no longer asymptotically normal around the MLE with variance equal to the inverse of the Fisher Information Matrix. Understanding the properties of these estimators under non-standard conditions is the subject for future research.