Importance sampling squared for Bayesian inference in latent variable models

Introduction

A wide range of statistical models, such as nonlinear and non-Gaussian state space models and generalized linear mixed models, lead to analytically intractable likelihood functions. When the density of the observations conditional on the parameters and a vector of latent variables is available in closed form, we can use importance sampling (IS) and, more generally, sequential importance sampling to estimate the likelihood unbiasedly. Our article considers importance sampling for Bayesian inference when working with an estimated likelihood. We call this procedure importance sampling squared (IS2).

A key motivation for the IS2 method is that it estimates the marginal likelihood (and the standard error of this estimator) accurately and efficiently in cases where only an estimate of the likelihood is available. The marginal likelihood is a fundamental tool for Bayesian model comparison (Kass and Raftery,, 1995), but it is challenging to use the output of Markov Chain Monte Carlo (MCMC) methods to estimate it, especially for models with an intractable likelihoods (Chib and Jeliazkov,, 2001; Perrakis et al.,, 2014). The IS2 method obtains the marginal likelihood automatically when estimating expectations with respect to the posterior. Moreover, the method can be used as a tool for marginal likelihood estimation even if we estimate the posterior distribution by MCMC methods. In this case, it is typically straightforward to use the MCMC output to form a proposal density for the parameters to use in the IS2 method to estimate the marginal likelihood.

Our article shows that IS is still valid for estimating expectations with respect to the posterior when the likelihood is estimated unbiasedly, and prove a law of large numbers and a central limit theorem for these estimators. This analysis relates directly to the results in Fearnhead et al., (2008) and Fearnhead et al., (2010), who considered random weight importance sampling in the context of particle filtering. Our results allow us to analyze how much asymptotic efficiency is lost when working with an estimated likelihood by comparing the asymptotic variance obtained under IS2 to the standard case where the likelihood is known. We show that the ratio of the asymptotic variance of the IS2 estimator to the asymptotic variance of the IS estimator, which we call the inflation factor, is greater than or equal to 1, and is equal to 1 if and only if the likelihood is known. The inflation factor increases exponentially with the variance of the estimator of the log-likelihood.

A critical implementation issue is the choice of the number of particles $N$ for estimating the likelihood. A large $N$ gives a more accurate estimate of the likelihood at greater computational cost, while a small $N$ can lead to an estimator with a very large variance. We provide theoretical and practical guidelines on how to select $N$ to obtain an optimal tradeoff between accuracy and computational cost. Our results show that the efficiency of $\text{\rm IS}^{2}$ is weakly sensitive to the number of particles around the optimal value. Moreover, the loss of efficiency decreases at worst linearly when we choose $N$ higher than the optimal value, whereas the efficiency can deteriorate exponentially when $N$ deviates appreciably below its optimal value. We therefore advocate a conservative choice of $N$ in practice. We propose two approaches for selecting the number of particles. The first approach is static because it selects the same number of particles for all parameter values. The second approach is dynamic because it selects the optimal number of particles depending on the parameter values. We show both theoretically, in section S3 of the supplementary material, and empirically, in section S8 of the supplementary material, that the dynamic approach can be much more efficient than the static approach.

Our method relates to alternative approaches to Bayesian inference for models with intractable likelihoods. Beaumont, (2003) develops a pseudo marginal Metropolis Hastings (PMMH) scheme to carry out Bayesian inference with an estimated likelihood. Andrieu and Roberts, (2009) formally study Beaumont’s method and give conditions under which the chain converges. Andrieu et al., (2010) use MCMC for inference in state space models where the likelihood is estimated by the particle filter, and Pitt et al., (2012) and Doucet et al., (2015) discuss the issue of the optimal number of particles to be used in likelihood estimation. The SMC2 method of Chopin et al., (2013) is a sequential Monte Carlo procedure for inference in space state models in which intractable likelihood contributions are estimated by random weight particle filters. The initialisation step in their algorithm corresponds to IS2 for an initial set of observations using the prior as a proposal for the parameters.

Given a statistical model with an intractable likelihood, we can choose to carry out off-line Bayesian inference using either the IS2 or MCMC as in Andrieu et al., (2010) and Pitt et al., (2012). There can be several advantages in following the IS2 approach. First, IS2 estimates the marginal likelihood. Second, it is straightforward to obtain the MC standard error of IS2 estimators since they are based on independent draws. In contrast, it can be difficult to obtain precise standard errors for estimators based on MCMC output due to autocorrelation in the Markov Chain. Third, it is straightforward to fully parallelize the IS2 procedure so that it can be more computationally attractive to use IS2 than PMMH when it is expensive to estimate the likelihood. Fourth, it is simple to implement variance reduction methods such as antithetic sampling, stratified mixture sampling and control variates for IS2, as well as to use randomized Quasi-Monte Carlo techniques to improve numerical efficiency; see the applications in section 5. Fifth, MCMC requires computationally expensive burn-in draws and assessing whether the Markov chain has converged and loses the information from rejected values (in standard implementations). IS2 uses all the draws from the proposal for the parameters. Sixth, high variance likelihood estimates can lead the PMMH Markov Chain to get trapped in certain regions of the extended sampling space, while isolated importance weights in IS2 can be directly diagnosed and stabilised using methods such as Pareto smoothed importance sampling (Vehtari and Gelman,, 2015). Finally, IS2 can deal with a multimodal posterior while PMMH methods can often get trapped in local modes.

Although the SMC2 method is designed for sequential updating, it can also be used effectively for off-line inference to deal with multimodal posteriors, for example. However, for off-line inference IS2 has advantages over SMC2. First, the sequential nature of the SMC2 algorithm can lead to a large implementation and computational effort compared to IS2. Second, IS2 provides estimates of the Monte Carlo standard errors for the posterior estimates, while it is not easy to do so with SMC2. Third, the IS2 method can take advantage of efficient off-line methods to estimate the likelihood, e.g. in state space models, which can outperform online likelihood estimation, as in SMC2, by orders of magnitude (e.g. Scharth and Kohn,, 2016). Third, because it is sequential, SMC2 can perform poorly in time series models with structural breaks or other interventions, unlike the PMMH and IS2 methods. However, we view IS2 and SMC2 as being complementary when used for prediction: we can use SMC2 for real time updating after starting from off-line IS2 estimates.

