{'fraction_non_alphanumeric': 0.04941998668821907, 'fraction_numerical': 0.01887420367024817, 'mean_word_length': 4.090634075508229, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 19, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On the existence of moments for high dimensional importance sampling

Theoretical results for importance sampling rely on the existence of certain moments of the importance weights, which are the ratios between the proposal and target densities. In particular, a finite variance ensures square root convergence and asymptotic normality of the importance sampling estimate, and can be important for the reliability of the method in practice. We derive conditions for the existence of any required moments of the weights for Gaussian proposals and show that these conditions are almost necessary and sufficient for a wide range of models with latent Gaussian components. Important examples are time series and panel data models with measurement densities which belong to the exponential family. We introduce practical and simple methods for checking and imposing the conditions for the existence of the desired moments. We develop a two component mixture proposal that allows us to flexibly adapt a given proposal density into a robust importance density. These methods are illustrated on a wide range of models including generalized linear mixed models, non-Gaussian nonlinear state space models and panel data models with autoregressive random effects.


Introduction
This paper develops robust high-dimensional importance sampling methods for estimating the likelihood of statistical and econometric models with latent Gaussian variables. We propose computationally practical methods for checking and ensuring that the second (and possibly higher) moments of the importance weights are finite. The existence of particular moments of the importance weights is fundamental for establishing the theoretical properties of estimators based on importance sampling. A central limit theorem applies to the importance sampling estimator if the second moment of the weights is finite (Geweke, 1989). The Berry-Esseen theorem (Berry, 1941;Esseen, 1942) specifies the rate of convergence to normality when the third moment exists. In a more general case where the weights can have local dependence and may not be identically distributed, if the n th moment of the weights exists, Chen and Shao (2004) show that the rate of convergence to normality is of the order of O(S −(n−2)/4 ) with S the number of importance samples. This implies that the higher the moments that exist, the faster the rate of convergence. Even though infinite variance may in some cases be due to a region of the sampling space that has no practical impact for Monte Carlo sampling (Owen and Zhou, 2001), this problem can result in highly inaccurate importance sampling estimates in a variety of applications; see Robert and Casella (2005) and Richard and Zhang (2007) for examples.
In many models, e.g. non-Gaussian nonlinear state space models and generalized linear mixed models, the likelihood involves analytically intractable integrals. These integrals are often estimated by importance sampling. A popular approach for estimating the parameters in such models is the simulated maximum likelihood (SML) method (see, e.g., Gourieroux and Monfort, 1995). SML first estimates the likelihood by importance sampling and then maximizes the estimated likelihood. A second approach for estimating such models is by Bayesian inference using Markov chain Monte Carlo simulation. If the likelihood is estimated unbiasedly, then the Metropolis-Hastings algorithm with the likelihood replaced by its unbiased estimator is still able to sample exactly from the posterior. See, for example, Andrieu and Roberts (2009) and Flury and Shephard (2011). Pitt et al. (2012) show that it is crucial for this approach that the variance of the estimator of the log-likelihood exists, which is guaranteed by a finite second moment of the weights when the likelihood is estimated by importance sampling. Both the SML and MCMC methods require an efficient importance sampling estimator of the likelihood.
This paper introduces methods that directly check the existence of any particular moment of the importance weights when using a possibly high-dimensional Gaussian importance density. If the Gaussian importance density does not possess the required moments, our results make it straightforward to modify the original Gaussian importance density in such a way that the desired moment conditions hold. This approach is extended by developing a mixture importance density which is a combination of that robust modified Gaussian density with any other importance density. We prove that this mixture importance density satisfies the desired moment conditions, while providing substantial flexibility to design practical and efficient importance densities. The only requirement for the validity of these methods is that the log measurement density is concave or bounded linearly from above as a function of the latent variables. Our method contrasts with previous approaches which rely on statistical testing. See for example Monahan (1993) and Koopman et al. (2009), who develop diagnostic tests for infinite variance based on extreme value theory.
We also develop specific results that allow us to efficiently check and impose the required moment conditions for nonlinear non-Gaussian state space models. Gaussian importance samplers are extensively applied in this setting. Some examples include non-Gaussian unobserved components time series models as in Durbin and Koopman (2000), stochastic volatility models in Liesenfeld and Richard (2003), stochastic conditional intensity models in Bauwens and Hautsch (2006), stochastic conditional duration models in Bauwens and Galli (2009), stochastic copula models in Hafner and Manner (2012) and dynamic factor models for multivariate counts in Jung et al. (2011).
We illustrate the new method in a simulation study and in an empirical application to Bayesian inference for a Poisson panel data model with autoregressive random effects. We consider the Shephard and Pitt (1997) and Durbin and Koopman (1997) (SPDK) importance sampling algorithm, a commonly used approach based on local approximation techniques. We find that the SPDK method leads to infinite variance for this problem. The results show that the mixture approach for imposing the finite variance condition can provide an insurance against poor behavior of the SPDK method for estimating the likelihood at certain parameter values (for example, at the tail of the posterior density), even though the standard method performs nearly as well in most settings. This result suggests that the infinite variance problem typically has low probability of causing instability for this example. We also show that the mixture importance sampler leads to more efficient estimates compared to Student's t importance sampler with the same mean and covariance matrix as the SPDK Gaussian sampler.
Section 2 provides the background to importance sampling and establishes the notation. Section 3 presents the main theoretical results. Section 4 discusses the robust importance sampling methods for three popular classes of models including generalized linear mixed models, non-Gaussian nonlinear state space models and panel data models. Section 5 provides several illustrative examples.

Importance sampling
In many statistical applications we need to estimate an analytically intractable likelihood of the form where y is the observed data, α is a d × 1 vector of latent variables, and ψ is a fixed parameter vector at which the likelihood is evaluated. The density p(y|α; ψ) is the conditional density of the data given the latent variable α and p(α; ψ) is the density of α. These densities usually depend on a parameter vector ψ. In some classes of models such as state space models, (1) can be a very high-dimensional integral. In some other classes such as generalized linear mixed models and panel data models, the likelihood decomposes into a product of lower dimensional integrals. Since the parameter vector ψ has no impact on the mathematical arguments that follow, we omit to show this dependence, except that we still write L(ψ) to explicitly indicate the dependence of the likelihood on ψ.
. If E g [ω(α) 2 ] < ∞, a central limit theorem holds, i.e. W d → N(0, 1) (see, e.g., Geweke, 1989). Furthermore, if E g [ω(α) 3 ] < ∞, the rate of convergence of L to normality can be determined by the Berry-Esseen theorem (Berry, 1941;Esseen, 1942), which states that C is a finite constant and Φ(z) is the standard normal cdf. Our article first obtains conditions for the existence of the required moments of the weights ω(α) when using a Gaussian importance density g(α) = N(α|µ * , Q * −1 ). The results are then extended to the case of a mixture importance density.

General properties of importance weights
We are concerned with the existence of the n th moment E g [ω(α) n ] of the importance weights ω(α) when using a Gaussian importance density g(α) = N(α|µ * , Q * −1 ). In this paper, n can be any positive number. Of special interest is cases where n is a positive integer and we refer to n = 1 as the fist moment, n = 2 as the second moment and so forth. Let l(α) = log p(y|α). We write the log of the importance weights as where the constant c does not depend on α. In this paper, for a square matrix A, the notation A > 0 (A < 0) means that A is a positive definite (negative definite). We obtain the following general result on the moments of ω(α).
Proposition 1. Suppose that there exists a constant scalar k, a vector ξ and a symmetric matrix Q > 0 such that Proof. See Appendix A.
When n = 0 condition (4) requires that the inverse covariance matrix Q * of the importance density g(α) is positive definite. For n = 1, condition (4) holds because Q > 0 by assumption. The requirement for a finite variance of the weights, corresponding to the case n = 2, requires that 2Q−Q * > 0. In practice the bounding assumption (3) can be difficult to verify for general models. However, for models in which the density of the latent vector p(α) is Gaussian, verifying this assumption is straightforward. Also, when p(α) is Gaussian it is relatively straightforward to check the existence condition (4) and to impose conditions on Q * to ensure the existence of specific moments of the weights. The rest of this paper will therefore focus on models for which the latent density p(α) is itself multivariate Gaussian N(α|µ, Q −1 ) with Q > 0.
(ii) Suppose that for at least one j (j = 1, ..., d), Then Proof. See Appendix A.
Part (i) of Proposition 2 indicates that when the log measurement density is concave in the latent variable α, a sufficient condition for the existence of the first n moments is that Q * − n(Q * − Q) > 0. Under the additional assumption of part (ii), the condition is also necessary. For many models, it is possible to verify these Assumptions. Assumption (5) covers all exponential family models with a canonical link to the latent variable α, i.e. the density p(y|α) is of the form where η(α) = z ′ α+c, z is a vector of covariates, c is a constant that does not depend on α, and ̟ > 0 is a dispersion parameter (see Section 4). Becauseb(η) = ∂ 2 b(η)/∂η 2 > 0 by the property of the exponential family, It follows that l(α) is a concave function in α, because a differentiable function is concave if and only if its Hessian matrix is negative definite (see, e.g., Bazaraa et al., 2006, Chapter 3). Similarly, we can show that the stochastic volatility model (see, e.g., Ghysels et al., 1996) satisfies (5) with Gaussian or Student t errors in the observation equation. The observation equation for the univariate stochastic volatility model with Gaussian errors is y t = exp(α t /2)ǫ t , t = 1, ..., T , with ǫ t ∼ N(0, 1). Hence, It is straightforward to show that the Hessian matrix of l(α) is negative definite, thus l(α) is concave. The concavity of l(α) in the case with t errors can be shown similarly. Assumption (6) is also satisfied by many popular models. For example, in the Poisson model where b(η(α)) = exp(z ′ α+c) and in the binomial model where b(η(α)) = log(1 + e z ′ α+c ), we can easily check that (6) holds. For the stochastic volatility model with l(α) given in (7), (6) holds because lim α j →+∞ (l(α)/α 2 j ) = 0 for any j.

A general method for checking and imposing the existence condition
This section presents a general method for checking and imposing the existence condition (4). A more efficient method that exploits the structure of non-Gaussian nonlinear state space models is presented in Section 4.2. Suppose that g(α) = N(α|µ * , Q * −1 ) is an importance density available in the literature, e.g. one obtained by the Laplace method. It is straightforward to check the positive definiteness of the matrix Q * − n(Q * − Q) by verifying that all its eigenvalues are positive or by using Sylvester's criterion, which is discussed in Section 4.2.1.
If Q * fails the existence condition, we modify Q * to construct a matrix Q such that Q − n( Q − Q) > 0 by first decomposing nQ = AA ′ using the Cholesky decomposition.
The new importance density g(α) = N(α|µ * , Q −1 ) then satisfies the existence condition (4). For SML, it is essential that the likelihood estimator L(ψ) is continuous, and desirably, differentiable with respect to the model parameters ψ. The function λ j = λ j (λ j ) in (8) is not a continuous function of λ j , so λ j (ψ) = λ j (λ j (ψ)) is not continuous in ψ. Appendix C presents a modification of (8) in which λ j is a continuous and differentiable function of λ j .
The principle behind the above method is that the difference Q − Q * is often small, so that g(α) still provides a good fit to the target density, while the condition for the existence of the n th moment is guaranteed to hold. Intuitively, only directions of Q * along which the density N(α|µ * , Q * −1 ) has light tails are modified; the other directions are unchanged. The resulting density N(α|µ * , Q −1 ) has heavier tails than the original density N(α|µ * , Q * −1 ).
However, we have seen that the observation density p(y|α) does not matter for the existence of moments under the assumptions of Proposition 2, so that the importance density g(α) may be a poor approximation to p(y|α)p(α) near its mode. To overcome this problem, we establish the following result. Let f (α) = p(y|α)p(α) and sup(f ) = {α : f (α) = 0} be the support of f . Proposition 3. Suppose that g 1 (α) is a density such that sup(g 1 ) ⊇ sup(f ), and for some n ≥ 0, For any density g 2 (α), consider the mixture importance density g(α) = πg 1 (α) + (1 − π)g 2 (α) with 0 < π < 1. Then, for all m ≤ n, The proposition suggests an approach to combine the importance density provided by the standard methods with a second importance density which by itself ensures the existence of the required moments. This result gives substantial flexibility in that it is useful in practice to put very little weight on the heavy density component g 1 (α), while leaving the rest of the mass for the lighter tailed component. The existence of the moments will still be entirely governed by the heavier term g 1 (α), but the use of the mixture proposal may lead to lower variance since g 2 (α) is designed to approximate the target density accurately.
In all the models and examples considered in this article, Assumption (9) is satisfied with g 1 (α) = N(α|µ * , Q −1 ) and it is obvious that sup(g 1 ) ⊇ sup(f ) as sup(g 1 ) = R d . The density g 2 (α) can be any density that we can sample from. In practice, we would like to choose g 2 such that the weight ω(α) has a small variance, while the desired moments are theoretically guaranteed to exist. In this paper we choose g 2 (α) = N(α|µ * , Q * −1 ), and set π = 0.1 based on some experimentation. We refer to this approach as the n th -moment constrained mixture importance sampler or n th -IS for short.

Generalized linear mixed models
A popular class of models that typically needs importance sampling for likelihood estimation is generalized linear mixed models (GLMM) (see, e.g., Jiang, 2007). Consider a GLMM with data y i = (y i1 , ..., y in i ) ′ , i = 1, ..., m. Conditional on the random effects α i , the observations y ij are independently and exponentially distributed as .., m, where x ij , z ij are p-and u-vectors of covariates. For simplicity, we assume that the dispersion parameter ̟ is known; the cases with an unknown ̟ require only a small modification to the procedure described below. The random effects are assumed to have a normal distribution α i iid ∼ N(0, Q −1 ). The parameters of interest are ψ = (β, Q). The density of y conditional on ψ and α is so that the likelihood is It is often difficult to estimate the model parameters ψ because the likelihood (10) involves analytically intractable integrals L i (ψ). A popular approach for estimating ψ is the SML method (see, e.g., Gourieroux and Monfort, 1995). SML first estimates the integrals in (10) by importance sampling and then maximizes the estimated likelihood over ψ. It is necessary to use common random numbers so that the resulting estimator L(ψ) is smooth in ψ. Another approach for estimating ψ is by using MCMC with the likelihood (10) replaced by its unbiased estimator (Andrieu and Roberts, 2009;Flury and Shephard, 2011;Pitt et al., 2012). Both the SML and MCMC methods require an efficient estimator of the likelihood (10). The robust importance sampling approach for estimating the integrals L i (ψ) can be carried out as follows.
where C is the constant term independent of α i . Then whereḃ(η) andb(η) are the first and second derivatives of b(η) and are understood com- It is straightforward to use Newton's method to find a maximizer α * i of F i (α i ), because both the first and second derivatives of We can now construct the n th -moment constrained mixture importance sampler as described in Section 3.1 to ensure the existence of the first n moments of the weights. It is therefore straightforward to carry out importance sampling in GLMM with the required moments of the weights guaranteed to exist. We note that l i (α i ) = log p(y i |β, α i ) is concave in α i for all GLMMs with the canonical link (see Section 3), so that Proposition 2(i) applies.
A common practice is to use a t importance density with mean α * i , scale matrix Q * i −1 and small degrees of freedom ν, although the moments of the weights are not theoretically guaranteed to exist. Section 5.4 compares empirically the performance of such a t importance density to our constrained mixture importance density.

Nonlinear non-Gaussian state space models
Consider the following class of nonlinear non-Gaussian state space models where y t is the k ×1 observation vector, θ t is the p×1 signal vector, α t is the m×1 state vector, and Z t is the p × m selection matrix; the m × 1 constant vector d, the m × m transition matrix Φ and the m × m covariance variance matrix Σ jointly determine the dynamic properties of the model. The parameter vector ψ contains the unknown coefficients in the observation density and in the system matrices. A thorough overview of the models and methodology is provided in Durbin and Koopman (2001).
The likelihood for the model is with p(α) = p(α 1 ) T t=2 p(α t |α t−1 ) and p(y|α) = T t=1 p(y t |α t ). The latent states are Gaussian, with p(α) = N(α|µ; Q −1 ), where Q is a block tridiagonal matrix because the states α t follow a first order autoregressive process. A closed form representation of Q is derived in Appendix D.
A suitable Gaussian importance density for evaluating the likelihood (12) is g(α) = g(y|α)p(α)/g(y), where g(y) = g(y|α)p(α)dα and . . , C T ). Shephard and Pitt (1997) and Durbin and Koopman (1997) propose an efficient method for selecting the importance parameters B and C, which is referred to as the SPDK method. Appendix B provides the details. We apply the SPDK method in Section 5. It is simple to verify that g(α) is a multivariate normal density with inverse covariance Q * = C + Q and mean µ * = Q * −1 (B + Qµ). Shephard and Pitt (1997) and Durbin and Koopman (1997) show that we can efficiently sample from g(α) and calculate its integration constant by interpreting g(α) as an approximating linear Gaussian state space model with observations y t = C −1 t b t and the linear Gaussian measurement equation The simulation smoothing methods of de Jong and Shephard (1995) and Durbin and Koopman (2002) sample α from g(α| y) in O(T ) operations. The Kalman filter calculates the integration constant g( y) by evaluating the likelihood function for the linear state space model (14). Jungbacker and Koopman (2007) show that it is only necessary that the individual matrices C t are non-singular for the Kalman filter based procedure to be valid.

