Indirect likelihood inference
Abstract
Given a sample from a fully specified parametric model, let Zn be a given finite-dimensional statistic - for example, an initial estimator or a set of sample moments. We propose to (re-)estimate the parameters of the model by maximizing the likelihood of Zn. We call this the maximum indirect likelihood (MIL) estimator. We also propose a computationally tractable Bayesian version of the estimator which we refer to as a Bayesian Indirect Likelihood (BIL) estimator. In most cases, the density of the statistic will be of unknown form, and we develop simulated versions of the MIL and BIL estimators. We show that the indirect likelihood estimators are consistent and asymptotically normally distributed, with the same asymptotic variance as that of the corresponding efficient two-step GMM estimator based on the same statistic. However, our likelihood-based estimators, by taking into account the full finite-sample distribution of the statistic, are higher order efficient relative to GMM-type estimators. Furthermore, in many cases they enjoy a bias reduction property similar to that of the indirect inference estimator. Monte Carlo results for a number of applications including dynamic and nonlinear panel data models, a structural auction model and two DSGE models show that the proposed estimators indeed have attractive finite sample properties.
Full text
INDIRECT LIKELIHOOD INFERENCE MICHAEL CREEL AND DENNIS KRISTENSEN ABSTRACT. Given a sample from a fully specified parametric model, let Znbe a given finite-dimensional statistic - for example, an initial estimator or a set of sample moments. We propose to (re-)estimate the parameters of the model by maximizing the likelihood of Zn. We call this the maximum indirect likelihood (MIL) estimator. We also propose a computationally tractable Bayesian version of the estimator which we refer to as a Bayesian Indirect Likelihood (BIL) estimator. In most cases, the density of the statistic will be of unknown form, and we develop simulated versions of the MIL and BIL estimators. We show that the indirect likelihood estimators are consistent and asymptotically normally distributed, with the same asymptotic variance as that of the corresponding efficient twostep GMM estimator based on the same statistic. However, our likelihood-based estimators, by taking into account the full finite-sample distribution of the statistic, are higher order efficient relative to GMM-type estimators. Furthermore, in many cases they enjoy a bias reduction property similar to that of the indirect inference estimator. Monte Carlo results for a number of applications including dynamic and nonlinear panel data models, a structural auction model and two DSGE models show that the proposed estimators indeed have attractive finite sample properties. Keywords: indirect inference; maximum-likelihood; simulation-based methods; bias correction; Bayesian estimation. JEL codes: C13, C14, C15, C33. Date: May 2011. We wish to thank M. Arellano, S. Bonhomme, C. Bos, F. Crudu, U. Müller, P.C.B. Phillips, E. Sentana and participants at seminars at Columbia University, Groningen University, Singapore Management University and at the Greater New York Area Econometrics Colloquium 2010 at NYU for helpful comments and suggestions. This work was supported by grants MICINN-ECO2009-11857, SGR2009-578, and NSF grant no. SES-0961596. 1
INDIRECT LIKELIHOOD INFERENCE 2 1. INTRODUCTION Suppose we have a fully specified and thus simulable model, indexed by a parameter θ∈Θ⊂Rk. We have observed a sample Yn=(y1, ..., yn)generated at the unknown true parameter value θ0about which we wish to learn. A natural tool to this end is the likelihood function, f(Yn|θ), and the associated maximum likelihood estimator (MLE), which has a number of attractive large sample optimality properties. However, the MLE is in some situations difficult to compute due to the complexity of the model, and it may require numerical approximations that can deteriorate the performance of the resulting approximate MLE. For example, if the model involves latent variables, they must be integrated out in order to obtain the likelihood in terms of observables. Moreover, even if the MLE is easily computed, it may suffer from significant biases in finite samples with the resulting precision being rather poor, which complicates finite-sample inference. Wellknown examples are the biases of least-squares estimators in autoregressive models (Andrews, 1993) and in dynamic and nonlinear panel data models (Hahn and Kuersteiner, 2002; Hahn and Newey, 2004). To deal with the issue of computational complexity, researchers often resort to GMMtype methods where a statistic Zn=Zn(Yn)is used to draw inference regarding the parameter of interest. Suppose for example, that Znis a set of sample moments: Then a natural way to estimate parameters is to minimize the distance between sample and model-implied moments. When the form of the population moments are unknown, simulations may be used, and one obtains the simulated method of moments (SMM; McFadden, 1989; Duffie and Singleton, 1993). The indirect inference estimator (II; Gouriéroux, Monfort, Renault, 1993; Smith, 1993) proposes an alternative choice for Zn, namely as an extremum estimator based on an auxiliary model. The efficient method of moments (EMM; Gallant and Tauchen, 1996) sets Znto be the score vector of an auxiliary model. Similarly, there exist numerous methods designed to reduce biases in estimators such as bootstrap (Everaert and Pozzi, 2007; Hall and Horowitz, 1996), jackknife (Hahn and Newey, 2004; Kezdi, Hahn, and Solon, 2001), analytical methods (Hahn and Kuersteiner, 2002; Hahn and Newey, 2004) and II (Gouriéroux, Phillips and Yu, 2010; Gouriéroux, Renault and Touzi, 2000). Alternatively, one can adjust the estimator to obtain medianunbiased estimators; see e.g. Andrews (1993). With Znchosen as the initial estimator, one can think of these methods as a type of GMM procedure where the sample statistic is matched against its model implied version, e.g. its finite-sample mean or median to obtain a new, improved estimator. This is in particular the case with the II estimator when the auxiliary model is chosen as the actual model. We here propose a method that offer finite-sample improvements over the aforementioned estimation methods. As with all the above GMM-type estimators1, we take as starting point some statistic Zn, which, for example, could be an initial estimator of θ0, a set of sample moments, or an auxiliary model statistic as used in II. However, rather than minimizing some L2-distance, we propose to (re-)estimate the parameters of interest by 1We use the term “GMM-type estimators” to refer to GMM, MSM, II or EMM estimators based upon a statistic Zn, as described in the text.
INDIRECT LIKELIHOOD INFERENCE 3 maximizing the likelihood implied by Zn. This leads to a maximum-likelihood type estimator which we call the maximum indirect likelihood estimator (MIL), since we operate though the statistic rather than on the sample directly. As a computationally attractive alternative to the MIL, we also propose a Bayesian version of our estimator which is termed a Bayesian indirect likelihood (BIL) estimator. These IL estimators offer finite-sample improvements over the corresponding GMM-type estimators based on the same statistic as we will argue in the following. We derive the asymptotic distributions of the IL estimators and find that they are firstorder equivalent to the GMM estimator that is based on the same auxiliary statistic and uses an optimal weighting matrix. However, the former will in general enjoy better small sample performance compared to the latter for two reasons: First, while GMM estimators only utilize the first and second moment of the statistic, IL estimators are based on a full description of its finite sample distributional characteristics. As such, we expect them to be superior to the GMM estimator in terms of higher-order optimality criteria such as the “large deviations” principle (Bahadur, Zabell and Gupta, 1980), higher-order efficiency (Pfanzagl and Wefelmeyer, 1978), and large deviation probabilities of type II errors (Zeitouni and Gutman, 1991). Second, the first-order equivalence results rely on the GMM estimator being computed using the optimal weighting matrix. Since this in general is unknown, it has to be estimated in order for the efficient GMM estimator to be feasible. This is particularly difficult in time series models where HAC-type estimators have to be employed. In contrast, for our estimators there is no need to estimate an optimal weighting matrix since the likelihood function already embodies the information inherent in the optimal weight matrix. This eliminates an important source of imprecision that can adversely affect the small sample performance of over identified GMM-type estimators (Altonji and Segal, 1996; Doran and Schmidt, 2006; Hansen, Heaton and Yaron, 1996). To justify the above claims of higher-order optimality of the MIL over the corresponding GMM estimator, we provide a higher-order asymptotic analysis of both estimators. In particular, we demonstrate that while the competing estimators have same leading variance components, and so are first-order equivalent, the MIL estimator is third-order efficient in the sense that it has a smaller higher-order variance relative to the GMM estimators. The implementation of the indirect likelihood estimators depends on the likelihood function of the statistic being available on closed form, which will normally not be the case. However, if the model is fully specified, our ability to learn about the likelihood of the statistic is limited only by willingness to do simulations. In particular, we formulate feasible versions of the MIL and BIL estimators by combining simulations with nonparametric density and regression techniques respectively as in, for example, Creel and Kristensen (2009), Fermanian and Salanié (2004), and Kristensen and Shin (2008). The simulated versions are shown to be asymptotically first-order equivalent to the infeasible MIL and BIL estimators as the number of simulations increases. The above mentioned theoretical arguments for improved finite-sample performance of our indirect likelihood estimators over GMM-type estimators are supported by Monte
INDIRECT LIKELIHOOD INFERENCE 4 Carlo results. We investigate the performance of the proposed estimators using a wide range of models, including time series, dynamic and nonlinear panel data, structural auction, and dynamic stochastic general equilibrium models. In terms of root mean squared error and bias, we find that the simulated version of the BIL estimator exhibits performance that is almost always as good, and in most cases better, than the corresponding GMM-type estimators. In particular, BIL is found to inherit the automated bias-correction feature of the standard Indirect Inference (II henceforth) estimators discussed above. When this paper was nearly completed, we became aware of so-called Approximate Bayesian Computation (ABC) or likelihood-free Bayesian inference (see, e.g., Tavaré et al., 1997; Marjoram et al., 2003; Sisson, Fan and Tanaka, 2007) which are used in the biological sciences, including genetics, epidemiology and population biology. One form of ABC (Beaumont, Zhang and Balding, 2002) directly implements what we call the simulated BIL (SBIL) estimator. While the ABC literature is quite mature from an empirical point of view, no theoretical results are available for ABC estimators and their simulated versions, and so this paper offers a number of contributions in this direction. Moreover, the ABC literature only contains rather limited results on the BIL’s finite-sample performance; we provide extensive Monte Carlo examples investigating this. As such, this paper provides an asymptotic theory and finite-sample analysis that has been missing for this literature. The remains of the paper is organized as follows: Section 2 presents the indirect likelihood estimators, and Section 3 discusses their implementation. Firstand higher-order theory of the estimators are developed in Sections 4 and 5 respectively. Section 6 contains the simulation studies, while Section 7 concludes. All proofs have been relegated to the Appendix. 2. INDIRECT LIKELIHOOD INFERENCE We consider the setting described in the introduction, where we wish to learn about a parameter θ∈Θ⊂Rkdescribing a model. Given a sample Yn=(y1, ..., yn)from the model, we choose to make inference on θthrough a d-dimensional statistic of the sample, Zn=Zn(Yn)∈Rd. We can think of Ynas a (random) mapping taking a parameter value into the corresponding observed sample, Yn=Yn(θ). This in turn implies that the statistic also implicitly is a function of θthrough the data, and write Zn(θ)≡Zn(Yn(θ)). In particular, the observed statistic is this random mapping evaluated at the true parameter value which we denote θ0,Zn=Zn(θ0). Let fn(Zn|θ)be the likelihood of the statistic for a given value of the parameter. Suppose for now that the likelihood of the statistic is known on closed form.2We then propose to estimate the parameters by maximizing the indirect likelihood defined through Zn: (1) ˆ θMIL =arg sup θ∈Θ log fn(Zn|θ). 2In general, this will not be the case; in the next section we therefore develop a simulated version of it.
INDIRECT LIKELIHOOD INFERENCE 5 The proposed estimator is indirect, because the sample data is filtered through a statistic, and we refer to the estimator as a maximum-indirect likelihood (MIL) estimator. Compared to the actual MLE based on the full sample, the MIL estimator will in general suffer from an information loss and will only obtain full maximum-likelihood efficiency if the statistic is sufficient in the sense that it spans the score of the full sample log-likelihood. On the other hand, the computation of the indirect likelihood is a lowerdimensional problem compared to the full likelihood (dim(Zn)<dim (Yn)). Moreover, even when the full MLE is computationally feasible, the IL estimator can be used to adjust for finite-sample biases as argued below. Finally, we note that the IL estimator in general will be more robust compared to the full MLE in that it can handle misspecified models and remains consistent as long as Znidentifies the parameter of interest. These are to some extent shared by GMM estimators based on the same statistic. However, in finite samples the two estimators will perform differently, and the MIL will in general exhibit higher-order improvements relative to the GMM estimator. Two leading examples illustrating this general phenomenon are the following: In the first example, suppose that we have available some initial estimator, say ˆ θ. Under suitable regularity conditions, this estimator will be asymptotically normally distributed centered around the true parameter value θ0. However, in finite samples the estimator will in general not be normally distributed and not be centered around θ0. It therefore appears sensible to try to learn about the estimator’s finite-sample distribution, and utilize this information to obtain a better estimate. By choosing our statistic as Zn=ˆ θ, the MIL estimator is an updated version of the initial estimator that takes into account the finite-sample characteristics of ˆ θ. In particular, we expect that the MIL estimator automatically adjusts for potential biases in the initial estimator. As such it is similar to the II bias adjustment mechanism reported in Gouriéroux, Renault and Touzi (2000) and Gouriéroux, Phillips and Yu (2010). However, since MIL estimator at the same time takes into account features of the distribution of ˆ θbeyond its first moment, it should be expected that it will in general dominate the II estimator. As a second example, suppose Znhas been chosen as a set of sample moments; this is for example the case with simulated method of moments. These are mean-unbiased estimators of the corresponding population means and so there is no need for bias adjustment. As such it would seem that a (two-step) GMM estimator based on Znwould suffice. However, the statistic may in finite samples still be non-Normally distributed and taking into account these features will improve the estimator. Furthermore, in the over identified case where d>k, the efficient GMM requires either knowledge or a preliminary estimator of the efficient weighting matrix. In contrast, the MIL automatically incorporates information about the efficient weight and as such is similar to the (generalized) empirical likelihood (GEL) estimator in that it utilizes the full distributional characteristics of the chosen statistic in the estimation of the parameters. As a consequence, the MIL estimator will share the higher-order optimality properties of the GEL (see Newey and Smith, 2004) and dominate the corresponding GMM estimator. In certain situations, the optimization problem defining the MIL estimator may be difficult to solve numerically. The likelihood function θ7→ fn(Zn|θ)may be non convex,
INDIRECT LIKELIHOOD INFERENCE 6 have multiple local maxima, flat spots, or discontinuities, in which case the global maximizer, ˆ θMIL, can be difficult to compute in practice. These features may be even more pronounced when the estimation is based on a simulated version of the likelihood function. This is particularly an issue when the parameter space Θis “large” since the search has to be done over a large-dimensional space. To circumvent these potential problems in the computation of ˆ θMIL, we introduce a Bayesian version of it as a computationally attractive alternative, since it does not require numerical optimization. In the simulation studies, we focus on the posterior mean of θgiven Zndefined as (2) ˆ θBIL =ZΘθfn(θ|Zn)dθ, where fn(θ|Zn)is the posterior distribution given by fn(θ|Zn):=fn(Zn,θ) fn(Zn)=fn(Zn|θ)π(θ) RΘfn(Zn|θ)π(θ)dθ for some density π(θ)on the parameter space Θ. We refer to this particular estimator as the Bayesian indirect likelihood (BIL) estimator. More generally, θ0could be estimated by: (3) ˆ θBIL =arg inf ζ∈ΘZΘρ√n(θ−ζ)fn(θ|Zn)dθ, for some penalty or loss function ρ(u). This includes the posterior mean which is obtained by specifying a quadratic loss ρ(u)=|u|2, while the τth quantile of the posterior follows from choosing the penalty function as the so-called “check” function, ρ(u)= ∑k i=1(τi−1{ui≤0}), where 1{•}denotes the indicator function. The posterior quantiles can be used to construct asymptotically valid confidence intervals as shown in the next section. It should be stressed that we do not give the BIL estimator a Bayesian interpretation and merely see it as a computational device to circumvent the numerical issues related to the maximization problem that has to be solved in order to compute ˆ θMIL. In particular, we do not interpret π(θ)as a prior density in the Bayesian sense, in that it does not necessarily reflect beliefs about the parameter. It is simply used to give weights to different parts of the parameter space, and in our examples, we alway use a uniform density. As such, ˆ θBIL is close in spirit to the class of Laplace type estimators (LTE’s) introduced in Chernozhukov and Hong (2003). 3. COMPUTATION OF FEASIBLE ESTIMATORS In most situations, it will not be possible to derive the exact finite-sample distribution of the statistic Znon closed form. Thus the likelihood fn(Zn|θ)will normally not be available, and one has to resort to numerical approximations instead. In the computation of the BIL estimator, it is in addition required to compute the integral RΘρn(θ−ζ)fn(θ|Zn)dθ. For the latter problem, one could follow the suggestions of Chernozhukov and Hong (2003) and compute the integral using Markov chain Monte Carlo (MCMC) methods.
INDIRECT LIKELIHOOD INFERENCE 7 However, we here opt for an alternative solution which handles the numerical approximation of fn(Zn|θ)and the integral in one step; the proposed method which we describe below is easy to implement and in general quite robust. First, for the implementation of the MIL, we have to be able to compute fn(Zn|θ)at any given trial value θ. Since the model is simulable and the mapping Zn(θ)≡Zn(Yn(θ)) is known (as chosen by the econometrician), we propose to estimate the density using kernel density methods: Draw Sindependent samples, Ys n(θ)for s=1, ..., S, from the model evaluated at the trial value θ, compute the associated statistic, Zs n(θ)≡Zn(Ys n(θ)), s=1, ...., S, and then estimate the density by kernel methods (see e.g. Li and Racine, 2007, Ch. 1 for an introduction): (4) ˆ fn,S(Zn|θ) = S ∑ s=1 Kh(Zs n(θ)−Zn), where Kh(z)=K(z/h)/h,K(z)is a kernel function and h>0 is a bandwidth. One then embeds the approximated density inside (1), and uses an optimization algorithm to obtain an estimator. This yields a simulated MIL (SMIL) estimator: (5) ˆ θSMIL =argsup θ∈Θ log ˆ fn,S(Zn|θ). The simulated version is akin to the nonparametric simulated maximum-likelihood estimator (NPSMLE) of Fermanian and Salanié (2004) and Kristensen and Shin (2008). The above kernel density estimator implicitly assumes that Zn(θ)has a continuous distribution. However, we show that even if this is not the case, the simulated version will still asymptotically behave as the MIL estimator. For the computation of the BIL estimator, we not only need to evaluate the likelihood but also the integral over the quasi-posterior density. Chernozhukov and Hong (2003) propose to handle the latter computational problem through MCMC, but this can be quite a delicate method which in some cases has unstable properties (see Kormiltsina and Nekipelov, 2009). Instead, we opt to also combine simulations and nonparametric techniques in the implementation of the BIL estimator. Suppose, to illustrate, that the penalty function is ρ(u) = |u|2. In this case, the Laplace-type estimator is the mean of the posterior density, ˆ θBIL =ZΘθfn(θ|Zn)dθ=E[θ|Zn]. Our idea is then to compute ˆ θBIL =E[θ|Zn]by combining nonparametric regression methods and simulations as follows: Make i.i.d. draws θs,s=1, ..., S, from the pseudoprior density π(θ), for each draw generate a sample Yn(θs)from the model at this parameter value, and then compute the corresponding statistic Zs n=Z(Yn(θs)),s=1, ..., S. Given the i.i.d. draws (θs,Zs n),s=1,...S, we can obtain a simulated version of the BIL (SBIL) through nonparametric regression techniques. One such is the kernel estimator (see Li and Racine, 2007, Ch. 2), (6) ˆ θSBIL =∑S s=1θsKh(Zs n−Zn) ∑S s=1Kh(Zs n−Zn),
INDIRECT LIKELIHOOD INFERENCE 8 while another one is the k-nearest neighbor (KNN) estimator (see Li and Racine, 2007, Ch. 14), where the bandwidth is chosen as h=dk(Zn)with dk(Zn)denoting the Euclidean distance between Znand the k-th nearest neighbor among the simulated values. As such the KNN estimator can be thought of as a kernel regression estimator with an adaptive bandwidth. A member in the general class of BIL estimators given in eq. (3) can be expressed as minimizing a conditional moment, ˆ θBIL =arg inf ζ∈ΘEρ√n(θ−ζ)|Zn, which can be approximated by replacing the exact moment by a simulated nonparametric version, ˆ θSBIL =arg inf ζ∈Θ ˆ ESρ√n(θ−ζ)|Zn, where ˆ ESρ√n(θ−ζ)|Zn, for example, can be computed by kernel regression, (7) ˆ ESρ√n(θ−ζ)|Zn=∑S s=1ρ√n(θ−ζ)Kh(Zs n−Zn) ∑S s=1Kh(Zs n−Zn), or nearest neighbor estimation where again h=dk(Zn). When the penalty function is chosen as the “check”-function, this leads to simulated versions of the posterior quantiles, which are used to compute confidence intervals. In this case, the kernel smoothed version becomes the kernel quantile regression estimator (Li and Racine, 2007, Sec. 6.4). For both the SMIL and SBIL, there are two sources of error in comparison with the exact MIL and BIL estimators (which only suffer from the sampling error in Zn). First, randomness is added due to the use of simulations, and there is also a bias component due to the use of nonparametric estimators. We treat the nonparametric fitting step as a computational tool used to find the value of the estimator, in the same way that Chernozhukov and Hong (2003) treat MCMC as a means of computing LTEs. As the number of simulated draws Sbecomes large, nonparametric density and regression estimators are consistent. Thus, both the randomness due to use of simulations and the bias due to use of nonparametric methods can be controlled for by choosing Ssufficiently large. We analyze the impact of simulations and kernel smoothing in Section 6. One may wish to explore different pseudo-priors. If a large body of simulations have been generated using the pseudo-prior π(θ), then one can obtain results for a different pseudo-prior without doing additional simulations by using importance sampling. The SBIL estimator based on π(θ) presented in equation 6can be written as ˆ θSBIL = ∑S s=1θswS(Zs n,Zn), where wS(Zs n,Zn)has an obvious definition. Given simulations {(θs,Zs n)}S s=1 based on π(θ), the SBIL estimator corresponding to the new pseudo-prior, say π∗(θ), can be computed by ˆ θSBIL = S ∑ s=1 θswS(Zs n,Zn)π∗(θs) π(θs). This may be useful when Sis very large or when it is costly to compute the auxiliary statistic, as in the case of the DSGE models presented later in this paper.
INDIRECT LIKELIHOOD INFERENCE 9 In the ABC literature or likelihood-free literature, discussed in the introduction, methods of computing estimators using likelihood-free Markov chain Monte Carlo and sequential Monte Carlo have been studied in some detail (Marjoram et al., 2003; Sisson, Fan and Tanaka, 2007; Beaumont et al. 2009). One could also employ so-called importance sampling to reduce variances due to simulations: For any conditional density gn(θ|z) with support Θ, we can rewrite ˆ θBIL as ˆ θBIL =Rθ{fn(θ|Zn)/gn(θ|Zn)}gn(θ|Zn)dθ, and so a generalized version of our proposed simulated version would be ˆ θBIL =1 S S ∑ s=1 θsˆ fn,S(θs|Zn) gn(θs|Zn)=∑S s=1θsπ(θs)/gn(θs|Zn)Kh(Zs n−Zn) ∑S s=1Kh(Zs n−Zn), where θs∼i.i.d.gn(θs|Zn). The optimal choice of gn(θ|z)in terms of variance reduction is gn(θ|z)=1{θ∈Θ}|θ|fn(θ|z) RΘ|θ|fn(θ|z)dθ. Unfortunately, it is not feasible to draw from this choice since RΘ|θ|fn(θ|z)dθis unknown, but approximate methods exist; see, for example, Zhang (1996). These methods should in principle be computationally more efficient compared to the basic sampling method in equation (6) to obtain a given level of precision. In our simulation study we focus on the basic sampler, and leave the implementation of importance samplers for future research.Application of these methods could provide savings in computational time when it is costly to sample from the model, but on the other hand require more careful implementation. Of the examples considered in this paper, only the DSGE models (below) present serious computational burden. 4. FIRST-ORDER ASYMPTOTICS As a first step towards a complete asymptotic analysis of the MIL and BIL estimators, we here derive their first-order asymptotic distribution. The asymptotic analysis of the MIL estimator proceeds along the standard steps for parametric extremum estimators, while the BIL estimator on the other hand requires a bit more care. Fortunately, since the BIL estimator can be regarded as a specific LTE, we can employ the general results of Chernozhukov and Hong (2003) to establish √n-consistency and asymptotic normality of our Bayesian estimator, as well as the equivalence with the MIL estimator when the penalty function ρis symmetric. We impose the following conditions on the parameter space and the weighting function Assumption 1. Assume that: (i) the parameter space Θ⊂Rkis compact with θ0being an interior point; (ii) the weighting function π(θ)is a continuous, uniformly positive density; and (iii) the penalty function is convex and satisfies ρ(u)=0⇔u=0, ρ(u)≤1+|u|pfor some p≥1, and φ(x)=Rρ(u−x)eu0audu is uniquely minimized at some x∗for any a >0. This set of assumptions is completely standard, and are identical to the conditions found in Chernozhukov and Hong (2003). It should be noted that (ii)-(iii) are only needed to develop theory for the BIL estimator, and the asymptotics of the MIL estimator only require (i).
INDIRECT LIKELIHOOD INFERENCE 16 we demonstrate that the bias properties of the MIL are very favorable and are as good as the ones of the CU GMM-type estimator. We now turn our attention to the higher-order efficiency of the GMM and IL estimators. As was shown in the previous section, their first-order asymptotic variances are identical. However, in finite samples, the IL estimators are expected to dominate for a number of reasons: First, the GMM estimators are only first-order equivalent to the IL estimators if Wn=Ω−1(θ0)+op(1). If not, the IL estimators are asymptotically more efficient than GMM. Moreover, the first-step estimation error contained in Wnin general has an adverse impact on the performance of the resulting two-step estimator which may perform poorly in small and moderate samples; see e.g. Altonji and Segal (1996), Hansen, Heaton and Yaron (1996) and Newey and Smith (2004). Furthermore, while increasing the dimension of the auxiliary statistic increases the asymptotic efficiency of the GMM estimator, it also increases the dimension of the weight matrix to be estimated and numerical singularities can appear making the inversion difficult. In contrast, our estimators do not require estimation of the optimal weighting matrix, and so increasing the dimension of the auxiliary statistic causes no difficulties with singular matrices. If on the other hand Ω(θ)is known, then we can estimate the parameters using the CU estimator which will remove the additional estimation errors due to the use of Wn; see Donald and Newey (2000) and Newey and Smith (2004). However, in finite samples, the CUE still only utilizes information contained in the first and second moments of Zn, while the MIL takes into account all distributional characteristics. This difference means that the indirect likelihood estimators in general will have better small sample performance than both two-step efficient GMM and CU based on the same auxiliary statistic. The formal proof of higher-order efficiency can be done by ranking the GMM-type and MIL estimators in terms of their higher-order MSE. If an estimator ˆ θsatisfies the expansion in eq. (14), we obtain (again ignoring Rn) MSE √n(ˆ θ−θ0)≃BnB0 n+Vn, where Bn=√nBias ˆ θand Vn=nVar ˆ θ. For each of the three estimators, the variance can be decomposed into Vn=J−1+Ξ/n+o(1/n), where J−1=J−1(θ0)is the leading variance component, while Ξis the higher-order variance. The expression of Ξfor each of the three estimators (CUE, GMM, MIL) is straightforward to obtain from the expansion, but it is rather complicated. This makes a direct ranking of the estimators in terms of their respective Ξ’s difficult. Instead, we first develop an Edgeworth expansion of the distribution of the MIL estimator. For standard maximum-likelihood estimators where the log-likelihood takes the form of a sample average over i.i.d. observations, Edgeworth expansions have been established; see, for example, Bhattacharya and Ghosh (1978). However, we can in general not write log fn(Zn|θ)as a sample average of i.i.d. variables and so the standard proof does not directly carry over to our setting. However, by importing some of the arguments of Bhattacharya and Ghosh (1978), we can still show that ˆ θMIL ≃H(Wn(Zn)) for some analytic function Hand with Wn(Zn)denoting the first rderivatives of log fn(Zn|θ)w.r.t.
INDIRECT LIKELIHOOD INFERENCE 17 θ. Since (the normalized version of) Znsatisfies an Edgeworth expansion, we can then apply the general results of Phillips (1977) on Edgeworth expansions of transformations of random sequences to obtain the desired result: Proposition 5. Under Assumptions 1-4with E |Zn|p<∞for all n,p≥1, and P|Zn−Z(θ0)|>c1qlog (n)/n=on−r/2, the MIL satisfies an rth order Edgeworth expansion: sup yP√nˆ θMIL −θ0J(θ0)≤y− y Z −∞ φ(x)"1+ r ∑ i=1 n−i/2 ˜ πi(x)#dx=on−r/2, where ˜ πi(x)is a polynomial of order 3i, i =1, ...,r. The assumption that Znhas moments of all orders is somewhat restrictive and rules out heavy tails. We conjecture that this assumption is not strictly necessary for the above result to holds. In particular, one might be able to show Proposition 5by using the results of Skovgaard (1981) where weaker moment restrictions are required. This would on the other hand complicate the proof and so for clarity we maintain the assumption of all moments existing. The tail probability condition is satisfied for most regular statistics; see, e.g., Bhattacharya and Ghosh (1978, Theorem 3). Once we have shown that the distribution of the MIL estimator can be approximated by an Edgeworth expansion, it now follows by standard results for maximum-likelihood estimators (see e.g. Ghosh, 1994 and Bickel, Götze, and van Zwet, 1985), that the biasadjusted MIL estimator is third-order efficient amongst all estimators relying on the statistic Zn. In particular, ΞGMM ≥ΞMIL and ΞCUE ≥ΞMIL. 6. PROPERTIES OF SIMULATED VERSIONS We analyze the impact of the use of simulations and nonparametric estimation in the implementation of the MIL and BIL estimators. For the simulated version of the MIL estimator, we combine the general results of Kristensen (2009) and Kristensen and Shin (2008) to show that it is first-order asymptotically equivalent to the infeasible MIL estimator. The analysis of the simulated version of the BIL estimator can be done directly since we can write it up on closed form. To utilize existing results on convergence rates of kernel estimators, we make the following assumptions regarding the kernel function used in the computation of the SMIL and SBIL defined in Section 3: Assumption 5. The kernel K satisfies: There exist C,L<∞such that either (i) K(u) = 0for kuk>L and |K(u)−K(u´)| ≤ Cku−u0k,or (ii) K(u)is differentiable with supu|K0(u)|< ∞. For some a >1,|K(u)| ≤ Ckuk−afor ||kuk>L„ and RK(z)dz =1,RzK (z)dz =0, Rz2K(z)dz <∞. The above assumptions imposed on the kernel are quite standard and are for example satisfied by the Gaussian kernel. We first restrict ourselves to the case where the likelihood is a density:
INDIRECT LIKELIHOOD INFERENCE 18 Assumption 6. The indirect likelihood fn(z|θ)is a density with respect to the Lebesgue measure and is twice continuously differentiable in z. Under this assumption on the kernel and the likelihood, the following result holds: Proposition 6. Assume that Assumptions 1-6hold. Then the SMIL and SBIL estimators are asymptotically first-order equivalent to the actual ones under the following conditions: For the kernel-smoothed versions, nh2→0and nlog (S)/Shd→0. For the nearest-neighbor versions, n [k/S]2/d→0, and nlog (S)/k→0. The restrictions on Sand hare fairly standard and require the number of simulations to grow at a slightly faster rate than the number of observations. In particular, standard bandwidth selectors will satisfy the above rates and so these can be used in the implementation of the simulated versions. We also note that the simulated versions of our estimators suffer from a curse of dimensionality. This appears explicitly in the conditions on Sand hgiven in Proposition 6 where we require nlog (S)/Shd→0. Thus, the larger d=dim (Zn)(which must be at least that of θ, and which is larger in most of the applications below), the more simulations are required for the simulations to have a negligible impact on the estimator. This is a well-known issue which is shared by most other simulation-based estimators: The larger the dimension of the space over which we need to integrate, the larger the number of simulations should be chosen to control the simulation error. The above result requires the likelihood to be a density. We now demonstrate that the SMIL and SBIL estimators enjoy the same asymptotic properties even if this is not the case. In fact, we will not even require that Assumption 4holds and as such allow for both continuous and discrete observations. To be more specific, we replace Assumptions 4and 6with the following one: Assumption 7. For some N ≥1: supn≥NEhsupθ∈Θ√n(Zn(θ)−Z(θ))2i<∞. This uniform integrability assumption is satisfied if, for example, Zn(θ)is a sample average with second moment. It imposes no smoothness restrictions on the finitesample likelihood and holds for both continuous and discrete underlying data. It is used in conjunction with Assumption 2to ensure that √nE [|Zn(θ)−Z∗ n(θ)|]→0, where Z∗ n(θ)∼N(Z(θ),Ω(θ)/n)is its Normal limit sequence. We use this to show that the kernel smoother based on simulations from the distribution of Zn(θ)converges towards the one based on simulations of Z∗ n(θ). Since Z∗ n(θ)satisfies Assumption 4and 6by construction, this in turn implies that SMIL and SBIL have the desired asymptotic properties: Proposition 7. Assume that Assumptions 1-3,5and 7hold, and the kernel K is uniformly Lipschtiz, |K(u)−K(v)|≤D|u−v|. Then the SMIL and SBIL have the same asymptotic properties as those stated in Proposition 1under the bandwidth conditions stated in Proposition 6 together with nh2→∞(kernel smoother) and n/k2→∞(nearest neighbor). The intuition behind the above result is the following: If the distribution of Zn(θ)cannot be described by a density, one can think of the kernel smoothing inherent in both the SMIL and SBIL as a type of regularization that generates a smooth objective function
INDIRECT LIKELIHOOD INFERENCE 19 which can be used instead of the more irregularly behaved true likelihood. As such the SMIL and SBIL estimators are similar in nature to the smoothed maximum score estimator proposed in Horowitz (1992) where a non-smooth estimator is regularized through smoothing. In practice, we choose the number of simulations Sso large, that the additional variance due to simulations is negligible. However, for completeness, we note that the simulated version of the BIL estimator satisfies ˆ θSBIL =ˆ θBIL +ES(Zn), for a stochastic function ES(z)which is independent of ˆ θBIL and satisfies either (in the case of kernel-smoothers), √ShdES(z)→dN0, kKk2σ2 n(z) fn(z), or (in the case of nearest-neighbor estimators), √kES(z)→dN0, kKk2σ2 n(z), where d=dim (Zn),kKk2=RK2(z)dz, and σ2 n(z)=Var [θ|Zn=z]. Thus, the variance estimator of the kernel-smoothed version of SBIL could be adjusted by adding kKk2σ2 n(z) fn(z)/Shdto J−1(θ0), and similarly for the nearest-neighbor version. A similar adjustment can be developed for the MIL estimator by using the arguments of Kristensen and Salanié (2010). 7. MONTE CARLO RESULTS In this section we explore the performance of the SMIL and SBIL estimators, comparing them to other estimators, using a variety of econometric models including simple time series models, a dynamic and nonlinear panel data models, a structural econometric model of an auction and two dynamic stochastic general equilibrium (DSGE) models. We focus on several issues. First, the SBIL estimator is considerably more convenient to use than is the SMIL estimator, from a computational point of view, so we would like to know if the two estimators perform similarly before focusing our attention on the SBIL estimator. Second, Proposition 5tells us that the exact MIL is higher-order more efficient than the GMM estimator that uses the optimal weight matrix. This leads us to hope that the SMIL and SBIL estimators have better small sample performance than GMM-type competitors. A factor that could undermine these potential gains is the need to use simulations and nonparametric fitting to implement the feasible versions (the feasible SMIL and SBIL versus the infeasible MIL and BIL). This section throws light on the actual performance of the feasible versions. A third purpose of this section is simply to give examples of how the SMIL and SBIL estimators may be implemented in practice. Examples of practical issues to deal with are the choice of the auxiliary statistic, and the specification of the parameter space in the case of the SBIL estimator. A fourth issue is the accuracy of confidence intervals computed using estimated quantiles of the pseudo-posterior. We find mixed results for confidence interval coverage: in
INDIRECT LIKELIHOOD INFERENCE 20 some cases coverage is very accurate, while in others the confidence intervals are too broad, so true size is smaller than the nominal size. Because the findings are mixed, we do not present tabular results, and we leave this issue for future research. It is perfectly feasible to use other means (asymptotic, bootstrap, Monte Carlo) of computing confidence intervals and standard errors for the SBIL estimator. For example, Li (2010) found that bootstrap confidence intervals are very accurate for the II estimator of the structural auction model discussed below. The same method could be used for the SBIL estimator. We do not pursue the issue further in this paper. To implement the SMIL and SBIL estimators, we use between S=106and S=107 simulated points drawn randomly from the parameter space, depending on the application. The auxiliary statistics we use are in most cases computationally inexpensive, so generating a large number of replications is not burdensome. The exceptions are the DSGE models, which requires approximately two days of time on a 32 core cluster per 106replications of the auxiliary statistic3. For all problems, we use at least 5000 Monte Carlo replications at each design point. The nonparametric fit is done using the knearest neighbors approach4, using the ANN library (Arya, Malamatos and Mount, 2009; http://www.cs.umd.edu/~mount/ANN/). Using this C++ library, the KNN nonparametric fitting step requires at most several minutes of time on a single core. It is also a simple matter to switch to using approximate nearest neighbors, which can speed up the nonparametric fitting step if one uses an extremely large number of simulations. The number of neighbors kused for the nonparametric fit is chosen (with one exception) as k=1.5 ×S0.25, rounded down to the nearest integer. More careful choice using methods such as cross validation might improve the results, but we do not explore this possibility in this paper. We report the SBIL estimator computed as the posterior mean. The version computed as the posterior median gives very similar results. For all applications the pseudo-prior π(θ)is a uniform distribution over the parameter space Θ, so the only remaining issue is specifying the bounds of parameter space. For some of the applications (the MA and dynamic panel data models), prior beliefs such as stationarity or invertibility lead directly to the specification of at least some of the bounds of parameter space. For others (the auction model and the DSGE models) we have less information available regarding plausible bounds on at least some of the parameters. The issue of setting the parameter space in such cases is addressed in the subsection presenting the auction model. 7.1. Dynamic panel data. Gouriéroux, Phillips and Yu (2010; henceforth GPY) investigate the performance of the II estimator using a linear dynamic panel model (17) yit =αi+φ0yit−1+eit 3Performing Monte Carlo on a cluster is quite straightforward. We use PelicanHPC (http://pelicanhpc. org/), a framework very similar to that described in Creel (2007). 4We also have used kernel regression, which gives very similar results to the KNN results reported here.
INDIRECT LIKELIHOOD INFERENCE 21 where eit ∼N(0,1),αi∼N(0,1),φ0=0, 0.3, 0.6, 0.9 and αiand eiare independently distributed. The initial condition is yi0|αi∼Nαi 1−φ0,1 1−φ2 0. GPY use the (inconsistent) ML “fixed effects” estimator as the auxiliary statistic. They find that the II estimator outperforms a number of alternative estimators, in terms of root mean squared error (RMSE). While our asymptotic results do not straightforwardly generalize to dynamic panel data models (where the theory normally requires the number of time periods, T, to grow with sample size), we conjecture that the higher-order efficiency results also hold in this context. We here investigate this claim by comparing the performance of the SMIL and SBIL estimators to the II results obtained by GPY. GPY also report results for other bias correction methods such as jackknife and analytical bias correction and find that their II estimator dominates those; we therefore focus on the II estimator and do not reproduce the results for the other estimators. The parameter space is set to the stationary region φ0∈(−1, 1). We consider two auxiliary statistics: the same ML estimator as used by GPY, and also the ML estimator augmented with the OLS estimator of the naive model yit =δyit−1+νit that ignores the presence of individual effects. The SMIL estimator requires a nonparametric density fit embedded inside an optimization problem, while the SBIL estimator eliminates the optimization. In the present case, the parameter to estimate is a scalar, so for this problem it is relatively easy to apply both the SMIL and SBIL estimators. By comparing the two in this relatively simple case, we can get an indication of whether focusing on the SBIL estimator in more computationally demanding cases is warranted by a comparable performance of the two estimators. To implement the SMIL, we use a different approach than what is outlined in equations (4) and (5). The reason for this to take advantage of the large set of replications of (θs,Zs n) that are already available after computing the SBIL estimator. Instead of operating on the conditional density fn(Zn|θ), we work with the joint density fn(Zn,θ). When θsis drawn from a uniform density, as is the case here, fn(Zn,θ)and fn(Zn|θ)are maximized at the same value of θ, because the marginal density of θdoes not depend upon θ. We of course do not know the joint density, so it must be fit nonparametrically. We use the simple KNN density estimator given in equation 14.2 of Li and Racine (2007) to fit fn(Zn,θ). This nonparametric fit to the joint density, ˆ fn(Zn,θ)is then maximized with respect to θ using a grid search, in order to deal with the rough, nondifferentiable nature of the KNN density estimator. Because θis a scalar in the present case, use of grid search does not present a significant computational burden. Table 1presents the bias of the estimators, and Table 2presents the root mean squared errors (RMSEs). In these Tables, the columns labeled II, SBIL and SMIL all refer to use of the auxiliary statistic Zn=b φML, while the columns labeled SBIL(OI) and SMIL(OI) refer to use of the overidentifying auxiliary statistic Zn=b φML,b δOLS. Results for the inconsistent ML estimator are also presented, for reference. We see that the II and SBIL estimators have very small biases in almost all cases. With an exactly identifying auxiliary statistic, the estimators (except ML) all have similar biases and RMSEs, especially for
INDIRECT LIKELIHOOD INFERENCE 22 larger sample sizes. For small sample sizes, the SBIL estimator performs somewhat better than the II estimator, overall. When the difference favors the II estimator, it is small, but when it favors the SBIL estimator, it is larger. For the SMIL and SBIL estimators, it is easy to use an overidentifying auxiliary statistic, because no covariance matrix need be estimated. Looking at the columns labeled SMIL(OI) and SBIL(OI), we see that there are gains from doing so: bias is essentially unchanged, but RMSE is reduced considerably, especially for smaller sample sizes. There seems to be no reason to prefer SMIL to SBIL, as the RMSEs of the two are essentially the same in the case of the exactly identifying auxiliary statistic, while SBIL almost uniformly dominates SMIL when the overidentifying auxiliary statistic is used. Based on the good performance of SBIL compared to SMIL in this example, and the fact that the two estimators are first order equivalent, we focus on SBIL in the remaining examples. Most of the remaining examples have parameter vectors of higher dimension, which would make a global maximization strategy such as grid search or simulated annealing more tedious to employ (recall that a nonparametric density fit must be done for each trial parameter value). The SBIL estimator does not require this optimization step, so it avoids this difficulty. 7.2. Moving average. The previous section compared the proposed estimators to a just identified II estimator. It is also desirable to compare to an overidentified II estimator, because this is the situation where it is necessary to estimate the efficient weight matrix in order to obtain an efficient II estimator, given the chosen auxiliary statistic. We would like to see if the SBIL estimator benefits from the fact that it does not require estimation of the efficient weight matrix. The first order moving average (MA(1)) model has been widely used to investigate the performance of the indirect inference estimator, and a pth-order autoregressive model is often used to generate the auxiliary statistic (see, for example, Gouriéroux,Monfort and Renault, 1993; Chumacero, 2001). In this section we estimate the MA(1) model yt=et+ψet−1 et∼i.i.d.N(0, σ2) using sample sizes of n=50, 100 and 200 observations. The parameter ψis one of the values {−0.95, −0.9, −0.5, 0, 0.5, 0.9, 0.95}, so the model is always invertible. The parameter σis always equal to 1. The parameter vector is θ= (ψ,σ). We set the parameter space to Θ=(−1,1)×(0, 2), which imposes invertibility, which is needed for the parameter to be identified. The statistic Znis the vector of estimated parameters ρ0,ρ1, ..., ρP,σ2 υof an AR(P) model yt=ρ0+∑P p=1ρpyt−p+υt, fit to the data using ordinary least squares. For simplicity, we hold the order of the AR(P)model constant at P=10 across the Monte Carlo replications. Thus, the dimension of Znis 12, while the dimension of θis 2, so we have considerable overidentification. We estimate θusing SBIL and II, where both are based on the auxiliary statistic defined in the last paragraph. The II estimator is computed using continuously updated
INDIRECT LIKELIHOOD INFERENCE 23 GMM (Hanson, Heaton and Yaron, 1996). The moment conditions that define the continuously updated indirect inference (CU-II) estimator are mn(θ) = Zn−¯ ZS,n(θ)where ¯ ZS,n(θ) = 1 S∑S s=1Zs n(θ), and the weight matrix at each iteration is the inverse of ΩS n(θ) = 1 S∑S s=1[Zs n(θ)−¯ ZS,n(θ)] [Zs n(θ)−¯ ZS,n(θ)]0, where S=100. For reference, we also estimate θusing the conditional maximum likelihood estimator (Gaussian MLE with e0set to zero). When a replication of the ML or CU-II estimator lies in the non-invertible part of the parameter space, we use the observationally equivalent invertible parameter value in its place. The need for doing this and the means of doing so are explained by Chumacero (2001). Table 3reports the results. In this Table, SBIL(AR) refers to the SBIL estimator that uses the AR(10) auxiliary statistic, while SBIL(ML) is the SBIL estimator that uses the ML estimator as the auxiliary statistic. We see that the SBIL(AR) and CU-II estimators have biases that are of comparable magnitudes, overall. Comparing RMSEs, the SBIL(AR) estimator performs better than the CU-II estimator, almost uniformly. This result is not unexpected, given the previous theoretical results for higher order efficiency of MIL compared to CU-II. These theoretical grounds for efficiency plus the avoidance of estimation of the weight matrix appear to lead to real small sample efficiency gains. Comparing to the ML estimator, for the smaller sample size, SBIL(AR) has a larger RMSE than does ML, which is no doubt an indication that an AR(10) auxiliary model is excessively parameterized when the sample size is only 50. When the sample size is 200, the SBIL(AR) estimator has bias and RMSE comparable to those of the ML estimator. When the ML estimator is used as the auxiliary statistic for SBIL, there is no benefit in terms of RMSE when the sample size is 50, but for samples of size 100 and 200, the SBIL(ML) estimator has an RMSE lower than that of the ML estimator. 7.3. Nonlinear panel model. Section 7.1 explores a linear panel data model with normally distributed errors. One might expect that a nonlinear model could lead to a larger difference between the SBIL and II estimators, especially for smaller sample sizes, as in such a case the small sample distribution of the auxiliary statistic, which characterizes the objective functions of the IL estimators, could be less well approximated by the corresponding normal limiting distribution, which characterizes the objective function of the II estimator. To investigate this conjecture, we use the static logit panel model that Arellano and Bonhomme (2009) used in some of their Monte Carlo work to compare a set of semi-parametric nonlinear panel data estimators. Their static logit Monte Carlo design (see their Section 7.1) is used here to compare the SBIL and CU-II estimators. The design of the experiment is yit =1[xitφ0+αi0+eit >0] where xit ∼N(0,1)and the individual effects αi0∼N(¯ xi,1), where ¯ xi=1 T∑T t=1xit. The eit are independent draws from the logistic CDF. The true value of φ0=1. We set N∈{30,100}and T=5. The first component of the auxiliary statistic is the estimator of the misspecified logit model that results from the above model, with the exception that, erroneously, it is assumed that the individual effects are all identical. To be precise, it is the quasi-ML estimator resulting from logit estimation of the misspecified model yit =
INDIRECT LIKELIHOOD INFERENCE 24 1[α+xitφ+eit >0]. The second component of the auxiliary statistic is the OLS estimator of the linear probability model yit =α+xitφ+ηit. The logit and OLS estimators of αand φtogether yield an auxiliary statistic of dimension 4, so we have overidentification for the estimator of the scalar φ0. We use SBIL and CU-II to estimate φ0, using this auxiliary statistic. SBIL uses 2 ×106simulations, and CU-II was implemented as described in the previous section. For both SBIL and CU-II, the parameter space for φis set to [0,2]and the pseudo prior for SBIL is a uniform distribution over the parameter space. Table 4presents the results for bias, RMSE and mean absolute error (MAE). For both sample sizes, the SBIL estimator is less biased and has smaller RMSE and MAE than the II estimator. For the smaller sample size, the RMSE of the SBIL estimator is 88.4% that of the CU-II estimator, while for the larger sample size the percentage is 91.6%. This result supports the conjecture that the SBIL estimator will have better small sample performance than that of GMM-type estimators based on the same auxiliary statistic. Comparing these results to those for the linear dynamic panel data model, it seems that the nonlinearity of the model also contributes to accentuate the difference in performance of the SBIL and GMM-type estimators. For the sample size N=100, T=5, the MAE and bias results may be compared with the first panel of Table I in Arellano and Bonhomme (2009). Both SBIL and CU-II have less bias and lower MAE than any of the estimators considered by Arellano and Bonhomme. This is to be expected, because those estimators are semi-parametric, in that the distribution of the individual effects is unknown. The SBIL and CU-II estimators, in contrast, are based on simulations that require knowledge of the distribution of the fixed effects. The assumption that the distribution of the individual effects be known is quite implausible in this example. Nevertheless, the example serves to illustrate how the SBIL and II estimators can achieve a good bias reduction in small samples, through use of a simple naive auxiliary model, when one is able to write a fully simulable model. 7.4. Structural model of an auction. Li (2010) proposes to use indirect inference for estimation of structural econometric models, and illustrates with a Monte Carlo example of estimation of the parameters of a Dutch auction, where only the winning bid is observed. The number of bidders is fixed at N=6, and the sample size is n=100, meaning that the outcomes of 100 auctions are observed. At each auction i=1,2,..., 100, the quality, xi, of the item being auctioned is the square of a uniform (0,2)random variable, to introduce heterogeneity in the values of the objects across the auctions. The 6 bidders draw their independent private values from a common exponential distribution with density f(v|xi) = 1 exp(θ0+θ1xi)exp −v exp(θ0+θ1xi) so that exp(θ0+θ1xi)is the mean valuation of the item, over the bidders. The equilibrium strategy for the winning bid is b∗ i=v∗ i−1 FN−1(v∗ i|xi)Zv∗ i 0FN−1(u|xi)du
INDIRECT LIKELIHOOD INFERENCE 25 where v∗ iis the highest private valuation, and F(·|xi)is the exponential distribution function. For a given value of N(6 in this case), symbolic computation software can be used to obtain an analytic solution for the winning bid, which facilitates simulation of the model. The observed data are the 100 values of {xi,b∗ i}, and we seek to estimate θ0and θ1. The true values are set to θ0=1 and θ1=0.5. Li presents results for indirect inference using two auxiliary statistics: the fitted coefficients of a pseudo ML estimator, and the OLS regression coefficients (b β0,b β1)obtained by fitting the model b∗ i=β0+β1xi+σei. To apply the SBIL estimator, we must specify the parameter space. The present application is interesting, because we have no clear a priori bounds for the two parameters θ0 and θ1. Outside of the Monte Carlo context, one would only have the sample data, but would not know the true parameter value. We discuss the issue of how the parameter space may be specified at some length, because it is a necessary step to apply the SBIL estimator. Our proposal is to start with a parameter space that seems conservatively large, and to check that it in fact contains elements that can generate simulated statistics Zs nthat differ in important respects from the Zngenerated by the sample data. To do this, one can generate a preliminary set of Zs nsetting Ssmall enough to be convenient. Then one may compute the distance between each simulated statistic and the statistic using the sample data, giving the Sdistances ds. Then one can sort the Sreplications of (θs,Zs n,ds)by ds and check that the θsthat generate relatively small distances are always comfortably far away from the bounds of the proposed parameter space. If this is not the case, the parameter space can be expanded, and the procedure repeated again. Conversely, one may find evidence that the proposed parameter space is excessively broad, in that regions of the parameter space never generate statistics close to Zn. Such simulations will not contribute to the nearest neighbors version of SBIL, and as such are wasted. This could be avoided by using importance sampling, but we here for simplicity take a brute force approach and simply choose initially a large parameter space and a moderate number of simulations, Sfor an intial exploration of the distribution of the statistic across different parameter values. We then shrink the parameter space removing parts whith little or no contribution to the posterior distribution. We initially set the parameter space to Θ=(−5,5)×(0, 5). We generate a single sample at the true parameter value, and a fairly small number (105)simulated samples from the proposed parameter space. Inspection of the distribution of the auxiliary statistic used by Li reveals that the auxiliary statistic when sampling from the proposed parameter space presents some extreme outliers. This is a problem that may not be detected when using the II estimator with a limited number of replications of the auxiliary statistics (Li uses only one draw), because the II estimator maintains the underlying random draws fixed over the iterations, to avoid the phenomenon of “chatter” when doing the minimization to compute the estimator. The chances of encountering an outlying value of the auxiliary statistic are small, because only rare random draws generate outliers, by definition, and a fairly small number of draws are used. However, when a large number of auxiliary statistics are generated, as is the case with the SBIL estimator, outliers will eventually appear if the distribution of the auxiliary statistic has outliers in its support.
INDIRECT LIKELIHOOD INFERENCE 32 achieved by working with a finite dimensional statistic rather than with the full sample converts potentially infinite dimensional problems (as the sample grows) into tractable finite dimensional problems. This is an important simplification when nonparametric estimation methods are used. Moreover, with a careful choice of auxiliary statistic, once can hope for approximate sufficiency. As we have seen in the DSGE examples, the SBIL estimator may be computed even when the auxiliary statistic is of fairly high dimension, at the cost of requiring more simulations. The possibility of using a fairly high (but finite) dimensional auxiliary statistic makes it reasonably hopeful that the statistic approximately spans the space of the efficient score, in which case the SBIL estimator will be approximately fully asymptotically efficient. Our Monte Carlo results for the dynamic panel and nonlinear panel examples can be compared to the results of other authors for other estimators, giving support to the good relative efficiency of the SBIL estimator. Additional support comes from our MA example, where the SBIL estimator often exhibits an RMSE smaller than that of the ML estimator. The fact that the SBIL estimator may have better small sample performance than the ML estimator may be relevant when one seeks to estimate complex DSGE models. The combination of particle filtering and MCMC discussed above seeks to compute the ML estimator or related Bayesian likelihood-based estimators. The filtering/MCMC technology is relatively complicated to implement, and is computationally extremely demanding. In comparison, the SBIL estimator is simple to implement. In addition, it is certainly possible that the SBIL estimator could have better small sample performance than the ML estimator of such complex and often nonlinear models. An interesting avenue to explore would be to compare our estimator with the the MLE based on particle filtering/MCMC alternatives In our implementation, we have focused on the basic sampler as given in equation (6) choosing the number of neighbors kthrough the simple rule k=1.5 ×S0.25. There is certainly scope for use of more sophisticated rules, such as cross-validation, or different kernels, which could lead to better performance. Similarly, more complicated samplers using importance sampling methods could be used to improve on the computation time. We leave these numerical issues for future research.
INDIRECT LIKELIHOOD INFERENCE 33 REFERENCES [1] Altonji, J. and L.M. Segal, 1996, “Small sample bias in GMM estimation of covariance structures,” Journal of Economic and Business Statistics 14, 353-366. [2] An, S. and F. Schorfheide, 2007, “Bayesian Analysis of DSGE Models”, Econometric Reviews, 26, 113-172. [3] Andrews, D.W.K., 1993, “Exactly median-unbiased estimation of first order autoregressive/unit root models,” Econometrica 61, 139–165. [4] Arellano, M. and S. Bonhomme, 2009, “Robust priors in nonlinear panel data models,” Econometrica, 77, 489-536. [5] Aruoba, S.B, J. Fernández-Villaverde and J. Rubio-Ramírez, 2006, “Comparing solution methods for dynamic equilibrium economies”, Journal of Economic Dynamics & Control, 30, 2477-2508. [6] Arya, S., T. Malamatos, and D.M. Mount, 2009, “Space-time tradeoffs for approximate nearest neighbor searching, Journal of the ACM 57, 1-54. [7] Bahadur, R., S. Zabell, and J. Gupta, 1980, “Large deviations, tests, and estimates,” in: I.M. Chaterabarli (Ed.), Asymptotic Theory of Statistical Tests and Estimation, pp. 33–64,. New York: Academic Press. [8] Beaumont, M., W. Zhang and D. Balding, 2002, “Approximate Bayesian computation in population genetics”, Genetics, 162, 2025-2035. [9] Beaumont, M., J.-M. Cornuet, J.-M. Marin and C. Robert, 2009, “Adaptive approximate Bayesian computation”, Biometrika, 96, 983-990. [10] Bhattacharya, R.N. and J.K. Ghosh, 1978, “On the validity of the formal Edgeworth Expansion,” Annals of Statistics 6, 434-451. [11] Bhattacharya, R.N. and R. R. Rao, 1976, Normal Approximations and Asymptotic Expansions. New York: Wiley. [12] Bickel, P.J., F. Götze, and W.R. van Zwet, 1985, “A simple analysis of third-order efficiency of estimates,” in L. Le Cam and R.A. Olshen (Eds.), Proceedings of the Berkeley Conference in Honor of Jerzy Neyman and Jack Kiefer. Wadsworth. [13] Canova, F. and L. Sala, 2009, “Back to square one: Identification issues in DSGE models,” Journal of Monetary Economics, 56, 231-249. [14] Chernozhukov, V. and H. Hong, 2003, “An MCMC approach to classical estimation,” Journal of Econometrics 115, 293-346. [15] Chumacero, R., 2001, “Estimating ARMA models efficiently”, Studies in Nonlinear Dynamics and Econometrics, 5, 103-114. [16] Collomb, G. and W. Härdle, 1986, “Strong uniform convergence rates in robust nonparametric time series analysis and prediction: kernel regression estimation from dependent observations,” Stochastic Processes and Their Applications 23, 77-89. [17] Creel, M., 2007, “I ran four million probits last night: HPC clustering with ParallelKnoppix,” Journal of Applied Econometrics, 22, 215-223. [18] Creel, M. and D. Kristensen, 2009, “Estimation of dynamic latent variable models using simulated nonparametric moments,” UFAE and IAE Working Paper 792.09, http://ideas.repec.org/p/aub/autbar/792.09.html. [19] Donald, S. and W.K. Newey, 2000, “A Jackknife Interpretation of the Continuous Updating Estimator,” Economics Letters 67, 239-243. [20] Doran, H.E. and P. Schmidt, 2006, “GMM estimators with improved finite sample properties using principal components of the weighting matrix, with an application to the dynamic panel data model,” Journal of Econometrics, 133, 387–409. [21] Duffie, D. and K. J. Singleton, 1993, “Simulated moments estimation of Markov models of asset prices,” Econometrica, 61, 929–952. [22] Everaert, G. and L. Pozzi, 2007, “Bootstrap-based bias correction for dynamic panels”, Journal of Economic and Dynamics Control 31, 1160-1184.
INDIRECT LIKELIHOOD INFERENCE 34 [23] Fermanian, J.-D. and B. Salanié, 2004, “A nonparametric simulated maximum likelihood estimation method,” Econometric Theory, 20, 701-734. [24] Fuh, C.-D., 2006, “Efficient likelihood estimation in state space models,” Annals of Statistics 34, 2026-2068. [25] Gallant, A. R. and G. Tauchen, 1996, “Which moments to match?” Econometric Theory 12, 657681. [26] Ghosh, J.K. (1994) Higher Order Asymptotics, Hayward: IMS. [27] Gouriéroux, C., A. Monfort, and E. Renault, 1993, “Indirect inference,” Journal of Applied Econometrics, 8, S85-S118. [28] Gouriéroux, C., P.C.B. Phillips and J. Yu, 2010, “Indirect inference for dynamic panel models,” Journal of Econometrics 157, 68-77. [29] Gouriéroux, C., E. Renault and N. Touzi, 2000, “Calibration by simulation for small sample bias correction,” in Mariano, R.S., Schuermann, T., Weeks, M. (Eds.), Simulation-Based Inference in Econometrics: Methods and Applications, pp. 328-358. Cambridge: Cambridge University Press. [30] Guerron, P., 2010, “What you match does matter: The effects of Data on DSGE Estimation,” Journal of Applied Econometrics 25, 774-804. [31] Gusev, S.I., 1975, “Asymptotic expansions associated with some statistical estimators in the smooth case I: Expansions of random variables,” Theory of Probability and Its Applications 20, 470-498. [32] Hall, P., 1992, The Bootstrap and Edgeworth Expansion, New York: Springer. [33] Hall, P. and J.L. Horowitz, 1996, “Bootstrap critical values for tests based on GeneralizedMethod-of-Moments estimators,” Econometrica 64, 891-916. [34] Hahn, J. and G. Kuersteiner, 2002, “Asymptotically unbiased inference for a dynamic model with fixed effects when both nand Tare large,” Econometrica 70, 1639-1657. [35] Hahn, J. and W.K. Newey, 2004, “Jackknife and analytical bias reduction for nonlinear panel models”, Econometrica 72, 1295-1319. [36] Hansen, L.P., J. Heaton and A. Yaron, 1996, “Finite-sample properties of some alternative GMM estimators”, Journal of Business and Economic Statistics 14, 262-280. [37] Horowitz, J. L., 1992, “A smoothed maximum score estimator for the binary response model,” Econometrica 60, 505-531. [38] Inoue, A. and M. Shintani, 2006, “Bootstrapping GMM estimators for time series,” Journal of Econometrics 133, 531–555. [39] Karagedikli, Ö., T. Matheson, C. Smith, C. and S.P. Vahey, 2010, “RBCs AND DSGEs: the computational approach to business cycle theory and evidence,” Journal of Economic Surveys, 24, 113–136. [40] Kezdi, G., J. Hahn and G. Solon, 2002, “Jackknife minimum distance estimation,” Economics Letters 76, 35-45. [41] Kormiltsina, A. and D. Nekipelov, 2009, “Numerical performance of MCMC algorithms for classical estimation,” working paper, UC Berkeley. [42] Kristensen, D., 2009, “Uniform convergence rates of kernel estimators with heterogeneous, dependent data,” Econometric Theory 25, 1433-1445. [43] Kristensen, D. and B. Salanié, 2010, “Higher order improvements for approximate estimators,” CAM Working Papers 2010-04, University of Copenhagen. [44] Kristensen, D. and Y. Shin, 2008, “Estimation of dynamic models with nonparametric simulated maximum likelihood,” CREATES Research Papers 2008-58, University of Aarhus. [45] Li, Q. and J. Racine, 2007, Nonparametric Econometrics: Theory and Practice. Princeton: Princeton University Press. [46] Li, T., 2010, “Indirect inference in structural econometric models,” Journal of Econometrics 157, 120-128.
INDIRECT LIKELIHOOD INFERENCE 35 [47] Mancini, T., 2010, “Dynare user guide: An introduction to the solution & estimation of DSGE models”, http://www.dynare.org/documentation-and-support/user-guide/. [48] Marjoram, P., J. Molitor, V. Plagnol and S. Tavaré, 2003, “Markov chain Monte Carlo without likelihoods”, Proceedings of the National Academy of Sciences, USA, 100, 15324-15328. [49] McFadden, D., 1989, “A method of simulated moments for estimation of discrete response models without numerical integration,” Econometrica, 57, 995–1026. [50] Newey, W.K. and D. McFadden, 1994, “Large sample estimation and hypothesis testing,” in: R. Engle and D. McFadden (Eds.), Handbook of Econometrics, Vol. IV, 2111-2245. Amsterdam: Elsevier Science. [51] Newey, W. K. and R. J. Smith, 2004, “Higher order properties of GMM and generalized empirical likelihood estimators,” Econometrica, 72, 219–255. [52] Pfanzagl, J. and W. Wefelmeyer, 1978, “A third-order optimum property of the maximum likelihood estimator,” Journal of Multivariate Analysis, 8, 1-29. [53] Phillips, P.C.B, 1977, “A General Theorem in the Theory of Asymptotic Expansions as Approximations to Finite Sample Distributions of Econometric Estimators,” Econometrica 45, 15171534. [54] Rilstone, P., V.K. Srivastava and A. Ullah, 1996, “The second-order bias and mean squared error of nonlinear estimators,” Journal of Econometrics 75, 369-395. [55] Rothenberg, T.J. (1984), “Approximating the distributions of econometric estimators and test statistics,” in: Z. Griliches and M.D. Intriligator (Eds.), Handbook of Econometrics, Vol. II, 881935. Amsterdam: Elsevier Science. [56] Ruge-Murcia, F., 2007, “Methods to estimate dynamic stochastic general equilibrium models, Journal of Economic Dynamics and Control, 31, 2599-2636. [57] Ruge-Murcia, F., 2010, “Estimating nonlinear DSGE models by the simulated method of moments”, working paper, Cahier 19-2010, CIREQ. [58] Sisson, S., Y. Fan and M. Tanaka, 2007, “Sequential Monte Carlo without likelihoods”, Proceedings of the National Academy of Science, USA, 104, 1760-1765. [59] Skovgaard, I., 1981, “Transformation of an Edgeworth expansion by a sequence of smooth functions,” Scandinavian Journal of Statistics 8, 207-217. [60] Skovgaard, I., 1986, “On multivariate Edgeworth expansions,” International Statistical Review 54, 169-186. [61] Smith, A., 1993, “Estimating nonlinear time series models using simulated vector autoregressions,” Journal of Applied Econometrics, 8, S63-S84. [62] Tavaré, S., D. Balding, R. Griffiths and P. Donnelly, 1997, “Inferring coalescence times from DNA sequence data”, Genetics, 145, 505-518. [63] Winschel, V. and Krätzig, M., 2010, “Solving, estimating, and selectingnonlinear dynamic models without the curse of dimensionality,” Econometrica, 78, 803–821. [64] Zeitouni, O. and M. Gutman, 1991, “On universal hypothesis testing via large deviations,” IEE Transactions on Information Theory 37, 285–290. [65] Zhang, P., 1996, “Nonparametric Importance Sampling,” Journal of the American Statistical Association 91, 1245-1253.
INDIRECT LIKELIHOOD INFERENCE 36 APPENDIX A: PROOFS Proof. [Proposition 1]We first investigate the MIL estimator: To this end, first note that by Lemma 1the log-likelihood satisfies 1 nlog f(Zn|θ) = 1 nlogφ∗ n(Zn|θ)+LRn(θ)=1 nlogφ∗ n(Zn|θ)+oP1/√n uniformly in θ. Thus, for the first-order analysis, we can treat Ln(θ):=logφ∗ n(Zn|θ)as the actual log-likelihood. To show consistency, note that uniformly in θ∈Θ: 1 nLn(θ) = −1 2nlog (|Ω(θ)|)−Tn(θ)0Tn(θ) 2n+oP(1) =−1 2(Z(θ0)−Z(θ))0Ω−1(θ) (Z(θ0)−Z(θ)) +oP(1) (21) =:L(θ)+oP(1), where L(θ)is a continuous function with a unique minimum at θ=θ0by Assumption 3. It now follows by standard results (see e.g. Newey and McFadden, 1994, Theorem 2.1), that the MLE is consistent. Next, we show asymptotic normality: With ˙ Z(i)(θ)=∂Z(θ)/(∂θi)and ˙ Ω(i)(θ)= ∂Ω(θ)/(∂θi), ∆n,i(θ):=∂Ln(θ) ∂θi =−1 2Ω−1(θ)˙ Ω(i)(θ)−√nTn(θ)0Ω−1/2 (θ)˙ Z(i)(θ)+1 2Tn(θ)0˙ Ω(i)(θ)Tn(θ) =−√nTn(θ)0Ω−1/2 (θ)˙ Z(i)(θ)+oP√n and with ¨ Z(i,j)(θ)=∂2Z(θ)/∂θi∂θjand ¨ Ω(i,j)(θ)=∂2Ω(θ)/∂θi∂θj, Jn,ij (θ):=1 n ∂2Ln(θ) ∂θi∂θj =1 nΩ−2(θ)˙ Ω(i)(θ)˙ Ω(j)(θ)−1 nΩ−1(θ)¨ Ω(i,j)(θ) +˙ Z(i)(θ)0Ω−1(θ)˙ Z(j)(θ)+Tn(θ)0Ω−1/2 (θ)¨ Z(i,j)(θ)/√n+oP(1), With J(θ)defined in Assumption 3, it now holds that (22) 1 √n∆n(θ0)=−Tn(θ0)0Ω−1/2 (θ0)˙ Z(θ0)+oP(1)→dN(0, J(θ0)) , and, uniformly in θ,Jn(θ)=J(θ)+oP(1). Since the score of the log-likelihood converges weakly towards a normal distribution while the Hessian converges uniformly towards a non-singular limit in probability, it now follows by a standard Taylor expansion of the score that the MILE is √n-asymptotically normally distributed with asymptotic variance J−1(θ0). Next, the properties of the BIL are established by verifying Assumptions 1-4 in Chernozhukov and Hong (2003), CH henceforth, with Ln(θ)chosen as above. First note that CH’s Assumptions 1-2 are satisfied by our Assumption 1. What remains is to verify their Assumption 3-4. But by combining their Lemmas 1-2 with the above derivations, these
INDIRECT LIKELIHOOD INFERENCE 37 are easily verified. We can now appeal to CH’s Theorem 2 which yields the desired result. Proof. [Proposition 2]This follows directly from Chernozhukov and Hong (2003, Theorem 3) since eqs. (22) and Jn(θ)=J(θ)+oP(1)imply that the generalized information equality holds. Proof. [Proposition 3]By assumption, ¯ Zn(θ)=Z(θ)+o1/√n, while Zn=Zn(θ0)→P Z(θ0). Thus, Dn(θ) = D(θ) + op(1), where the limit is given by D(θ) = 1 2(Z(θ0)−Z(θ))0Ω−1(θ0) (Z(θ0)−Z(θ)) . By Assumption 3in conjunction with standard arguments, it now follows that ˆ θGMM is consistent. To derive its asymptotic distribution, first note that ˆ θGMM solves 0=∂Dn(θ) ∂θ0=−∂Zn(θ) ∂θ 0Wn(Zn−¯ Zn(θ)) =−˙ Z(θ)0Wn(Zn−Z(θ)) +oP1/√n, where, by Assumption 2, Zn−Z(θ)=Zn−Z(θ0)−˙ Zθ(θ−θ0), where θlies on the line between θand θ0. Combining these two equations, 0=−˙ Zˆ θGMM0WnZn−Zˆ θGMM+oP1/√n =−˙ Zˆ θGMM0Wn{Zn−Z(θ0)}+˙ Zˆ θGMM0Wn˙ Zθ(ˆ θGMM −θ0) + oP1/√n. The result now follows by Assumption 2together with Wn→PΩ−1(θ0). Proof. [Proposition 4]First, consider the two-step GMM estimator, ˆ θGMM. With mn(θ)=(Zn−¯ Zn(θ))0Wn ∂¯ Zn(θ) ∂θ , we can apply Lemma 2. The first and second order derivatives are given by ∂mn(θ0) ∂θ =−∂¯ Zn(θ0)0 ∂θ Wn ∂¯ Zn(θ0) ∂θ +(Zn−¯ Zn(θ))0Wn ∂2¯ Zn(θ0) ∂θ2, and ∂2mn(θ0) ∂θ2=−3∂2¯ Zn(θ0)0 ∂θ2Wn ∂¯ Zn(θ0) ∂θ +(Zn−¯ Zn(θ))0Wn ∂3¯ Zn(θ0) ∂θ3 With D¯ mn=−∂¯ Zn(θ0)0 ∂θ Ω−1 n(θ0)∂¯ Zn(θ0) ∂θ , D2¯ mn=−3∂2¯ Zn(θ0) ∂θ2Ω−1 n(θ0)∂¯ Zn(θ0) ∂θ ,
INDIRECT LIKELIHOOD INFERENCE 38 where Ω−1 n(θ0)denotes the variance of Zn−¯ Zn(θ), and ∆ndefined in the proposition, An:=∂mn(θ0) ∂θ −D¯ mn =∂2¯ Zn(θ)0 ∂θ2Ω−1 n(θ0) (Zn−¯ Zn(θ)) −∂¯ Zn(θ0)0 ∂θ ∆n ∂¯ Zn(θ0) ∂θ +∂2¯ Zn(θ)0 ∂θ2∆n(Zn−¯ Zn(θ)) =∂2¯ Zn(θ)0 ∂θ2Ω−1 n(θ0) (Zn−¯ Zn(θ)) −∂¯ Zn(θ0)0 ∂θ ∆n ∂¯ Zn(θ0) ∂θ +OP(1/n). Thus, E[Anmn(θ0)] =∂2¯ Zn(θ)0 ∂θ2Ω−1 n(θ0)Eh(Zn−¯ Zn(θ)) (Zn−¯ Zn(θ))0iΩ−1 n(θ0)∂¯ Zn(θ) ∂θ −∂¯ Zn(θ0)0 ∂θ E∆n ∂¯ Zn(θ0) ∂θ (Zn−¯ Zn(θ))0Ω−1 n(θ0)∂¯ Zn(θ) ∂θ =1 n ∂2¯ Zn(θ)0 ∂θ2Ω−1 n(θ0)∂¯ Zn(θ) ∂θ −1 n ∂¯ Zn(θ0)0 ∂θ E∆n ∂¯ Zn(θ0) ∂θ (Zn−¯ Zn(θ))0Ω−1 n(θ0)∂¯ Zn(θ) ∂θ ≃1 n1 3D2¯ m+BW,n with BW,ndefined in the proposition. The other bias component can be written as: Em2 n(θ0)=∂¯ Zn(θ)0 ∂θ Ω−1 n(θ0)Eh(Zn−¯ Zn(θ)) (Zn−¯ Zn(θ))0iΩ−1 n(θ0)∂¯ Zn(θ) ∂θ =1 n ∂¯ Zn(θ)0 ∂θ Ω−1 n(θ0)∂¯ Zn(θ) ∂θ ≃ −1 nJ(θ0). Thus, by Lemma 2, Eˆ θGMM−θ0≃ −J−2(θ0)E[Anmn(θ0)] −1 2 D2¯ mn D¯ mn Em2 n(θ0). ≃1 nJ−2(θ0)1 6D2¯ m+BW,n Next, consider the CU estimator: It is easily checked that the expansion goes through with mn(θ):=2∂¯ Zn(θ)0 ∂θ Ω−1 n(θ) (Zn−¯ Zn(θ)) +(Zn−¯ Zn(θ))0∂Ω−1 n(θ) ∂θ (Zn−¯ Zn(θ)) , and D¯ mnand D2¯ mngiven as before. However, in the case of CU, An:=∂mn(θ0) ∂θ −D¯ mn=∂2¯ Zn(θ)0 ∂θ2Ω−1 n(θ0) (Zn−¯ Zn(θ)) +OP(1/n) and so the bias term due to the first-step estimation of the weighting matrix vanishes and we obtain the claimed result.
INDIRECT LIKELIHOOD INFERENCE 39 Finally, consider the MIL estimator: Since LR (θ)=oP1/n2, we can choose mn(θ)= n−1∂logf∗ n(Zn|θ)/(∂θ)such that ∂mn(θ) ∂θ =1 n ∂2logf∗ n(Zn|θ) ∂θ2,∂2mn(θ) ∂θ2=1 n ∂3logf∗ n(Zn|θ) ∂θ3. From the definition of f∗ n(Zn|θ),mn(θ)=mn,1 (θ)+mn,2 (θ), where the first term is the Gaussian component, mn,1 (θ)≃˙ Z(θ)0Ω−1(θ) (Zn−Z(θ)) , while the second one is due to the higher-order component, mn,2 (θ)≃1 n3/2 ∂π1(Tn(θ)|θ)/∂θ 1+π1(Tn(θ)|θ)/√n ≃1 n3/2 ∂π1(Tn(θ)|θ) ∂θ ≃ −1 nπ(1) 1(Tn(θ)|θ)Ω−1/2 (θ)˙ Z(θ). The derivatives satisfy ∂m1,n(θ) ∂θ ≃1 2¨ Z(θ)0Ω−1(θ) (Zn−Z(θ)) −˙ Z(θ)0Ω−1(θ)˙ Z(θ), ∂m2,n(θ) ∂θ ≃ −1 nπ(1) 1(Tn(θ)|θ)Ω−1/2 (θ)¨ Z(θ) +1 √n˙ Z(θ)0Ω−1/2 (θ)π(2) 1(Tn(θ)|θ)Ω−1/2 (θ)˙ Z(θ), and ∂2m1,n(θ) ∂θ2≃... Z(θ)0Ω−1(θ) (Zn−Z(θ)) −3¨ Z(θ)0Ω−1(θ)˙ Z(θ), ∂2m2,n(θ) ∂θ2≃ −1 nπ(1) 1(Tn(θ)|θ)Ω−1/2 (θ)... Z(θ) +2 √n˙ Z(θ)0Ω−1/2 (θ)π(2) 1(Tn(θ)|θ)Ω−1/2 (θ)˙ Z(θ) +∑ i ˙ Z(θ)0Ω−1/2 (θ)˜ π(3) 1,i(Tn(θ)|θ)Ω−1/2 (θ)˙ Z(θ), where ˜ π(3) 1,i(Tn(θ)|θ)=∂π(2) 1(Tn(θ)|θ)/(∂ti)˙ Tn,i(θ)/√n. Since ¯ Tn(θ0)=O1/√n, we can choose D¯ mnand D¯ m2 nas for the GMM and CU estimators except that Z(θ)replaces ¯ Zn(θ). Next, in order to obtain an expression of the bias, we Taylor-expanding w.r.t. the statistic: With ¯ Z0,n:=¯ Zn(θ0),¯ f∗ n:=¯ f∗ n(¯ Z0,n|θ0)and ¯ Tn(θ):=√nΩ−1/2 (¯ Z0,n−Z(θ)), ∂imn(θ) ∂θi≃1 n ∂i+1log ¯ f∗ n(θ) ∂θi+1 n ∂i+2log ¯ f∗ n(θ) ∂θi∂z(Zn−¯ Z0,n), for i=0,1,2, where 1 n ∂2log ¯ f∗ n(θ) ∂θ∂z≃˙ Z(θ)0Ω−1(θ)+1 √nΩ−1/2 (θ)π(2) 1(¯ Tn(θ)|θ)Ω−1/2 (θ)˙ Z(θ),
INDIRECT LIKELIHOOD INFERENCE 40 1 n ∂3log ¯ f∗ n(θ) ∂θ2∂z≃¨ Z(θ)0Ω−1(θ)−1 √nΩ−1/2 (θ)π(2) 1(¯ Tn(θ)|θ)Ω−1/2 (θ)¨ Z(θ) +∑ i ˙ Z(θ)0Ω−1/2 (θ)¯ π(3) 1,i(¯ Tn(θ)|θ)Ω−1/2 (θ)˙ Z(θ), where ¯ π(3) 1,i(t)is defined in the proposition, and 1 n ∂4log ¯ f∗ n(θ) ∂θ3∂z≃... Z(θ)0Ω−1(θ)+1 √nΩ−1/2 (θ)π(2) 3(¯ Tn(θ)|θ)Ω−1/2 (θ)... Z(θ) +2∑ i ˙ Z(θ)0Ω−1/2 (θ)¯ π(3) 3,i(¯ Tn(θ)|θ)Ω−1/2 (θ)˙ Z(θ). We note that ¯ Tn(θ0)=O1/√nsuch that π(i) 1(¯ Tn(θ)) ≃π(i) 1(0). Thus, E[Anmn(θ0)] ≃1 n2 ∂2log ¯ f∗ n(θ0) ∂θ∂z0Eh(Zn−¯ Zn(θ0)) (Zn−¯ Zn(θ0))0i∂3log ¯ f∗ n(θ0) ∂θ2∂z ≃1 3nD2¯ m+1 nBπ, and Em2 n(θ0)≃1 n2 ∂2log ¯ f∗ n(θ0) ∂θ∂z0Eh(Zn−¯ Zn(θ0)) (Zn−¯ Zn(θ0))0i∂2log ¯ f∗ n(θ0) ∂θ∂z≃1 nJ(θ0). Lemma 2now yields the claimed result. Proof. [Proposition 5]As usual, we can treat f∗(Zn|θ)as the actual likelihood due to Lemma 1. By an rth order Taylor expansion of the corresponding score equation w.r.t. θ, (23) 0 =Wn,1 (Zn)+ r ∑ i=1 1 i!Wn,i(Zn)ˆ θMIL −θ0i+Rn,=:AWn(Zn),ˆ θMIL+Rn, where Wn(z)=(Wn,1 (z), ...,Wn,r(z)) with Wn,i(z)=n−1∂ilogf∗ n(z|θ0)/∂θi 0, and Rn=n−1|∂rlogf∗ n(z|θ)/(∂θr)|θ=θˆ θMIL −θ0r. First, ignore Rnand redefine ˆ θMIL as the solution to AWn(Zn),ˆ θMIL=0. From the expression of f∗(Zn|θ), it is easily seen that Wn(Z(θ0)) =W∞(Z(θ0)) + r ∑ i=1 1 ni/2 Mi+on−r/2, where Miare constants depending on derivatives of the polynomials π1, ..., πrand W∞,i(Z(θ0)) is the leading term of n−1∂ilog f∗(Z(θ0)|θ0)/∂θi 0. In particular, the limiting score and Hessian satisfy ¯ W∞,1 (Z(θ0)) =0 and ¯ W∞,2 (Z(θ0)) =−J(θ0). Thus, A(W∞(Z(θ0)) ,θ0)= 0, and ∂A(W∞(Z(θ0)) ,θ)/∂θ|θ=θ0=−J(θ0)has full rank. Hence, by the implicit function theorem, there exists an analytic function H(w)in a neighborhood of W∞(Z(θ0)) such that θ0=H(W∞(Z(θ0))). Moreover, for all nlarge enough, the solution θ0,nto A(Wn(Z(θ0)) ,θ0,n)=0, can be expressed as θ0,n=H(Wn(Z(θ0))) since Wn(Z(θ0))
INDIRECT LIKELIHOOD INFERENCE 41 lies in a neighborhood of W∞(Z(θ0)) for all nlarge enough. The sequence θ0,nsatisfies θ0,n−θ0=H(Wn(Z(θ0))) −H(W∞(Z(θ0))) = r ∑ i=1 ∂iH(W∞(Z(θ0))) ∂wi[Wn(Z(θ0)) −W∞(Z(θ0))]i+on−r/2 =:r ∑ j=1 1 nj/2 ˜ Mj+on−r/2, where ˜ Mjis a constant depending on M1,..., Mrand the first rderivatives of H(W∞(Z(θ0))), j=1,...,r. We obtain an Edgeworth expansion ofˆ θMIL −θ0,n=H(Wn(Zn)) −H(Wn(Z(θ0))) by applying the general result of Phillips (1977) for Edgeworth expansions of transformations of random sequences: We define the following sequence of functions en(q):=H(Wn(q+Z(θ0))) −H(Wn(Z(θ0))) , such that en:=ˆ θMIL −θ0,n=en(qn), where qn:=Zn−Z(θ0), and verify Phillips (1977, Assumptions 3-5): First, since the distribution of the normalized statistic Tn(θ0)=√nqn satisfies an Edgeworth expansion by Assumption 4, Phillips (1977, Assumption 3) holds. Next, the two function Hand Wnare both rtimes continuously differentiable and the derivatives of Wn(z)converges towards those of W∞(z). Thus, en(q)is rtimes differentiable with its derivatives uniformly bounded in a neighborhood around 0. Finally, we know from the implicit function theorem that ∂H(W∞(Z(θ0))) /(∂w)has full rank while it is easily checked that ∂W∞,1 (Z(θ0)) /(∂z)=Ω−1/2 (θ0)˙ Z(θ0). Hence, by the chain rule, |∂en(q)/∂q|is bounded away from zero as n→∞. This shows that Phillips (1977, Assumptions 4-5) hold. We have shown that √nenadmits an Edgeworth expansion, say f∗ en(x)=φ(x)"1+ r ∑ i=1 n−i/2 ¯ πi(x)#. This in turn implies that the distribution of en:=√nˆ θMIL −θ0=√nen+bn, where bn=√n(θ0,n−θ0)=∑r j=1n−j/2Mj+on−r/2, can be approximated by f∗ en(x)=φ(x−bn)"1+ r ∑ i=1 n−i/2 ¯ πi(x−bn)#. Expanding around f∗ en(x)and rearranging terms, we then obtain the desired result where the coefficients of the polynomial ˜ πi(x)depend on the ones of ¯ πj(x)and the coefficients Mj,j=1,...,r. Finally, we have to verify that we are allowed to ignore the remainder term Rnin the Taylor expansion. By the arguments in Rothenberg (1984, p. 898), this will follow if P(|Rn|>logcn)=on−r/2. This will in turn hold if P|Wn(Zn)−Wn(Z(θ0))|>c1qlog (n)/n=on−r/2,
INDIRECT LIKELIHOOD INFERENCE 48 TABLE 6. Fully observed DSGE model with monopolistic competition Bias RMSE Parameter Lower Bound Upper bound True values Prior mean SBIL Prior mean SBIL α0.15 0.4 0.33 -0.055 -0.002 0.091 0.006 β0.95 0.999 0.99 -0.016 -0.000 0.021 0.001 δ0.005 0.06 0.023 0.010 0.001 0.019 0.001 ψ1 3 1.75 0.250 0.004 0.629 0.014 ρ0.85 0.99 0.95 -0.030 -0.017 0.050 0.024 σ0.005 0.04 0.01 0.012 0.000 0.016 0.001 e9 13 10 1.000 0.002 1.529 0.033 TABLE 7. Partially observed DSGE model with habit formation, first design Bias RMSE Parameter Lower Bound Upper bound True values Prior mean SBIL Prior mean SBIL α0.25 0.4 0.36 -0.035 -0.003 0.056 0.008 β0.93 0.99 0.95 0.010 0.001 0.020 0.004 δ0.02 0.04 0.025 0.005 0.000 0.008 0.001 η0 0.5 0.2 0.050 -0.024 0.153 0.044 γ1 4 2 0.500 0.220 1.000 0.283 ρ0.8 0.99 0.85 0.045 -0.005 0.071 0.017 σ0.01 0.08 0.04 0.005 0.001 0.021 0.004 ψNA NA 3.197 9.854 0.356 22.529 0.530 TABLE 8. Partially observed DSGE model with habit formation, second design Bias RMSE Parameter Lower Bound Upper bound True values Prior mean SBIL Prior mean SBIL α0.25 0.4 0.36 -0.035 0.001 0.056 0.006 β0.93 0.99 0.95 0.010 -0.003 0.020 0.004 δ0.02 0.04 0.025 0.005 0.001 0.008 0.002 η0 0.5 0.4 -0.150 -0.033 0.208 0.056 γ1 4 3 -0.500 0.018 1.000 0.185 ρ0.8 0.99 0.85 0.045 -0.003 0.071 0.016 σ0.01 0.08 0.04 0.005 0.000 0.021 0.003 ψNA NA 13.562 -0.511 0.792 22.266 3.052
INDIRECT LIKELIHOOD INFERENCE 49 FIGURE 1. Fully observed DSGE model. Pseudo-priors, true parameter values, and density of SBIL (A)α(B)β(C)δ (D)ψ(E)ρ(F)σ (G)e FIGURES UNIVERSITAT AUTÒNOMA DE BARCELONA AND MOVE COLUMBIA UNIVERSITY AND CREATES (CENTER FOR RESEARCH IN ECONOMETRIC ANALYSIS OF TIME SERIES, UNIVERSITY OF AARHUS).