We illustrate the IS2 method in empirical applications for the generalized multinomial logit (GMNL) model of Fiebig et al., (2010) and a two factor stochastic volatility (SV) model with leverage effects. Our results for the GMNL model show that the IS2 approach is numerically more efficient than PMMH. In the standard implementation, the IS2 method leads 65-89% lower Monte Carlo mean-squared errors (MSEs) for estimating the posterior means compared to PMMH. When incorporating Quasi-Monte Carlo techniques, we obtain 87% to 99% reductions in MSE over the PMMH method. We show that the method for optimally selecting the number of particles in estimating the GMNL model leads to improved performance. The IS2 method accurately estimates the posterior distribution under our optimal implementation.

The SV application is based on daily returns of the S&P 500 index between 1990 and 2012. We show that the variance of the log-likelihood estimates based on the particle efficient importance sampling method of Scharth and Kohn, (2016) is small for this problem, despite the long time series. Hence, IS2 leads to highly accurate estimates of the posterior expectations for this example in a short amount of computing time. As suggested by the theory, the efficiency of the IS2 method is insensitive to the choice of $N$ in this example. As few as two particles for estimating the likelihood (including an antithetic draw) lead to an efficient procedure for this model, highlighting the practical convenience of the approach.

There is an online supplement to the article which contains results that complement those in the main paper. References to the main paper are of the form section 2, equation (2) and assumption 2, etc, whereas for the supplementary material we use section S2, equation (S2) and assumption S2, etc.

Latent variable models

This section sets out the main class of models with latent variables that we wish to estimate. However, the application of the method presented in this article is not limited to this class of models and section 5.2 applies the IS2 method to a time series stochastic volatility model. The method can also be applied in a variety of other contexts, including panel data estimation for large datasets using the subsampling approach of Quiroz et al., (2015).

Suppose there are $n$ individuals. Individual $i$ has $T_{i}$ observations $y_{i}=\{y_{i1},...,y_{iT_{i}}\}$, and for each $y_{it}$, there is a vector of covariates $x_{it}$. Let $y=\{y_{1},\dots,y_{n}\}$. We associate a latent vector $\alpha_{i}$ with individual $i$ and denote the vector of unknown parameters in the model by $\theta$. Assuming that the individuals are independent, the joint density of the latent variables and the observations is

By assigning a prior $p(\theta)$ for $\theta$, standard MCMC approaches can be used to sample from the joint posterior $p(\alpha,\theta|y)\propto p(y,\alpha|\theta)p(\theta)$ by cycling between $p(\alpha|\theta,y)$ and $p(\theta|\alpha,y)$. Typically, such models are estimated by introducing auxiliary latent variables as in Albert and Chib, (1993); Polson et al., (2013) for example. These methods require computationally expensive burn-in draws and the assessment of the convergence of the Markov chain. More importantly, they also suffer from the problem of slow mixing and are inefficient when the latent vector $\alpha$ is high dimensional or when the data is unbalanced (see Johndrow et al.,, 2016).

In this article, we are interested in inference when the likelihood is analytically intractable and given by

Section 5 considers more specific structures. We will use an unbiased IS estimator $\widehat{p}_{N}(y|\theta)$ of $p(y|\theta)$ based upon a simulation sample size of $N$, which we shall call the number of particles. Let $h_{i}(\alpha_{i}|y,\theta)$ be an importance density for $\alpha_{i}$. The density $p_{i}(y_{i}|\theta)$ is estimated unbiasedly by

Hence, $\widehat{p}_{N}(y|\theta)=\prod_{i=1}^{n}\widehat{p}_{N,i}(y_{i}|\theta)$ is an unbiased estimator of the likelihood $p(y|\theta)$. We note that it is possible to use a different number of particles $N_{i}$ for each individual $i$ and we shall do so in the empirical examples in section 5. We also consider the particle filter applied to state space models in section 5.2, which again yields an unbiased estimator of the likelihood $\widehat{p}_{N}(y|\theta)$.

When the likelihood cannot be evaluated, standard IS cannot be used because the weights $w(\theta_{i})$ in (3) are unavailable. Algorithm 1 presents the $\text{\rm IS}^{2}$ scheme for estimating the integral (2) when the likelihood $p(y|\theta)$ is estimated unbiasedly by $\widehat{p}_{N}(y|\theta)$, i.e. assuming that

Generate $\theta_{i}\stackrel{{\scriptstyle iid}}{{\sim}}g_{\text{\rm IS}}(\theta)$ and compute the likelihood estimate $\widehat{p}_{N}(y|\theta_{i})$.

Compute the weight $\widetilde{w}(\theta_{i})={p(\theta_{i})\widehat{p}_{N}(y|\theta_{i})}/{g_{\text{\rm IS}}(\theta_{i})}$.

with $\widetilde{g}_{\text{\rm IS}}(\theta,z)=g_{\text{\rm IS}}(\theta)g_{N}(z|\theta)$ an importance density on $\widetilde{\Theta}$. Let $(\theta_{i},z_{i})\stackrel{{\scriptstyle iid}}{{\sim}}\widetilde{g}_{\text{\rm IS}}(\theta,z)$, i.e. generate $\theta_{i}\sim g_{\text{\rm IS}}(\theta)$ and then $z_{i}\sim g_{N}(z|\theta_{i})$. It is straightforward to see that the estimator $\widehat{\varphi}_{\text{\rm IS}^{2}}$ defined in (4) is exactly an IS estimator of the integral defined in (5), with importance density $\widetilde{g}_{\text{\rm IS}}(\theta,z)$ and weights

This formally justifies Algorithm 1. Theorem 1 gives some asymptotic properties (in $M$) of $\text{\rm IS}^{2}$ estimators. Its proof is in the Appendix.

is finite for $h(\theta)=\varphi(\theta)$ and $h(\theta)=1$ for all $N$, then

where the asymptotic variance in $M$ for fixed $N$ is given by