Checking the existence of moments
If the first n moments are required, then by Proposition 2 it is only necessary to check that It is possible to exploit the special structure of the non-Gaussian state space model to obtain a fast O(T ) method for checking and imposing this condition. This method is more computationally efficient than the eigenvalue method described in Section 3.1 when T is large. Suppose now that α t is scalar. The main principles and procedures in this section are best illustrated and analytically explored for the case with scalar α t . We discuss the case with multivariate α t below. The stationary univariate AR(1) model for α t is for t = 1, ..., T − 1 with |φ| < 1 and initial condition α 1 ∼ N(µ, σ 2 α ), where the unconditional variance of the states is σ 2 α = σ 2 /(1 − φ 2 ). The observation equation for the approximating linear state space model (14) is To check that this matrix is positive definite, we use Sylvester's criterion (see, e.g., Horn and Johnson, 1990, Chapter 7) which states that a symmetric real matrix is positive definite if and only if all of the leading principal minors are positive. This means that all the determinants of the square upper left sub-matrices, in our case, of the T × T matrix Q − (n − 1)C must be positive. Denote the determinant of the top left 1 × 1 sub-matrix as Λ 1 = 1 − σ 2 (n − 1)v −1 1 . Setting Λ 0 = 1, it can be checked that the top left t × t sub-matrix has determinant Λ t which satisfies the recursion for t = 2, .., T − 1 and It is therefore necessary that all the determinants Λ t , t = 1, ..., T −1 and Λ T are positive for the matrix Q − (n − 1)C to be positive definite, which gives a straightforward and fast approach for checking whether the required moments exist.