If the conditions in (ii) hold, then $\widehat{\sigma^{2}_{\text{\rm IS}^{2}}(\varphi)}\stackrel{{\scriptstyle a.s.}}{{\longrightarrow}}\sigma^{2}_{\text{\rm IS}^{2}}(\varphi)$ as $M\to\infty$, for given $N$.

We note that both $\sigma^{2}_{\text{\rm IS}^{2}}(\varphi)$ and $\widehat{\sigma^{2}_{\text{\rm IS}^{2}}(\varphi)}$ depend on $N$, but, for simplicity, we do not show this dependence explicitly in the notation. Here, all the probabilistic statements, such as the almost sure convergence, must be understood on the extended probability space that takes into account the extra randomness occurring when estimating the likelihood. The result (iii) is practically useful, because (9) allows us to estimate the standard error of the IS2 estimator.

These convergence results are well known in the literature; see, e.g. Geweke, (1989).

The marginal likelihood $p(y)=\int_{\Theta}p(\theta)p(y|\theta)d\theta$ is a fundamental tool for Bayesian model choice (see, e.g., Kass and Raftery,, 1995). Estimating the marginal likelihood and its associated standard error accurately has proved difficult with standard MCMC techniques (Chib and Jeliazkov,, 2001; Perrakis et al.,, 2014). This problem is even more severe when the likelihood is intractable. This section and section 3.3 discuss how to estimate $p(y)$ and the associated standard error optimally, reliably and unbiasedly by $\text{\rm IS}^{2}$ when the likelihood is intractable.

The $\text{\rm IS}^{2}$ estimator of $p(y)$ is

with the samples $\theta_{i}$ and the weights obtained from Algorithm 1. The estimator of the variance of $\widehat{p}_{\text{\rm IS}^{2}}(y)$ is

The following theorem shows some properties of these $\text{\rm IS}^{2}$ estimators. Its proof is in the Appendix.

Let $M$ be the number of particles for $\theta$ and $N$ be the number of particles for estimating the likelihood. Under the assumptions in Theorem 1

2 The effect on importance sampling of estimating the likelihood

The results in the previous section show that it is straightforward to use importance sampling even when the likelihood is intractable but unbiasedly estimated. This section addresses the question of how much asymptotic efficiency is lost when working with an estimated likelihood. We follow Pitt et al., (2012) and Doucet et al., (2015) and make the following idealized assumption to make it possible to develop some theory. Proposition S2 of section S2.1 of the supplementary material justifies this assumption for panel data.

There exists a function $\gamma^{2}(\theta)$ such that the density $g_{N}(z|\theta)$ of $z$ is ${\cal N}(-\frac{\gamma^{2}(\theta)}{2N},\frac{\gamma^{2}(\theta)}{N})$, where ${\cal N}(a,b^{2})$ is a univariate normal density with mean $a$ and variance $b^{2}$.

If Assumption 2 holds for a fixed $\sigma^{2}$, then (7) becomes

These are the standard conditions for IS (Geweke,, 1989). The proof of this lemma is straightforward and omitted.

Recall that $\sigma^{2}_{\text{\rm IS}}(\varphi)/M$ and $\sigma^{2}_{\text{\rm IS}^{2}}(\varphi)/M$ are respectively the asymptotic variances of the IS estimators we would obtain when the likelihood is available and when the likelihood is estimated. We refer to the ratio $\sigma^{2}_{\text{\rm IS}^{2}}(\varphi)/\sigma^{2}_{\text{\rm IS}}(\varphi)$ as the inflation factor. Theorem 3 obtains an expression for the inflation factor, shows that it is independent of $\varphi$, greater than or equal to 1 and increases exponentially with $\sigma^{2}$. Its proof is in the Appendix.

Under Assumption 2 and the conditions in Theorem 1,

3 Optimally choosing the number of particles N𝑁N

From Theorems 1 and 3, the variance of the estimator $\widehat{\varphi}_{\text{\rm IS}^{2}}$ based on $M$ importance samples from $g_{\text{\rm IS}}(\theta)$ is approximated by

We now define the measure of the computing time of the IS2 method relative to that of the hypothetical IS method, to achieve the same precision, as

which is the product of the relative variance from (17) and the expected computational effort. It is straightforward to check that $\text{\rm CT}(\sigma^{2})$ is convex and minimized at

Section 4 discusses practical guidelines for selecting an optimal number of particles.

4 Optimal N𝑁N for estimating the marginal likelihood

Consequently, the relative variance of the two schemes is given by

Let $\sigma^{2}_{\min}(v)$ minimize $CT_{ML}(\sigma^{2})$ for a given $v$. The following proposition summarizes some properties of $\sigma^{2}_{\min}(v)$, where $\sigma^{2}_{\text{opt}}$ below is given by (19). Its proof is in section S4 of the supplementary material.

For any value of $v$, $\text{\rm CT}_{ML}(\sigma^{2})$ is a convex function of $\sigma^{2}$; therefore $\sigma^{2}_{\min}(v)$ is unique.

$\sigma^{2}_{\min}(v)$ increases as $v$ increases and $\sigma^{2}_{\min}(v)\longrightarrow\sigma^{2}_{\text{opt}}$ in (19) as $v\longrightarrow\infty$.

It can be readily checked that $\sigma^{2}_{\min}(v)$ is insensitive to $v$ in the sense that it is a flat function of $v$ for large $v$, with $\sigma^{2}_{\min}\to\sigma^{2}_{\text{opt}}=0.17$ as $v$ increases. See also Table 7 in section S6 of the supplementary material. Based on these observations, we advocate using $\sigma^{2}_{\text{opt}}$ in (19) as the optimal value of the variance of the log likelihood estimates in estimating the marginal likelihood.

A potential drawback with $\text{\rm IS}^{2}$, when estimating the posterior of $\theta$, but not for estimating the marginal likelihood, is that its performance depends on the proposal density $g_{\text{\rm IS}}(\theta)$ for $\theta$, which may be difficult to obtain in complex models. This section outlines the Mixture of $t$ by Importance Sampling Weighted Expectation Maximization (MitISEM) of Hoogerheide et al., (2012) approach we used in the article for designing efficient and reliable proposal densities for $\text{\rm IS}^{2}$. Section S5.1 of the supplementary material outlines a second approach based on annealed importance sampling for models with latent variables.

MitISEM constructs a mixture of $t$ densities for approximating the target distribution by minimizing the Kullback–-Leibler divergence between the target and the $t$ mixture, and can handle target distributions that have non-standard shapes such as multimodality and skewness. When the likelihood is available, the MitISEM can be used to design proposal densities that approximate the posterior accurately. It is natural to use an estimated likelihood when the likelihood is unavailable. We write $\widehat{p}_{N}(y|\theta)=\widehat{p}_{N}(y|\theta,u)$, with $u$ a fixed random number stream for all $\theta$. The target distribution in the MitISEM is $p(\theta)\widehat{p}_{N}(y|\theta,u)/p(y)$, which can be considered as the posterior $p(\theta|y,u)$ conditional on $y$ and the common random numbers $u$. Our procedure is analogous to using common random numbers $u$ to obtain simulated maximum likelihood estimates of $\theta$ (see, e.g., Gourieroux and Monfort,, 1995), except that we obtain a histogram estimate of the ‘posterior’ $p(\theta|u,y)\propto p(y|\theta,u)p(u)$. This ‘posterior’ is biased, but sufficiently good to obtain a good proposal density.

Practical guidelines for selecting an optimal N𝑁N

In the applications in Section 5, we use the time normalized variance (TNV) of an $\text{{IS}}^{2}$ estimator $\widehat{\varphi}_{\text{{IS}}^{2}}$ as a measure of its inefficiency, which we define as

Applications

The generalized multinomial logit (GMNL) model of Fiebig et al., (2010) specifies the probability of individual $i$ choosing alternative $j$ on occasion $t$ as

where $\beta_{i}=(\beta_{0i1},\ldots,\beta_{0iJ},\beta_{1i},\ldots,\beta_{Ki})^{\prime}$ and $X_{it}=(x_{1i1t},\ldots,x_{Ki1t},\ldots,x_{1iJt},\ldots,x_{KiJt})^{\prime}$ are the vectors of utility weights and choice attributes respectively. The GMNL model specifies the alternative specific constants as $\beta_{0ij}=\beta_{0j}+\eta_{0ij}$ with $\eta_{0ij}\sim\textrm{\rm N}(0,\sigma_{0j}^{2})$ and the attribute weights as

with $\eta_{ki}\sim\textrm{\rm N}(0,\sigma_{k}^{2})$ and $\zeta_{i}\sim\textrm{\rm N}(0,1)$. The expected value of the scaling coefficients $\lambda_{i}$ is one, implying that $\textrm{\rm E}(\beta_{ki})=\beta_{k}$.

When $\delta=0$ (so that $\lambda_{i}=1$ for all individuals) the GMNL model reduces to the mixed logit (MIXL) model, which we also consider in our analysis. The MIXL model captures heterogeneity in consumer preferences by allowing individuals to weight the choice attributes differently. By introducing taste heterogeneity, the MIXL specification avoids the restrictive independence of irrelevant alternatives property of the standard multinomial logit model (Fiebig et al.,, 2010). The GMNL model additionally allows for scale heterogeneity through the random variable $\lambda_{i}$, which changes all attribute weights simultaneously. Choice behavior in this model can therefore be more random for some consumers than others. The $\gamma$ parameter weights the specification between two alternative ways of introducing scale heterogeneity into the model.

The parameter vector is $\theta=(\beta_{01},\ldots\beta_{0J},\sigma_{0}^{2},\beta_{1},\ldots,\beta_{K},\sigma_{1}^{2},\ldots,\sigma_{K}^{2},\delta^{2},\gamma)^{\prime}$, while the vector of latent variables for each individual is $x_{i}=(\eta_{0i},\ldots,\eta_{Ki},\lambda_{i})$. The likelihood is therefore

where $y_{it}$ is the observed choice, $y=(y_{11},\ldots,y_{1T},\ldots,y_{I1},\ldots,y_{IT})^{\prime}$ and $p(y_{it}|x_{i})$ is given by the choice probability (23).

We apply our methods to the Pap smear data set considered by Fiebig et al., (2010), who used simulated maximum likelihood estimation. In this data set, $I=79$ women choose whether or not to have a Pap smear test ($J=2$) on $T=32$ choice scenarios. We let the observed choice for individual $i$ at occasion $t$ be $y_{it}=1$ if the woman chooses to take the test and $y_{it}=0$ otherwise. Table 1 lists the choice attributes and the associated coefficients. We set $\sigma_{5}^{2}=0$ in our analysis since we found no evidence of heterogeneity for this attribute beyond the scaling effect. The model (23) is identified by setting the coefficients as zero when not taking the test.

We specify the priors as $\beta_{01}\sim\textrm{\rm N}(0,100)$, $\sigma_{01}\propto(1+\sigma_{01}^{2})^{-1}$, $\beta_{k}\sim\textrm{\rm N}(0,100)$, $\sigma_{k}\propto(1+\sigma_{k}^{2})^{-1}$, for $k=1,\ldots,K$, $\delta\propto(1+\delta/0.2)^{-1}$, and $\gamma\sim\textrm{U}(0,1)$. The standard deviation parameters have half-Cauchy priors as suggested by Gelman, (2006).

1.2 Implementation details

We estimate the likelihood (24) by integrating out the vector of latent variables for each individual separately using different approaches for the MIXL and GMNL models. For the MIXL model, we combine the efficient importance sampling (EIS) method of Richard and Zhang, (2007) with the defensive sampling approach of Hesterberg, (1995). The importance density is the two component defensive mixture

where $h^{\text{EIS}}(x_{i}|y_{i1},\ldots,y_{iT})$ is a multivariate Gaussian importance density obtained using EIS. Following Hesterberg, (1995), including the natural sampler $p(x_{i})$ in the mixture ensures that the importance weights are bounded. We set the mixture weight as $\pi=0.5$. For the GMNL model, we follow Fiebig et al., (2010) and use the model density $p(x_{i})$ as an importance sampler. We implement this simpler approach for the GMNL model because the occurrence of large values of $\lambda_{i}$ causes the defensive mixture estimates of the log-likelihood to be pronouncedly right skewed in this case.