Imposing the existence of moments
We now examine a simple condition on v t which ensures that the first n moments exist. We want to find a constant value for v = v t in (16) for which the resulting matrix Q − (n − 1)C v is positive definite, with C v = diag(1/v, ..., 1/v) . This corresponds to a constant variance in the measurement density of the approximating state space form, Proposition 4 shows that Q − (n − 1)C v > 0 if v = (n − 1)σ 2 α 1+|φ| 1−|φ| when φ = 0 and v = (n − 1)σ 2 α + ǫ when φ = 0, for some small ε > 0. We set ε = 10 −5 in this paper.
Proposition 4. Consider the non-Gaussian nonlinear state space model given in (11) with scalar α t . Suppose there exists a constant scalar k and a vector δ such that Suppose that the proposal density is Proof. See Appendix A.   (18), and therefore sinusoidal paths that cross the horizontal axis. If the model is less persistent, i.e. φ is smaller, then it is unnecessary for v to be as large. For example, suppose that n = 2, with σ 2 α = 0.5 as before and φ = 0.6. Then it is only necessary that v ≥ 2 for the second moment to be finite.
Using the approximating state space model (20) with constant variance v as the importance density may be inefficient because it does not depend on the v t which we estimate from the data. To take greater account of the computed v t , let v T ). We suggest the following algorithm for finding the smallest k ≥ 0 such that Q − (n − 1)C (k) > 0. Algorithm 1. Set k = 0.
1. Check the positive definiteness of Q − (n − 1)C (k) as in Section 4.2.1.
2. If it is positive definite then stop. Otherwise, set k = k + 1 and go back to Step 1.
It is straightforward to see that the algorithm always converges, as with a large enough k, all the v ≥ v after convergence. The motivation for Algorithm 1 is to adjust the data-based variances v t as little as possible while ensuring the existence of the required moments.
After the algorithm converges, let C * = C (k) to obtain Q − (n − 1)C * > 0. In practice, it is more efficient to use a mixture of two approximating Gaussian models in the form of (16), one with ε t ∼ N(0, v * t ) and the other with ε t ∼ N(0, v t ), as justified by Proposition 3. We put a small weight π = 0.1 on the former in all the examples below.
2. If λ min > 0 then stop. Otherwise, set k = k + 1 and go back to Step 1.
After convergence, as a by-product for checking the existence condition, a resulting k > 0 means that the original C fails the condition. Computing λ min will be fast because Q − (n − 1)C (k) is symmetric and block tridiagonal (Saad, 2011). We now construct a mixture of two approximating Gaussian models (14) using the matrices C and C * = C (k) .