To obtain the parameter proposals $g_{\text{\rm IS}}(\theta)$, we use the MitISEM approach as described in section 3.5, which approximates the posterior of the two models as two component mixtures of multivariate Student’s $t$ distributions.

We implement two standard variance reduction methods at each IS stage: stratified mixture sampling and antithetic sampling. The first consists of sampling from each component at the exact proportion of the mixture weights. For example, when estimating the likelihood for the MIXL model we generate exactly $\pi N$ draws from the efficient importance density $h^{\text{EIS}}(x_{i}|y_{i1},\ldots,y_{iT})$ and $(1-\pi)N$ draws from $p(x_{i})$. The antithetic sampling method consists of generating pairs of perfectly negatively correlated draws from each mixture component, see, e.g., Ripley, (1987).

1.3 Posterior analysis

Table 11 presents selected posterior statistics estimated by the IS2 method under the optimal schemes (targeting $\sigma^{2}_{\text{opt}}=0.17$ for the MIXL model and $\sigma^{2}_{\text{opt}}=1$ for the GMNL model). We estimated the posterior distribution using $M=50,000$ importance samples for the parameters. We also estimated the Monte Carlo standard errors by bootstrapping the importance samples. We highlight two results. The MC standard errors are low in all cases, illustrating the high efficiency of the IS2 method for carrying out Bayesian inference for these two models. The estimates of the log of the marginal likelihood allow us to calculate a Bayes factor of approximately 20 for the GMNL model over the MIXL model, so that there is some evidence to suggest the presence of scale heterogeneity for this dataset.

1.4 Comparing IS2, pseudo marginal Metropolis-Hastings and SMC2

The IS2, PMMH, and SMC2 methods can tackle the same set of problems. This section compares the performance of IS2, PMMH and SMC2. We evaluate the first two methods under the same independent proposal for the parameter vector $\theta$ and likelihood estimation method, which leads to a direct comparison. To accurately measure the Monte Carlo efficiency of the three methods for the MIXL model, we run 500 independent replications of the IS2 algorithm using $M=5,000$ importance samples, 500 independent PMMH Markov chains with 5,000 iterations (after a burn-in of 1,000 iterations), and 500 independent replications of the SMC2 method with 5,000 particles. We focus on estimates of the posterior mean. We obtain the initial draw for the PMMH Markov chain by using the discrete IS2 approximation of the posterior. We implement the SMC2 method by sequentially updating the posterior for each individual. We estimate the likelihood contribution from each individual in the same way as for the IS2 for the PMMH method. We initialise the method implementing IS2 for the first patient in the dataset. Due to the large prior variances, we found it necessary to fit a Hoogerheide et al., (2012) proposal for this step. For the sequential updates, we follow Chopin et al., (2013) and use independent Metropolis-Hastings rejuvenation steps by fitting a multivariate Gaussian proposal to the current particle system.

In addition, we consider two improvements to the basic IS2 method: the Pareto smoothing (PS) method of Vehtari and Gelman, (2015) and Quasi-Monte Carlo (QMC) sampling of the parameters and trimmed means. The motivation for using the Pareto smoothing method is to stabilize the right tail of the distribution of IS2 weights, increasing the robustness of the method when there is noise in the weights because the likelihood is estimated. We implement QMC by replacing the standard normal draws used to generate samples from $q(\theta|y)$ by scrambled Sobol sequences converted to quasi-random $N(0,1)$ draws through the inverse CDF method.

Table 2 presents the MC mean squared errors for the posterior mean estimates and the MC standard error estimates, normalised by the performance of the PMMH method for the MIXL model. We approximate the posterior mean by combining the particles generated in all replications for the standard IS2 method. The results support the argument that the IS2 method leads to improved efficiency due to the use of independent samples and variance reduction methods. The table shows that the standard IS2 method leads to approximately 70% average reductions in the MC MSE compared to the PMMH method and the SMC2 implementation. We find that Pareto smoothing dramatically increases the accuracy of the MC standard error estimates, although it does not seem to have a large impact on the MSE of the posterior mean estimates for this model. In results not reported in the table, we found that Pareto smoothing leads to better behaved, near Gaussian, posterior mean estimates, so that we recommend this approach in practice to ensure robustness. The use of QMC leads to a modest incremental improvement of approximately 15% in the MSE of the posterior mean estimates. Regarding the computational costs, the IS2 and PMMH methods are equivalent in our setting, except for the additional cost of the burn-in period for the latter. Each replication of IS2 and PMMH required approximately 50 seconds in a machine equipped with an Intel Core i7-4790 CPU with 4 cores. In contrast, each SMC2 replication required approximately six minutes, highlighting the larger computational efficiency of IS2 and PMMH for off-line estimation.

Table 3 extends the analysis to the more challenging GMNL specification. We consider only the IS2 and PMMH methods in this case, as the the SMC2 method was prohibitively costly. The results are similar to the ones in Table 3: the IS2 method leads to 78% reductions in MSE on average compared to PMMH. When considering only the variations of the $\text{\rm IS}^{2}$ method, the results show that for this model the use of Pareto smoothing substantially improves the performance of the method, leading to 57% lower MSEs compared to the standard implementation. As before, the results also suggest that the estimates of the MC variance within each of the replications are substantially more accurate under IS2, especially with Pareto smoothing.

2 Stochastic volatility

We consider the univariate two-factor stochastic volatility (SV) model with leverage effects

where the return innovations are i.i.d. and have standardized Student’s $t$ distributions with $\nu$ degrees of freedom and $1>\phi_{1}>\phi_{2}>-1$ for stationarity and identification. This SV specification incorporates the most important empirical features of a volatility series, allowing for fat-tailed return innovations, a negative correlation between lagged returns and volatility (leverage effects) and quasi-long memory dynamics through the superposition of AR(1) processes (see, e.g., Barndorff-Nielsen and Shephard,, 2002). The two volatility factors have empirical interpretations as long term and short term volatility components.

We use the following priors: $\nu\sim\textrm{U}(2,100)$, $c\sim\textrm{\rm N}(0,1)$, $\phi_{1}\sim\textrm{U}(0,1)$, $\sigma_{1}^{2}\sim\textrm{IG}(2.5,0.0075)$, $\rho_{1}\sim\textrm{U}(-1,0)$, $\phi_{2}\sim\textrm{U}(0,\phi_{1})$, $\sigma_{2}^{2}\sim\textrm{IG}(2.5,0.045)$, and $\rho_{2}\sim\textrm{U}(-1,0)$, and estimate the two factor SV model using daily returns from the S&P 500 index obtained from the CRSP database. The sample period runs from January 1990 to December 2012, for a total of $5,797$ observations. As in the previous example, we obtain the parameter proposal $g_{\text{\rm IS}}(\theta)$ using the MitISEM method of Hoogerheide et al., (2012), where we allow for two mixture components. To estimate the likelihood, we apply the particle efficient importance sampling (PEIS) method of Scharth and Kohn, (2016), which is a sequential Monte Carlo procedure for likelihood estimation that uses the EIS algorithm of Richard and Zhang, (2007) to obtain a high-dimensional importance density for the states. We implement PEIS using antithetic variables as described in Scharth and Kohn, (2016).

2.2 Optimal number of particles

To investigate the accuracy of the likelihood estimates for the stochastic volatility model, we replicate the Monte Carlo exercise of section S8 for this example. Table 4 displays the average variance, skewness and kurtosis of the log-likelihood estimates based on $N=10$ and $N=20$ particles, as well as the computing times. We base the analysis on 1,000 draws from the proposal for $\theta$ and 200 independent likelihood estimates for each parameter value.

The results show that the PEIS method is highly accurate for estimating the likelihood of the SV model, despite the large number of observations. When $N=10$, the inflation factor of (16) is approximately 1.01, indicating that performing IS with the PEIS likelihood estimate has essentially the same efficiency as if the likelihood was known. On the other hand, the log-likelihood estimates display positive skewness and excess kurtosis and are clearly non-Gaussian, suggesting that the theory of section 3.3 only holds approximately for this example.

We approximate the optimal number of particles by assuming that the variance of log-likelihood estimates is constant across different values of $\theta$, as it is not possible to use the jackknife method for estimating the variance of log-likelihood estimators based on particle methods. Based on the average variance of the log-likelihood estimate with $N=20$, we estimate the asymptotic variance of the log-likelihood estimate to be $\widehat{\gamma^{2}}(\theta)=0.1$ on average. The last row of the table allows us to calculate that the overhead of estimating the likelihood is $\tau_{0}=1.051$ seconds, while the computational cost of each particle is $\tau_{1}=0.018\times 10^{-1}$ seconds. Assuming normality of the log-likelihood estimate and that $\gamma^{2}(\theta)=\bar{\gamma}^{2}=0.1$ for all $\theta$, the optimal number of particles is $N_{\textrm{opt}}=8$.

Figure 1 plots the predicted relative time normalized variance $\text{TNV}(M,N)/\text{TNV}(M,N_{\text{opt}})$ as a function of $N$. The figure suggests that the efficiency of $\text{\rm IS}^{2}$ for this example is fairly insensitive to the number of particles at the displayed range $N\leq 40$, with even the minimal number of particles $N=2$ (including an antithetic draw) leading to a highly efficient procedure.

Table 5 displays estimates of the relative variances for the posterior means based on $M=50,000$ importance samples for the parameters. We estimate the Monte Carlo variance of the posterior statistics by bootstrapping the importance samples. The results are consistent with Figure 1 and indicate that the efficiency of the IS2 method is insensitive to the number of particles in this example, with no value of $N$ emerging as clearly optimal in the table. Based on this result, we recommend a conservative number of particles in practice as there is no important efficiency cost in choosing $N$ moderately above the optimal theoretical value.

2.3 Posterior analysis

Table 6 presents estimates of selected posterior distribution statistics estimated by the IS2 method. We estimated the posterior distribution using $M=50,000$ importance samples for the parameters and $N=20$ particles to estimate the likelihood. As before, we estimate the Monte Carlo standard errors by bootstrapping the importance samples. The results show that IS2 leads to highly accurate estimates of the posterior statistics and that the short term volatility component is almost entirely driven by leverage effects.

Conclusions

This article proposes the IS2 method for Bayesian inference when the likelihood is intractable but can be estimated unbiasedly. Its advantages are that it gives accurate estimates of posterior moments and the marginal likelihood and estimates of the standard errors of these estimates. We note that if the model is estimated by IS2, then the marginal likelihood estimate and the standard error of the estimator are obtained automatically. However, even if the model is estimated by MCMC, IS2 can be used effectively to estimate the marginal likelihood and the corresponding standard error. The article studies the convergence properties of the $\text{\rm IS}^{2}$ estimators. It examines the effect of estimating the likelihood on Bayesian inference and provides practical guidelines on how to optimally select number of particles to estimate the likelihood in order to minimize the computational cost. The applications illustrate that the $\text{\rm IS}^{2}$ method can lead to fast and accurate posterior inference when optimally implemented, and demonstrate that the theory and methodology presented in this paper are useful for practitioners who are working on models with an intractable likelihood.

Supplementary material

The paper has an online technical supplement that contains: (a) Some large sample results for panel data about the estimator of the likelihood and static and dynamic estimators of the optimal variance of the log-likelihood; (b) proofs of all propositions; (c) Discussion of an alternative importance sampling density $g_{IS}(\theta)$ and the definition of equivalent sample size (ESS).

Appendix A Proofs

Proof of the unbiasedness of $\widehat{p}_{\text{\rm IS}^{2}}(y)$ is straightforward. The assumptions in Theorem 1 ensure that $\widetilde{w}(\theta_{i})$’s are i.i.d with a finite second moment. The results of (i) and (iii) then follow. To prove (ii), we have

Appendix S1 Prelminaries

This supplement contains results that complement those in the main paper. References to the main paper are of the form section 2, equation (2) and assumption 2, etc, whereas for the supplementary material we use section S2, equation (S2) and assumption S2, etc.