Panel data models with an AR(1) latent process
The random effects models considered in Section 4.1 do not take into account timevarying individual effects. A possible way to overcome this is to include in the model a time-varying individual-specific effect α it . Let y i = (y i1 , ..., y iT i ) ′ , i = 1, ..., m, be panel data, which are modeled as where the x it are p-vectors of covariates. The random effects α i = (α i1 , ..., α iT i ) ′ are assumed to follow the AR(1) process with |φ| < 1 and initial value α i1 ∼ N(µ, σ 2 /(1 − φ 2 )), u it iid ∼ N(0, 1). The parameter vector ψ consists of β, φ and σ 2 . See, e.g., Bartolucci et al. (2012).

The likelihood is
which decomposes into a product of lower dimensional integrals, where Each panel y i has the form of a non-Gaussian non-linear state space model (24) except that T i is often small, therefore the SPDK method can be used to obtain the importance parameters µ * i and Q * i −1 . Because the lengths T i of the panels are often small, we can use the general method in Section 3.1 to impose the condition for the existence of moments.

Illustrative example
Let y = (y 1 , ..., y N ) be observations from a Bernoulli distribution with success probability α ∈ (0, 1). The likelihood is with k the total number of successes. We use a normal distribution N(0.5, Q −1 ), with Q = 0.1, truncated to the interval (0, 1) as a relatively uninformative prior on α. That is, p(α) = Cκ(α)1 (0,1) (α), where κ(α) = exp − 1 2 Q(α − 0.5) 2 , with C the normalizing constant such that 1 0 p(α)dα = 1. Suppose that we are interested in estimating the posterior mean of α, i.e. we wish to estimate the integral I = 1 0 αp(α|y)dα using importance sampling. Let g(α) be the importance density. Then the unnormalized weight is given by w(α) = 1 (0,1) (α)κ(α)p(y|α)/g(α). It is possible to use the Laplace method to construct a Gaussian importance density, but the constraint on the interval (0, 1) makes it difficult to find the mode of p(α|y) ∝ p(α)p(y|α) when the mode is near either boundary. We instead use a second order Taylor expansion of l(α) = log p(y|α) at the MLE estimate α mle = k/N of α. Then, for all α ∈ (0, 1), Suppose that g(α) = N(α|α * , Q * −1 ) is the importance density. Then the second moment of the unnormalized weight E g [w(α) 2 ] is approximated as with α = (Q − Q * α * )/(2Q − Q * ) and a finite constant K. In theory, unlike the result in Proposition 1, this second moment exists even if 2Q − Q * = Q + D < 0, because the integral is taken over the interval (0, 1) rather than the real line as in (2). However, in practice, this integral can be very large if 2Q − Q * = Q + D ≪ 0. This is likely to happen as D is often very small, especially in extreme cases where the underlying α is close to 0 or 1, while we need to keep Q small enough in order to have a flat prior. We refer to this as the normal-IS case.
To demonstrate the theory presented in Section 3, we use a 2 nd -moment constrained mixture importance sampler g(α) = πg 1 (α) + (1 − π)g 2 (α) constructed as in Section 3.1. The heavy-tailed component g 1 (α) with the weight π = 0.1 is a normal density N(α * , Q −1 ) with Q modified from Q * such that 2Q − Q > 0, and g 2 (α) is the original normal density N(α * , Q * −1 ). We refer to this case as the 2 nd -IS. Using a t importance density rather than a normal importance density is sometimes recommended in the literature (Geweke, 1989). As the third importance density, we use a t density with location α * , scale Q * −1 and small degrees of freedom ν. We set ν = 5 in this example after some experimentation and refer to this as the t-IS case.
The first row of Figure 2 plots log(w(α)) against α for the three importance densities described above for two cases: a 'hard' case where N = 100, k = 7 and an 'easy' case where N = 100, k = 50. In the hard case the weight function w(α) of the normal-IS behaves badly for large α, because the posterior is highly skewed to the right and the tail of the normal importance density declines much faster than the posterior. The figure shows that t-IS and 2 nd -IS both avoid this problem. In the easy case, where the posterior is symmetric and is well approximated by a normal distribution, all the three weight functions appear well behaved. The panels in the second row plot the three importance densities, which show that imposing the condition leads to an importance density with the necessary heavy tails.
Let α s iid ∼ g(α), s = 1, ..., S = 10 6 be S samples generated from the importance density g(α). The IS estimator of the integral I and its asymptotic variance are (Geweke, 1989) We implemented the statistical test for finite variance by Koopman et al. (2009) (KSC). Define a sequence of r random variables Z 1 , . . . , Z r created from the importance weights that exceed a given threshold u as Z i = ω(α i ) − u. The threshold u is set to the 90 th percentile. Their test is based on the fact as S and u increase, the distribution of the excesses Z 1 , . . . , Z r converges to the generalized Pareto distribution The shape parameter ξ characterizes the thickness of the tails of the distribution. It can be shown that when ξ > 0, E g [ω(α) n ] = ∞ for n ≥ 1/ξ. Hence, the variance of the importance weights is infinite when ξ > 1/2. The Koopman et al. (2009) tests consists of obtaining a maximum likelihood estimate of ξ and testing the null hypothesis that ξ ≤ 1/2. We implement the Wald test version of their procedure and reject the hypothesis that the variance is finite if the returning p-value is smaller than 0.01. Table 1 reports the following items averaged over 100 replications: (i) the estimate of the posterior mean I, (ii) the ratio of the asymptotic variance of I for each IS approach divided by the asymptotic variance for the normal-IS, (iii) the proportion of the 100 replications in which the KSC test rejects the hypothesis of a finite variance, (iv) and the CPU time in seconds. In the "hard" case N = 100 and k = 7, the asymptotic variances of the t-IS and 2 nd -IS are both small compared to that of the normal-IS, and the KSC test almost always rejects the hypothesis of a finite variance of the normal-IS. The improvement of the 2 nd -IS over the normal-IS and t-IS is bigger when the ratio k/N is closer to either 0 or 1. In the "easy" case N = 100 and k = 50, the normal IS has the smallest asymptotic variance, the KSC test indicates that all the variances are finite. However, we note that the loss in efficiency when using the 2 nd -IS in the "easy" case is not as severe as the loss for the normal-IS in the "hard" case. That is, in general the 2 nd -IS should be used if we do not have much information about the target density.
We found that the KSC test is very sensitive to the selection of the threshold u.
In the examples follow we therefore do not compare our results to the test.  Table 1: Bernoulli example: The table reports the estimate of the posterior mean, the ratio of the asymptotic variances, the KSC test rejections and the CPU times for the normal-IS, the t-IS and the 2 nd -IS. All the values are averaged over 100 replications.

Likelihood evaluation for a non-Gaussian non-linear state space model
We consider a time series {y t , t = 1, ..., T } generated from the dynamic Poisson model We generate the data using the parameter values β = −1.4, φ = 0.8 and σ 2 = 0.5(1 − φ 2 ). The objective is to determine whether the importance weights from the SPDK importance sampling method for estimating the likelihood in model (24) have a finite variance. We also investigate the potential benefits of imposing a finite second moment using the 2 nd -moment constrained mixture IS as in Section 4.2.
The simulation study generates 100 time series of length T = 500 from model (24). Given the model parameter ψ = (β, φ, σ 2 ), for each time series, we perform 100 likelihood evaluations at ψ using S = 10 4 importance samples for each likelihood evaluation. Table 2 reports the following performance measures averaged over the 100 realizations of the time series {y t , t = 1, ..., T } and 100 evaluations: (i) the ratio of the variance of the weights; (ii) the ratio of the Monte Carlo standard errors (MCE) of the likelihood evaluations, with the SPDK sampler that does not impose the existence condition as the benchmark in both cases; (iii) the computing time per one likelihood evaluation; and (iv) the proportion of the 100 replications in which the original SPDK sampler satisfies the criterion (4) for the existence of the second moment.  with the SPDK sampler without imposing the existence condition as the benchmark, the CPU time per likelihood evaluation, and the proportion of the 100 replications in which the original SPDK sampler satisfies the existence condition.
We estimate the likelihood for each simulated dataset at two different values of the model parameters ψ. The first parameter vector is ψ = (−1.4, 0.99, 1), which represents an extreme case away from the DGP values. From Proposition 4, higher values of φ and σ 2 lead to a stricter restriction on the approximating linear state space model in order for the variance of the importance sampler to exist. In this case, imposing the existence condition leads to a substantial improvement in terms of both variance of the weights and the Monte Carlo errors of the likelihood evaluation. The original SPDK method always fails the criterion of Proposition 2 for the existence of the second moment. In the second case, we evaluate the likelihood at the true parameter values. In this case imposing the existence condition does not improve much on the original SPDK method. This result holds despite the theoretical infinite variance of the SPDK importance sampler.