Appendix S2 Large sample results for panel data

This section studies the large sample (in $n$) properties of the estimator of the likelihood for panel data which is described in section 2, and obtains two results assuming that the number of particles is $N=O(n)$. First, we justify Assumption 2 by showing the asymptotic normality of the error $z$ in the log-likelihood estimator. Second, we obtain the convergence rates of two estimators of the optimal variance, $\sigma^{2}_{\rm opt}$, of the error $z$ in the log likelihood estimator. We show that the dynamic estimator of $\sigma^{2}_{\rm opt}$ that chooses the number of particles or samples $N$ depending on $\theta$ is much more efficient for large $n$ than the static estimator of $\sigma^{2}_{\rm opt}$ that chooses a single number of particles $N$ for all $\theta$. Using the notation in section 2, the error in the log-likelihood estimator is

$\omega_{i}(\alpha_{i},\theta)=p(y_{i}|\alpha_{i},\theta)p(\alpha_{i}|\theta)/h_{i}(\alpha_{i}|y,\theta)$, and $\alpha_{i}^{(j)}\stackrel{{\scriptstyle iid}}{{\sim}}h_{i}(\alpha_{i}|y,\theta)$.

All equations, sections etc in this supplement are prefixed by S, e.g. equation (S1). Equations, sections etc that are not so prefixed refer to the main paper.

Assumption S3 is used to prove Proposition S2.

We assume that for all $\theta$ and $i=1,\dots,n$,

$K_{2}^{-1}<p(y_{i}|\theta)\leq K_{2}$, where $K_{2}>1$, so that $p(y_{i}|\theta)$ is bounded and also bounded away from zero.

Proposition S2 gives the following asymptotic results for the error $z$. Its proof is in the section S4.

Proposition S2 justifies Assumption 2, with $\gamma^{2}(\theta)=n\psi(n,\theta)$.

The convergence of $z$ to normality is fast, because $N=O(n)$ and from Part (vii) of Lemma S6 and (S25), $z$ is a sum of $n$ terms, each of which converges to normality.

Appendix S3 Estimating the optimal number of particles statically and dynamically

Suppose the following five conditions hold.

where $\psi(\theta,n)$ is $O(1)$ in $n$ uniformly in $\theta\in\Theta$.

A consistent estimator $\widehat{\psi}_{N}(\theta,n)$ (in $N$) of $\psi(\theta,n)$ is available such that, uniformly in $\theta$,

We have a consistent estimator (in $n$) $\widehat{\theta}_{n}$ of $\overline{\theta}$.

$\psi(\theta,n)$ is a continuously differentiable function of $\theta\in\Theta$.

We note that Conditions (i) and (ii) are justified for the IS estimator by Parts (i) and (ii) of Proposition S2. The dynamic estimator is preferred because it converges faster and does not require any assumption on the behavior of the posterior $\pi(\theta)$. The efficiency of the dynamic estimator is shown empirically in section S8 of the article. For part (iv), we note that any consistent estimator, such as a method of moments estimator, may be employed. However, the dynamic estimator is more expensive than the static estimator as it also requires computing $\widehat{\psi}_{N_{\widehat{\theta}_{n}}}(\theta,n)$.

Appendix S4 Proofs of the Propositions

Denote $a=\tau_{0}$, $b=\tau_{1}\bar{\gamma}^{2}$, and $x=\sigma^{2}$. Then

We write $f(x)=f_{1}(x)+bf_{2}(x)$ with $f_{1}(x)=a((v+1)e^{x}-1)$ and $f_{2}(x)=((v+1)e^{x}-1)/x$. As $f_{1}(x)$ is convex, it is sufficient to show that $f_{2}(x)$ is convex. To do so, we will prove that $f_{2}^{\prime\prime}(x)>0$ for $x>0$, because a differentiable function is convex if and only if its second derivative is positive (see, e.g., Bazaraa et al.,, 2006, Chapter 3).

The first term is clearly positive. We consider the second term in the square brackets

Note that $1-e^{-x}>x-\frac{x^{2}}{2}$ for $x>0$ and so

This establishes that $f_{2}^{\prime\prime}(x)>0$ for all $x>0$.

To prove (ii), for any fixed $v$, let $x_{\min}(v)$ be the minimizer of $f(x)$. By writing

with $\text{\rm CT}(x)$ defined by (18), we can see that $f(x)$ is driven by the factor $(v+1)\text{\rm CT}(x)$ as $v\to\infty$. Hence $x_{\min}(v)$ tends to the $\sigma^{2}_{\text{opt}}$ in (19) that minimizes $\text{\rm CT}(x)$, as $v\to\infty$.

Because $f^{\prime}(x_{\min}(v))=0$ for any $v>0$,

By taking the first derivative of both sides of (S34), we have

This follows that $x_{\min}(v)$ is an increasing function of $v$. Furthermore, $x_{\min}^{\prime}(v)\to 0$ as $v\to\infty$, which establishes (iii). ∎

To obtain Proposition S2 we first obtain some preliminary results. Define,

so that $\varepsilon_{N,i}$ is a sum of the i.i.d. random variables.

Suppose Assumptions 1 and S3 hold. Then, for all $i=1,\dots,n$, and $\theta$,

$1/K_{5}<\sigma_{i}(\theta)^{2}<K_{5}$, where $K_{5}>1$ is a constant.

$\varepsilon_{N,i}\overset{d}{\rightarrow}\mathcal{N}(0,1)$ as $N\rightarrow\infty$.

We write the error $z$ in (S25) in terms of the $\varepsilon_{N,i}$ and expand as a third order Taylor series approximation with a remainder term,

The remainder term is $R_{N,i}(\theta)\leq\sigma_{i}^{4}(\theta)\varepsilon_{N,i}^{4}/4N^{2}$.

Part (ii) follows from Lemma S6. Part (iii) follows from Part (ii). To obtain Part (iv), it is sufficient to prove the central limit theorem for $z_{1}$, ,which is a sum of independent random variables. Now, by Parts (iv) and (v) of Lemma S6,