Likelihood estimation
We investigate the potential benefit of imposing the condition that the importance weights have a finite variance in the context of likelihood evaluation at a particular value of the model parameters. As noted in Section 4.3, the likelihood of this panel data model decomposes into a product of lower dimensional integrals. We estimate each of these integrals by first using the SPDK method to obtain the importance parameters, then construct the 2 nd -moment constrained mixture IS as described in Section 3.1. For the reasons discussed in Section 5.1, we compare this 2 nd -IS to a t-IS with ν = 5 degrees of freedom, noting again that using such a t importance density does not guarantee the existence of the second moments of the weights.  We draw 100 data sets from model (25) and perform 100 likelihood evaluations at a given model parameter vector ψ. We consider two different sets of the model parameters: in the first case ψ is the vector of the true parameters that generate the data, and in the second case ψ is set to (0, 0, 0, 1) which is far from the true values. We base each likelihood evaluation on S = 10 3 importance samples. Common random numbers (CRN) are used. Similarly to Section 5.2, Table 3 reports the ratio of the variance of the weights, the ratio of the Monte Carlo standard errors (MCE) of the likelihood estimates with the t-IS as the benchmark, and computing time per one likelihood evaluation. Note that each likelihood evaluation consists of estimating m separate integrals, so the variance reported here has been averaged over the m panels. The simulation results suggest that it is beneficial to impose the condition that the importance weights have a finite variance.
Table 3 also shows that although imposing the existence of the second moments takes extra CPU time, the 2 nd -IS method overall takes less time than the t-IS because working with a t distribution takes more time than working with the constrained Gaussian mixture density.

Bayesian inference
This section explores the potential benefit of imposing the existence condition in the context of Bayesian inference. Bayesian inference for the panel data model (25) can be carried out using the MCMC approach to sample from the posterior p(ψ|y), in which the likelihood used within the Metropolis-Hastings algorithm is estimated unbiasedly by importance sampling. We put a normal prior N(0, τ 0 I p ) on β, a reference prior p(φ) ∝ 1/ 1 − φ 2 on the correlation coefficient φ (Berger and Yang, 1994) and p(σ 2 ) ∝ 1/σ 2 , and set τ 0 = 100.
We consider two MCMC samplers. The first sampler uses the 2 nd -IS to estimate the likelihood to ensure the finite second moment of the importance weights, and leads to a finite variance of the log likelihood estimator which is important for MCMC simulation based on an estimated likelihood (Pitt et al., 2012). The second MCMC sampler uses the t-IS as in Section 5.3.1 to estimate the likelihood. Each sampler is run for 50,000 iterations after 50,000 burn-in iterations using S = 200 importance samples to estimate each likelihood. The Markov chains are drawn by the adaptive random walk Metropolis-Hastings algorithm in Haario et al. (2001).
We use the integrated autocorrelation time (IACT) and the acceptance rate of the Metropolis-Hastings algorithm as the performance measures. For a scalar parameter θ the IACT is defined as (Liu, 2001)  where ρ j is the jth lag autocorrelation of the iterates of θ in the MCMC scheme. We estimate the IACT by where ρ j is the jth lag sample autocorrelation and L * = min{1000, L}, with L the first index j such that | ρ j | ≤ 2/ √ K where K is the sample size used to estimate ρ j . When summarizing the performance of the MCMC with d parameters estimated, we report, for simplicity, the average IACT over the d parameters.
We investigate the performance of the two MCMC samplers for four combinations of T and σ 2 . Figures 3 plot the iterates and the autocorrelations for the four parameters β (first two rows), φ (third) and σ 2 (last) for the case with σ 2 α = 1 and T = 80, which shows that the sampler which imposes the existence condition works better. Table 4 summarizes the acceptance rates and the IACT ratios averaged over 5 replications. The table also reports the CPU time taken for each MCMC iteration. The results suggest that the new method is beneficial in cases with large σ 2 , after taking into account the higher computing time required to impose the variance existence condition.