The central limit theorem now follows from Theorem 5, p. 194, of Stirzaker and Grimmett, (2001) ∎

We note that as $N$ is $O(n)$, the central limit operates very quickly in $n$ upon the expression for $z_{1}$ in (S40) as the summation is over an increasing number of terms, each of which satisfies a central limit.

As $N_{S}$ is chosen proportional to $n$ then consequently, $N_{\widehat{\theta}_{n}}$ is of order $n$. Using assumption (i), the achieved variance at any ordinate $\theta$ with $N_{\widehat{\theta}_{n}}$ particles is

Using Assumptions (iii) to (v), we obtain from a Taylor expansion of $\psi(\theta,n)$ around $\widehat{\theta}_{n}$ in (S41) that

This establishes part A. For part B, we may write

since $N_{\widehat{\theta}_{n}}$ is of order $n$. Hence, the achieved variance at a given $\theta$ and $n$ is given by

Appendix S5 Discussion of complementary methodology

The second approach for constructing the density $g_{IS}$ is based on the annealed importance sampling for models with latent variables (AISEL) method proposed by Duan and Fulop, (2015) for state space models and independently by Tran et al., (2014) for general models with intractable likelihoods, which is based on the annealed importance sampling (AIS) approach of Neal, (2001), but where the likelihood is now estimated. We advocate taking a small number of steps at a high temperature to obtain $\theta$ particles to train the $g_{\rm IS}(\theta)$ proposal. The overall AISEL proposal is far more expensive than IS2, but taking a few steps of AISEL at high temperatures requires relatively few particles and so is relatively cheap. The AISEL method produces a sequence of weighted samples of $\theta$’s which are first drawn from an easily-generated distribution, and then moved towards the target distribution through Markov kernels. A very good proposal density $g_{\text{\rm IS}}(\theta)$ that is heavier tailed than the target density can then be obtained by fitting a mixture of $t$ densities to the resulting empirical density for $\theta$.

S5.2 Effective sample size

The efficiency of the proposal density $g_{\text{\rm IS}}(\theta)$ in the standard IS procedure (3) is often measured by the effective sample size defined as (Liu,, 2001, p.35)

If the number of particles $N(\theta)$ is tuned to target a constant variance $\sigma^{2}$ of $z$, then (S43) enables us to estimate $\text{ESS}_{\text{\rm IS}}$ as if the likelihood was given. This estimate is a useful measure of the efficiency of the proposal density $g_{\text{\rm IS}}(\theta)$ in the $\text{\rm IS}^{2}$ context.

This table is referred to in section 3.4 and shows that the ratio $\text{\rm CT}_{\widehat{p}_{\text{\rm IS}^{2}}(y)}(\sigma^{2}_{\text{opt}})/\text{\rm CT}_{\widehat{p}_{\text{\rm IS}^{2}}(y)}(\sigma^{2}_{\min})$ is insensitive to $v$ and so justifies using $\sigma^{2}_{\text{opt}}$ instead of $\sigma^{2}_{\text{min}}$ when estimating the marginal likelihood.

Appendix S7 Empirical performance of IS2 for panel data.

Appendix S8 Empirical comparison of the static and dynamic estimators of the optimal number of particles

This section investigates the empirical performance of the static and dynamic estimators of the optimal number of particles defined in section S3, for the panel data models discussed in section 5.1 of the main paper.

We generated 1,000 draws for the parameter vector $\theta$ from $g_{\text{\rm IS}}(\theta)$. For each parameter combination, we obtained 100 independent likelihood estimates using the same fixed number of particles for all individuals in the panel and by targeting a log-likelihood variance of 1 and 0.5 for the MIXL model using the static approach (fixed $N$) to determining $N$ and targeting a variance of 0.5 for each $\theta$ using the dynamic approach. Similarly, for the GMNL model, the static approach with fixed $N=8000$ and $N=16000$ particles targets a variance of 2.0 and 1.0 respectively, while the dynamic approach targets a variance of 1.0 for each $\theta$.

Table 8 also provides estimates of the required quantities for calculating the optimal $\sigma^{2}$ in (19). For the MIXL model, since the average variance when $N=24$ is 1.068, we estimate $\overline{\gamma^{2}}=1.068\times 24=25.63$. As the likelihood evaluation time using $N$ particles is determined as $\tau_{0}+\tau_{1}N$, from the computing times at $N=24$ and $N=48$, we obtain $\tau_{0}=0.067$ and $\tau_{1}=8.97\times 10^{-5}$ seconds. Using (19), we conclude that $\sigma^{2}_{\text{opt}}$ for the MIXL model is approximately 0.17. The optimal variance $\sigma^{2}_{\text{opt}}=1$ for the GMNL case because $\tau_{0}=0$ as there is no overhead cost in obtaining the state proposal.

To compare different implementations of IS2, we generate $M=50,000$ draws from the parameter proposal $g_{\text{\rm IS}}(\theta)$ and run the IS2 by targeting different variances of the log-likelihood estimates. The target variances are $\sigma^{2}=0.05$, 0.2, 0.75, 1 for the MIXL model and $\sigma^{2}=0.5$, 1, 1.5, 2 for the GMNL model. We estimate the Monte Carlo (MC) variances of the posterior means under each implementation by the bootstrap method. Tables 9 and 10 show the relative estimated MC variances, the total actual computing time in minutes, and the relative estimated time normalized variance (22) for the MIXL and GMNL models respectively. We also report for reference the theoretical relative MC variances and TNV (see (17) and (22)).

The results show that the relative estimated MC variances are on average close to their theoretical values, even though there is a large degree of variability in the estimates across the parameters (especially for higher values of $\sigma^{2}$). The estimated relative TNV are approximately consistent with their theoretical values for the optimal precision of the log-likelihood estimates. Finally, $\widehat{\text{\rm ESS}}_{\text{\rm IS}^{2}}/\widehat{\text{\rm ESS}}_{\text{\rm IS}}\approx 0.84$ for the MIXL model so that IS2 is only 16% less efficient than IS. The corresponding figure is 0.37 for the GMNL model, so here IS2 is 63% less efficient than IS as we used a very simple importance sampler.

References