An application to the anti-epileptic drug dataset
The anti-epileptic drug longitudinal dataset (see, e.g., Fitzmaurice et al., 2011, p.346) consists of seizures counts on 59 epileptic patients over 5 time-intervals of treatment.
The objective is to study the effects of the anti-epileptic drug on the patients. Following Fitzmaurice et al. (2011), we consider the mixed effects Poisson regression model p(y ij |β, α i ) = Poisson(exp(η ij )), η ij = c ij + β 1 + β 2 time ij + β 3 treatment ij + β 4 time ij × treatment ij +α i1 + α i2 time ij , j = 0, 1, ..., 4, i = 1, ..., 59 and c ij is an offset. As in Fitzmaurice et al. (2011), the offset c ij = log(8) if j = 0 and c ij = log(2) for j > 0, time ij = j, treatment ij = 0 if patient i is in the placebo group and treatment ij = 1 if in the treatment group. The α i = (α i1 , α i2 ) ′ ∼ N(0, Q −1 ) are random effects that need to be integrated out. We consider a normal prior N p (0, τ 0 I p ) for β and a Wishart W(ν 0 , τ 0 I u ) prior for the precision Q. In order to have flat priors, we select τ 0 = 1000, ν 0 = u + 1. In the MCMC scheme, we first use the Leonard and Hsu transformation Q = exp(Σ), where Σ is an unconstrained symmetric matrix, and reparameterize Q by the lower-triangle elements θ Q of Σ, which is an one-to-one transformation between Q and θ Q . We then use the adaptive random walk Metropolis-Hastings algorithm in Haario et al. (2001) to sample from the posterior p(β, θ Q |y). The dimension of ψ = (β, θ Q ) is 4 + 2(2 + 1)/2 = 7. Similarly to Section 5.3.2, we run two MCMC samplers based on 2 nd -IS and t-IS. Each sampler is run for 50000 iterations after discarding 50000 burn-in iterations. The number of samples used in each likelihood estimation is S = 200. Figure 4 plots the iterates of the parameters as well as the autocorrelations generated by the two samplers. The IACT values of the samplers imposing and not imposing the existence condition are 19.6 and 29.2. That is, imposing the condition leads to a sampler that is roughly 1.5 times more efficient than not imposing. The acceptance rates averaged over replications of the sampler imposing and not imposing the existence condition are 22.3% and 18.8%, respectively. This suggests that it is beneficial to impose the finite second moment condition of the importance weight within a MCMC sampler. The two samplers respectively take 0.079 and 0.058 seconds to run each MCMC iteration. E g [ω(α) n ] can now be written as By (6), for any fixed α −j , when either α j → −∞ or α j → +∞. It follows that the integral over α j is infinite for any fixed α −j , thus E g [ω(α) n ] is infinite.
This completes the proof.