New Estimators of the Bayes Factor for Models with High-Dimensional Parameter and/or Latent Variable Spaces

Abstract

Formal Bayesian comparison of two competing models, based on the posterior odds ratio, amounts to estimation of the Bayes factor, which is equal to the ratio of respective two marginal data density values. In models with a large number of parameters and/or latent variables, they are expressed by high-dimensional integrals, which are often computationally infeasible. Therefore, other methods of evaluation of the Bayes factor are needed. In this paper, a new method of estimation of the Bayes factor is proposed. Simulation examples confirm good performance of the proposed estimators. Finally, these new estimators are used to formally compare different hybrid Multivariate Stochastic Volatility–Multivariate Generalized Autoregressive Conditional Heteroskedasticity (MSV-MGARCH) models which have a large number of latent variables. The empirical results show, among other things, that the validity of reduction of the hybrid MSV-MGARCH model to the MGARCH specification depends on the analyzed data set as well as on prior assumptions about model parameters.

1. Introduction

The Bayes factor, defined as a ratio of the marginal likelihoods for the two models being compared, is an important quantity in Bayesian model comparison and hypothesis testing, see, e.g., in [1]. The posterior odds ratio, used for comparing two competing models, is equal to the product of the prior odds and the Bayes’ factor. Therefore, one of the challenges for Bayesians is to accurately estimate the factor, especially in models with a large number of parameters and/or latent variables. Most popular methods of calculating the Bayes factor, based on estimation of the marginal likelihoods of each model separately, are very time-consuming (or even computationally infeasible) in high-dimensional cases. Therefore, methods for direct estimation of the Bayes factor, instead of estimating two marginal likelihoods, have been devised. There already exist different approaches, such as the product-space approach, reversible jump MCMC, and path sampling (see, e.g., in [2,3,4,5]), which make it possible to estimate the Bayes factor without calculation of the marginal densities of the data, but in many empirical cases (especially in the case of a large number of model parameters) they require extremely complicated and extensively time-consuming operations, with no general guarantees of success.

Note that for nested models with a large number of latent variables such as stochastic volatility models, Jacquier et al. [6] have shown how to estimate the Bayes factor directly using Markov Chain Monte Carlo (MCMC) simulations from the posterior distribution. They pointed out that the Bayes factor can be expressed as the expected value of the ratio of corresponding densities with respect to the posterior distribution of parameters and latent variables. This expected value can be estimated by averaging over the MCMC draws. Thus, using the special structure of stochastic volatility models and exploiting prior assumptions about parameters, they have estimated Bayes factors for comparing the basic stochastic volatility model to the models with fat tails and correlated errors. Unfortunately, the finite sample properties of the estimator have not been analyzed by the authors. In this paper, we will show that even in very simple models it does not perform well, but it can be improved by trimming the space of parameters and latent variables.

Now, the issue of a direct estimation of the Bayes factor will be described more formally. Let us consider a Bayesian model, in which

(1)

the space of parameters is denoted by $\Theta$ and $p\left(\theta \right)$ is a prior density function of the parameters collected in $\theta \in \Theta$,

(2)

the space of latent variables is denoted by H and $p\left(h\right|\theta )$ is a prior density function of the latent variables collected in $h\in H$, and

(3)

y is a vector of observations.

The marginal data density, $p\left(y\right)$, is defined as the integral (calculated over the whole space of parameters and latent variables, $\Theta \times H$) of the conditional data density given the vector of parameters and latent variables, $p\left(y\right|\theta ,h)$, with respect to the prior distribution of the parameters and latent variables:

$$p\left(y\right)=\underset{\Theta \times H}{\int }p(y,h,\theta )d\theta dh=\underset{\Theta \times H}{\int }p\left(y\right|h,\theta \left)p\right(h\left|\theta \right)p\left(\theta \right)d\theta dh.$$

(1)

In the majority of models it is impossible to analytically integrate out parameters and latent variables from the joint distribution for y, h, and $\theta$, $p(y,h,\theta )$. It very often results from the number of dimensions of the space of the latent variables and parameters being too large. Furthermore, a correct assessment of the marginal likelihood is computationally challenging (see, e.g., in [7,8,9], where summaries of various methods can be found). In the presence of a large number of latent variables (e.g., in stochastic volatility models) popular methods of Monte Carlo estimation of $p\left(y\right)$ for the observed vector y are computationally extremely intensive, and often they turn out to be infeasible. Fortunately, we can consider the ratio of marginal data densities instead of $p\left(y\right)$. In fact, the main criterion of comparison of two competing models $M_i$ and $M_j$ is the posterior odds ratio between the two models:

$$O_{ij}=\frac{p\left(M_i\right|y)}{p\left(M_j\right|y)}=\frac{p\left(y\right|M_i)}{p\left(y\right|M_j)}\times \frac{p\left(M_i\right)}{p\left(M_j\right)},$$

(2)

where $p\left(M_i\right|y)$ is the posterior probability of the model $M_i$, $p\left(M_i\right)$ is the prior probability of the model $M_i$, and $p\left(y\right|M_i)=\underset{{\Theta }_i\times H_i}{\int }p\left(y\right|h_i,{\theta }_i,M_i)p\left(h_i\right|{\theta }_i,M_i)p\left({\theta }_i\right|M_i)d{\theta }_{i}d{h}_i$ is the marginal data density value (the marginal likelihood) in the model $M_i.$ The ratio $B_{ij}=\frac{p\left(y\right|M_i)}{p\left(y\right|M_j)}$ is referred to as the Bayes factor in favor of the model $M_i$ against the model $M_j$, in turn, the ratio $\frac{p\left(M_i\right)}{p\left(M_j\right)}$ is the prior odds ratio. Thus, if the prior odds ratio is equal to 1 (i.e., both models are equally probable a priori, $p\left(M_i\right)=p\left(M_j\right))$, the posterior odds ratio between the two models is then equal to the Bayes factor. Moreover, if we assume that the set of models $\{M_1,\dots ,M_m\}$ is exhaustive, Bayes factors can be used to obtain the posterior model probabilities:

$$p\left(M_i\right|y)=\frac{\frac{p\left(M_i\right)p\left(y\right|M_i)}{p\left(y\right|M_k)}}{\frac{{\displaystyle \sum _{j=1}^m}p\left(M_j\right)p\left(y\right|M_j)}{p\left(y\right|M_k)}}=\frac{p\left(M_i\right)B_{ik}}{{\displaystyle \sum _{j=1}^m}p\left(M_j\right)B_{jk}},$$

(3)

for $i,k\in \{1,\cdots ,m\}$. By choosing $k=i$, we obtain

$$p\left(M_i\right|y)=\frac{p\left(M_i\right)}{{\displaystyle \sum _{j=1}^m}p\left(M_j\right)B_{ji}},$$

(4)

for $i\in \{1,\cdots ,m\}$.

Thus, the Bayes factors can be used instead of the marginal data density values. Note that in many situations it is easier to estimate the Bayes factor than the marginal likelihood. In this paper, we show how to estimate the Bayes factor in models with a large number of parameters and/or latent variables, in which calculation of the marginal data density value is numerically impossible. It will be shown that the Bayes factor is equal to the posterior mean, restricted to a certain subset D of the parameter and latent variable space, of the ratio of conditional densities of the corresponding quantities times the reciprocal of the posterior probability of the subset D. This fact leads the researcher to using arithmetic mean estimator of the ratio based on simulation from the posterior distribution restricted to the subset D.

It is well known that the Savage–Dickey density ratio and its generalizations (see in [10,11,12]) are relatively simple and widely used methods for computing Bayes factors for nested models. The Savage–Dickey method requires an estimate of the value of the posterior distribution at a single point. It is not reliable when this point lies in the tail of the posterior distribution. As has already been mentioned, the Savage–Dickey method can be used only for nested models. Our method can be applied for both nested and non-nested ones.

Note that although the posterior odds principle is fundamental, there are other Bayesian approaches to model comparison or hypothesis testing, e.g., the so-called Lindley type test (see [13]) or the Full Bayesian Significance Test (FBST), introduced by Pereira and Stern [14] as a Bayesian significance test for precise (sharp) hypotheses. A solid theoretical background of the FBST can be found in [15]. Detailed discussions and extensions of the Pereira–Stern test as well as of its applications are presented in, e.g., [16,17,18,19,20,21]. Our approach to Bayesian model comparison (or comparing hypotheses) is based on posterior probabilities associated with models (or hypotheses). This leads directly to the so-called posterior odds approach. Motivations for the formal Bayesian approach to model selection and for the use of Bayes factors are discussed in [22]. First of all, the Bayes factors and posterior model probabilities are easy to understand. Moreover, this approach is consistent, penalizes model complexity, remains conceptually the same regardless of the number of models (or hypotheses) under consideration, and does not require nested models or regular asymptotics. This approach allows one not only to test hypotheses, but also to compare them—“the process of revising prior probabilities associated with alternative hypotheses does not necessarily involve a decision to reject or accept these hypotheses” (see [1], p. 291). Furthermore, the Bayes factors make it possible to compare models whose prior distributions are different. Finally, the posterior probabilities of models or the Bayes factors can be used in the so-called Bayesian pooling approach (or Bayesian model averaging, see, e.g., in [23]). This paper is organized as follows. In Section 2, our new method is presented. Section 3 contains the simulation study. In Section 4, we present results of comparison of hybrid MSV-MGARCH models. The last section contains conclusions and direction of further research.

2. New Estimators of the Bayes Factor

It is easy to show that the Bayes factor in favor of the model $M_i$ against the model $M_j$ can be expressed as an integral:

$$\frac{p\left(y\right|M_i)}{p\left(y\right|M_j)}=\underset{{\Theta }_i\times H_i}{\int }\frac{p\left(y\right|\theta {}_i,h_i,M_i)p\left(h_i\right|\theta {}_i,M_i)p\left(\theta {}_i\right|M_i)}{p\left(y\right|\theta {}_j,h_j,M_j)p\left(h_j\right|\theta {}_j,M_j)p\left(\theta {}_j\right|M_j)}p(\theta {}_j,h_j|y,M_j)d{\theta }_{i}d{h}_i,$$

(5)

where $p\left(y\right|{\theta }_i,h_i,M_i)$, $p\left(h_i\right|{\theta }_i,M_i)$, and $p\left({\theta }_i\right|M_i)$ denote the conditional probability density function of the vector of observations (conditional sampling density), the conditional probability density function of latent variables, and the prior density of parameters in the model $M_i$, respectively. Note that when the two competing models $M_i$ and $M_j$ have the same parameter vector and the vector of latent variables (i.e., ${\theta }_i={\theta }_j=\theta$, $h_i=h_j=h$, ${\Theta }_i={\Theta }_j=\Theta$, $H_i=H_j=H)$, then the Bayes factor can be computed in a relatively straightforward way which does not require estimation of marginal data density values for each model separately. Of course, it is possible only if draws from one of the posterior distributions are available (see in [24]). We have

$$B{F}_{ij}=\underset{\Theta \times H}{\int }{r}_{ij}(\theta ,h)p(\theta ,h|y,M_j)d\theta dh,$$

(6)

where

$$r_{ij}(\theta ,h)=\frac{p\left(y\right|\theta ,h,M_i)p\left(h\right|\theta ,M_i)p\left(\theta \right|M_i)}{p\left(y\right|\theta ,h,M_j)p\left(h\right|\theta ,M_j)p\left(\theta \right|M_j)}.$$

Therefore, given the sample ${\left\{({\theta }^{\left(q\right)}{}^{\prime },h^{\left(q\right)}{}^{\prime })\right\}}_{q=1}^k$ from the posterior distribution $p(\theta ,h|y,M_j)$, an estimator of $B{F}_{ij}$ can be expressed as

$${\widehat{BF}}_{ij}=\frac{1}{k}\sum _{q=1}^k{r}_{ij}({\theta }^{\left(q\right)},h^{\left(q\right)}).$$

(7)

As was pointed out in [24], it is crucial that the variability of the ratio $r_{ij}(\theta ,h)$ is small under the posterior distribution $p(\theta ,h|y,M_j)$; otherwise, estimates of $B{F}_{ij}$ might be driven by few values of $h^{\left(q\right)}$ and ${\theta }^{\left(q\right)}$, and thus an extremely large simulation sample could be required to obtain adequate result. To deal with this problem, we propose a certain modification of the estimator (7). The idea of this modification is based on trimming the posterior sample to a certain subset of parameter and latent variable space, similarly to the idea of the correction of arithmetic mean estimator, proposed in [25]. Let us assume that $D\subseteq \Theta \times H$ and $0<Pr(D|y,M_i)<\infty$. The equality

$$\frac{p\left(y\right|M_i)}{p\left(y\right|M_j)}Pr\left(D\right|y,M_i)=\underset{\Theta \times H}{\int }{I}_D(\theta ,h)\frac{p\left(y\right|M_i)}{p\left(y\right|M_j)}p(\theta ,h|y,M_i)d\theta dh,$$

(8)

implies that

$$\frac{p\left(y\right|M_i)}{p\left(y\right|M_j)}Pr\left(D\right|y,M_i)=\underset{\Theta \times H}{\int }{I}_D(\theta ,h)\frac{p(y,\theta ,h|M_i)}{p\left(y\right|M_j)}d\theta dh.$$

We can see immediately that

$$\frac{p\left(y\right|M_i)}{p\left(y\right|M_j)}Pr\left(D\right|y,M_i)=\underset{\Theta \times H}{\int }{I}_D(\theta ,h)\frac{p\left(y\right|\theta ,h,M_i)p(\theta ,h|M_i)}{p\left(y\right|\theta ,h,M_j)p(\theta ,h|M_j)}p(\theta ,h|y,M_j)d\theta dh,$$

thus

$$\frac{p\left(y\right|M_i)}{p\left(y\right|M_j)}=\frac{1}{Pr\left(D\right|y,M_i)}\underset{\Theta \times H}{\int }{I}_D(\theta ,h)r_{ij}(\theta ,h)p(\theta ,h|y,M_j)d\theta dh.$$

(9)

Equality (9) means that the Bayes factor can be expressed as a product of the reciprocal of the posterior probability of the subset D in the model $M_i$, $Pr\left(D\right|y,M_i)$, and of the expected value of the indicator function of subset D times the ratio $r_{ij}(\theta ,h)$, $I_D(\theta ,h)r_{ij}(\theta ,h)$. This expected value is calculated with respect to the posterior distribution of the model parameters and latent variables in the model $M_j$. Therefore, given the sample ${\left\{(h^{\left(q\right)}{}^{\prime },{\theta }^{\left(q\right)}{}^{\prime })\right\}}_{q=1}^k$ from the posterior distribution $p(\theta ,h|y,M_j)$, an estimator of ${BF}_{ij}$ can be expressed as

$${\widehat{BF}}_{ij,D}=\frac{1}{\widehat{P}r\left(D\right|y,M_i)}\left[\frac{1}{k}\sum _{q=1}^k{I}_D({\theta }^{\left(q\right)},h^{\left(q\right)})r_{ij}({\theta }^{\left(q\right)},h^{\left(q\right)})\right],$$

(10)

where $\widehat{P}r\left(D\right|y,M_i)$ is an assessment of the posterior probability of the subset D in the model $M_i$. Unfortunately, this assessment requires also sampling from the posterior distribution in the model $M_i$.

Note that equality (9) can obviously be written as

$$\frac{p\left(y\right|M_i)}{p\left(y\right|M_j)}=\frac{Pr\left(D\right|y,M_j)}{Pr\left(D\right|y,M_i)}\underset{D}{\int }{r}_{ij}(\theta ,h)p(\theta ,h|D,y,M_j)d\theta dh,$$

(11)

provided that $0<Pr\left(D\right|y,M_j)<\infty$. This equality suggests that the Bayes factor can be estimated by the product of the ratio of posterior probabilities of the subset D and the sample arithmetic mean of the ratio $r_{ij}(\theta ,h)$:

$${\widehat{BF}}_{ij,D}=\frac{\widehat{P}r\left(D\right|y,M_j)}{\widehat{P}r\left(D\right|y,M_i)}\left[\frac{1}{k}\sum _{q=1}^k{r}_{ij}({\theta }^{\left(q\right)},h^{\left(q\right)})\right],$$

(12)

based on ${\left\{(h^{\left(q\right)}{}^{\prime },{\theta }^{\left(q\right)}{}^{\prime })\right\}}_{q=1}^k$ drawn from the posterior distribution given by $p(\theta ,h|y,M_j)$, truncated to the subset D, that is, $p(\theta ,h|D,y,M_j)$.

As a by-product of our considerations, we have just shown that if a Monte Carlo sampler only visits a subset of the support of the posterior distribution, the correction of the sample arithmetic mean of the ratio $r_{ij}(\theta$,$h)$ is needed.

Now, we will extend our analysis to models which contain two groups of parameters: the parameters common to both models and parameters specific to only one of them. Moreover, additional assumptions pertaining to conditional probability density functions will be discussed.

Let us assume that

(i)

${\theta }_j=\theta {}_U={({\theta }_A^{\prime },{\theta }_R^{\prime })}^{\prime }\in \Theta {}_U=\Theta {}_A\times {\Theta }_R$ denotes all the parameters of the model $M_j$ (since called $M_U)$, while ${\theta }_i={\theta }_R$ contains the parameters common to both models: $M_U$ and $M_i$ (since called $M_R)$, and the vector ${\theta }_A$ denotes specific parameters of $M_U$;

(ii)

$p\left({\theta }_U\right|M_U)=p({\theta }_A,{\theta }_R|M_U)=p\left({\theta }_A\right|M_U)p\left({\theta }_R\right|M_U)$, i.e., the random vectors ${\theta }_R$ and ${\theta }_A$ are a priori independent;

(iii)

$\underset{{\Theta }_A}{\int }p\left({\theta }_A\right|M_U)d{\theta }_A=1$, i.e., the prior distribution for the vector of parameters ${\theta }_A$ is proper; and

(iv)

$h_i=h_j=h$; $H_i=H_j=H$, i.e., competing models have the same vector of latent variables.

2.1. Different Prior Distributions of Latent Variables in Both Models

In this section, we additionally assume that

(v)

$p\left(y\right|{\theta }_U,h_U,M_U)=p\left(y\right|{\theta }_R,h_R,M_R)$, i.e., competing models have the same conditional sampling density and

(vi)

$p\left({\theta }_R\right|M_U)=p\left({\theta }_R\right|M_R)$, i.e., in both models, the common parameters have the same prior distribution.

Under above assumptions the Bayes factor in favor of the model $M_R$ versus the model $M_U$ is given by the following integral:

$$\frac{p\left(y\right|M_R)}{p\left(y\right|M_U)}=\underset{{\Theta }_A\times {\Theta }_R\times H}{\int }\frac{p\left(h\right|{\theta }_R,M_R)}{p\left(h\right|{\theta }_A,{\theta }_R,M_U)}p({\theta }_A,{\theta }_R,h|y,M_U)d{\theta }_{A}d{\theta }_{R}dh.$$

(13)

Thus, the estimator of the Bayes factor takes the form

$${\widehat{BF}}_{R,U,h}=\frac{1}{k}\sum _{q=1}^k\frac{p\left(h^{\left(q\right)}\right|{\theta }_R^{\left(q\right)},M_R)}{p\left(h^{\left(q\right)}\right|{\theta }_A^{\left(q\right)},{\theta }_R^{\left(q\right)},M_U)},$$

(14)

where ${\left\{({\theta }_A^{\left(q\right)}{}^{\prime },{\theta }_R^{\left(q\right)}{}^{\prime },h^{\left(q\right)}{}^{\prime })\right\}}_{q=1}^k$ are drawn from the posterior distribution of parameters and latent variables in the model $M_U$, i.e., from the distribution given by $p({\theta }_A,{\theta }_R,h|y,M_U)$.

On the other hand, it is easy to show that the Bayes factor in favor of the model $M_U$ against the model $M_R$ is as follows:

$$\frac{p\left(y\right|M_U)}{p\left(y\right|M_R)}=\underset{{\Theta }_A\times {\Theta }_R\times H}{\int }\frac{p\left(h\right|{\theta }_A,{\theta }_R,M_U)}{p\left(h\right|{\theta }_R,M_R)}p({\theta }_R,h|y,M_R)p\left({\theta }_A\right|M_U)d{\theta }_{A}d{\theta }_{R}dh,$$

(15)

and consequently

$${\widehat{BF}}_{U,R,h}=\frac{1}{k}\sum _{q=1}^k\frac{p\left(h^{\left(q\right)}\right|{\theta }_A^{\left(q\right)},{\theta }_R^{\left(q\right)},M_U)}{p\left(h^{\left(q\right)}\right|{\theta }_R^{\left(q\right)},M_R)},$$

(16)

where ${\left\{({\theta }_R^{\left(q\right)}{}^{\prime },h^{\left(q\right)}{}^{\prime })\right\}}_{q=1}^k$ and ${\left\{{\theta }_A^{\left(q\right)}\right\}}_{q=1}^k$ are drawn from the posterior distribution $p({\theta }_R,h|y,M_R)$ and from the prior distribution, $p\left(\theta {}_A\right|M_U)$, respectively. However, if the dimension of the vector ${({\theta }^{\prime },h^{\prime })}^{\prime }$ is high, then the estimators ${\widehat{BF}}_{U,R,h}$ and ${\widehat{BF}}_{R,U,h}$ tend to suffer from the so-called “simulation pseudo-bias”, similarly to the harmonic mean estimator (see in [26]). This “simulation pseudo-bias” of the estimators will be illustrated on the basis of simulation studies in the next section. To deal with the problem of “pseudo-bias”, we propose drawing from the posterior and prior distributions restricted to a subset of the space of parameters and latent variables with non-zero posterior probability, and correcting the arithmetic mean of the ratio by the posterior probability of the subset visited by the sampler.

Let us assume that $D_R\subseteq {\Theta }_R\times H$ and $D_A\subseteq {\Theta }_A$, $0<Pr(D_A\times D_R|M_U)<\infty$, $0<Pr(D_A\times D_R|y,M_U)<\infty$.

Starting from the identity

$$\begin{array}{cc}\hfill \frac{p\left(y\right|M_R)}{p\left(y\right|M_U)}& Pr\left(D_R\right|y,M_R)Pr\left(D_A\right|M_U)=\hfill \\ \hfill & =\underset{D_A\times D_R}{\int }\frac{p\left(y\right|M_R)}{p\left(y\right|M_U)}p({\theta }_R,h|y,M_R)p\left({\theta }_A\right|M_U)d{\theta }_{A}d{\theta }_{R}dh,\hfill \end{array}$$

(17)

we obtain

$$\frac{p\left(y\right|M_R)}{p\left(y\right|M_U)}Pr\left(D_R\right|y,M_R)Pr\left(D_A\right|M_U)=$$

$$=\underset{{\Theta }_A\times {\Theta }_R\times H}{\int }{I}_{D_A\times D_R}({\theta }_A,{\theta }_R,h)\frac{p({\theta }_R,h,y|M_R)p\left({\theta }_A\right|M_U)}{p\left(y\right|M_U)p({\theta }_A,{\theta }_R,h|y,M_U)}p({\theta }_A,{\theta }_R,h|y,M_U)d{\theta }_{A}d{\theta }_{R}dh=$$

$$=\underset{{\Theta }_A\times {\Theta }_R\times H}{\int }{I}_{D_A\times D_R}({\theta }_A,{\theta }_R,h)\frac{p({\theta }_R,h,y|M_R)p\left({\theta }_A\right|M_U)}{p({\theta }_A,{\theta }_R,h,y|M_U)}p({\theta }_A,{\theta }_R,h|y,M_U)d{\theta }_{A}d{\theta }_{R}dh.$$

On the basis of the fact that

$$\frac{p({\theta }_R,h,y|M_R)p\left({\theta }_A\right|M_U)}{p({\theta }_A,{\theta }_R,h,y|M_U)}=\frac{p\left(y\right|{\theta }_R,h,M_R)p\left(h\right|{\theta }_R,M_R)p\left({\theta }_R\right|M_R)p\left({\theta }_A\right|M_U)}{p\left(y\right|{\theta }_A,{\theta }_R,h,M_U)p\left(h\right|{\theta }_A,{\theta }_R,M_U)p({\theta }_A,{\theta }_R|M_U)},$$

and under assumptions (i)–(vi), we have

$$\frac{p({\theta }_R,h,y|M_R)p\left({\theta }_A\right|M_U)}{p({\theta }_A,{\theta }_R,h,y|M_U)}=\frac{p\left(h\right|{\theta }_R,M_R)}{p\left(h\right|{\theta }_A,{\theta }_R,M_U)}.$$

Therefore,

$$\begin{array}{cc}\hfill \frac{p\left(y\right|M_R)}{p\left(y\right|M_U)}=& \frac{1}{Pr\left(D_R\right|y,M_R)Pr\left(D_A\right|M_U)}\times \hfill \\ \hfill & \times \underset{{\Theta }_A\times {\Theta }_R\times H}{\int }{I}_{D_A\times D_R}({\theta }_A,{\theta }_R,h)\frac{p\left(h\right|{\theta }_R,M_R)}{p\left(h\right|{\theta }_A,{\theta }_R,M_U)}p({\theta }_A,{\theta }_R,h|y,M_U)d{\theta }_{A}d{\theta }_{R}dh.\hfill \end{array}$$

(18)

Identity (18) naturally leads to the following estimator of the Bayes factor:

$$\begin{array}{cc}\hfill {\widehat{BF}}_{R,U,D,h}& =\frac{1}{\widehat{P}r\left(D_R\right|y,M_R)\widehat{P}r\left(D_A\right|M_U)}\times \hfill \\ \hfill & \times \sum _{q=1}^k\frac{p\left(h^{\left(q\right)}\right|{\theta }_R^{\left(q\right)},M_R)}{p\left(h^{\left(q\right)}\right|{\theta }_A^{\left(q\right)},{\theta }_R^{\left(q\right)},M_U)}{I}_{D_A\times D_R}({\theta }_A^{\left(q\right)},{\theta }_R^{\left(q\right)},h^{\left(q\right)}),\hfill \end{array}$$

(19)

where $\{\left({\theta }_A^{\left(q\right)}{}^{\prime }\right),{\theta }_R^{\left(q\right)}{}^{\prime },h^{\left(q\right)}{}^{\prime }){\}}_{q=1}^k$ are drawn from $p({\theta }_A,{\theta }_R,h|y,M_U)$, or equivalently

$$\begin{array}{cc}\hfill {\widehat{BF}}_{R,U,D,h}& =\frac{\widehat{P}r(D_A\times D_R|y,M_U)}{\widehat{P}r\left(D_R\right|y,M_R)\widehat{P}r\left(D_A\right|M_U)}\times \hfill \\ \hfill & \times \frac{1}{k}\sum _{q=1}^k\frac{p\left(h_{D_A\times D_R}^{\left(q\right)}\right|{\theta }_{R,D_A\times D_R}^{\left(q\right)},M_R)}{p\left(h_{D_A\times D_R}^{\left(q\right)}\right|{\theta }_{A,D_A\times D_R}^{\left(q\right)},{\theta }_{R,D_A\times D_R}^{\left(q\right)},M_U)},\hfill \end{array}$$

(20)

where ${\left\{({\theta }_{A,D_A\times D_R}^{\left(q\right)}{}^{\prime },{\theta }_{R,D_A\times D_R}^{\left(q\right)}{}^{\prime },h_{D_A\times D_R}^{\left(q\right)}{}^{\prime })\right\}}_{q=1}^k$ are drawn from $p({\theta }_A,{\theta }_R,h|y,M_U)$ truncated to the subset $D_A\times D_R$. Because the posterior simulation support is the subset $D_A\times D_R$ of the parameter and latent variable space, most of the ${\left\{({\theta }_{A,D_A\times D_R}^{\left(q\right)}{}^{\prime },{\theta }_{R,D_A\times D_R}^{\left(q\right)}{}^{\prime },h_{D_A\times D_R}^{\left(q\right)}{}^{\prime })\right\}}_{q=1}^k$ have “similar” values of the ratio of the conditional densities of latent variables, so that, having probabilities $\widehat{P}r(D_A\times D_R|y,M_U)$, $\widehat{P}r\left(D_R\right|y,M_R)$, $\widehat{P}r\left(D_A\right|M_U)$, the simulation process will be numerically more efficient than in the case of unrestricted space of parameters and latent variables. The estimates will not be dominated by few values of the ratio. Therefore, a smaller simulation sample will be required to obtain adequate precision in ${\widehat{BF}}_{R,U,D,h}$.

Now suppose that $D\subseteq {\Theta }_A\times {\Theta }_R\times H$ and $0<Pr\left(D\right|M_U)<\infty$, $0<Pr\left(D\right|y,M_U)<\infty$. This time we start with the identity

$$\frac{p\left(y\right|M_U)}{p\left(y\right|M_R)}Pr\left(D\right|y,M_U)=\underset{D}{\int }\frac{p\left(y\right|M_U)}{p\left(y\right|M_R)}p({\theta }_A,{\theta }_R,h|y,M_U)d{\theta }_{A}d{\theta }_{R}dh,$$

(21)

and we obtain

$$\frac{p\left(y\right|M_U)}{p\left(y\right|M_R)}Pr\left(D\right|y,M_U)=$$

$$=\underset{{\Theta }_A\times {\Theta }_R\times H}{\int }{I}_D({\theta }_A,{\theta }_R,h)\frac{p({\theta }_A,{\theta }_R,h,y|M_U)}{p\left(y\right|M_R)p({\theta }_R,h|y,M_R)}p({\theta }_R,h|y,M_R)d{\theta }_{A}d{\theta }_{R}dh=$$

$$=\underset{{\Theta }_A\times {\Theta }_R\times H}{\int }{I}_D({\theta }_A,{\theta }_R,h)\frac{p({\theta }_A,{\theta }_R,h,y|M_U)}{p({\theta }_R,h,y|M_R)}p({\theta }_R,h|y,M_R)d{\theta }_{A}d{\theta }_{R}dh.$$

In this case,

$$\frac{p({\theta }_A,{\theta }_R,h,y|M_U)}{p({\theta }_R,h,y|M_R)}=\frac{p\left(y\right|{\theta }_A,{\theta }_R,h,M_U)p\left(h\right|{\theta }_A,{\theta }_R,M_U)p({\theta }_A,{\theta }_R|M_U)}{p\left(y\right|{\theta }_R,h,M_R)p\left(h\right|{\theta }_R,M_R)p\left({\theta }_R\right|M_R)},$$

and under assumptions (i)–(vi)

$$\frac{p({\theta }_A,{\theta }_R,h,y|M_U)}{p({\theta }_R,h,y|M_R)}=\frac{p\left(h\right|{\theta }_A,{\theta }_R,M_U)p\left({\theta }_A\right|M_U)}{p\left(h\right|{\theta }_R,M_R)}.$$

Thus,

$$\begin{array}{cc}\hfill \frac{p\left(y\right|M_U)}{p\left(y\right|M_R)}& =\frac{1}{Pr\left(D\right|y,M_U)}\times \hfill \\ \hfill & \times \underset{{\Theta }_A\times {\Theta }_R\times H}{\int }{I}_D({\theta }_A,{\theta }_R,h)\frac{p\left(h\right|{\theta }_A,{\theta }_R,M_U)}{p\left(h\right|{\theta }_R,M_R)}p({\theta }_R,h|y,M_R)p\left({\theta }_A\right|M_U)d{\theta }_{A}d{\theta }_{R}dh.\hfill \end{array}$$

(22)

Given a sample ${\left\{({\theta }_R^{\left(q\right)}{}^{\prime },h^{\left(q\right)}{}^{\prime })\right\}}_{q=1}^k$, ${\left\{{\theta }_A^{\left(q\right)}\right\}}_{q=1}^k$ from the distributions $p({\theta }_R,h|y,M_R)$ and $p\left(\theta {}_A\right|M_U)$, respectively, then, as results from (22), an estimator of the Bayes factor for the model $M_U$ against the model $M_R$ can be

$${\widehat{BF}}_{U,R,D,h}=\frac{1}{\widehat{P}r\left(D\right|y,M_U)}\frac{1}{k}\sum _{q=1}^k\frac{p\left(h^{\left(q\right)}\right|{\theta }_A^{\left(q\right)},{\theta }_R^{\left(q\right)},M_U)}{p\left(h^{\left(q\right)}\right|{\theta }_R^{\left(q\right)},M_R)}{I}_D({\theta }_A^{\left(q\right)},{\theta }_R^{\left(q\right)},h^{\left(q\right)}).$$

(23)

Additionally, if $D=D_A\times D_R\subseteq {\Theta }_A\times {\Theta }_R\times H$, then

$${\widehat{BF}}_{U,R,D,h}=\frac{\widehat{P}r\left(D_R\right|y,M_R)\widehat{P}r\left(D_A\right|M_U)}{\widehat{P}r\left(D\right|y,M_U)}\frac{1}{k}\sum _{q=1}^k\frac{p\left(h_D^{\left(q\right)}\right|{\theta }_{A,D}^{\left(q\right)},{\theta }_{R,D}^{\left(q\right)},M_U)}{p\left(h_D^{\left(q\right)}\right|{\theta }_{R,D}^{\left(q\right)},M_R)},$$

(24)

where ${\left\{({\theta }_{R,D}^{\left(q\right)}{}^{\prime },h_D^{\left(q\right)}{}^{\prime })\right\}}_{q=1}^k$, ${\left\{{\theta }_{A,D}^{\left(q\right)}\right\}}_{q=1}^k$ are drawn from the distributions determined by $p({\theta }_R,h|y,M_R)$ and $p\left(\theta {}_A\right|M_U)$, respectively, restricted to the subset D.

2.2. Different Prior Distributions for Common Parameters in Both Models

Now, let us assume that

(v)

$p\left(y\right|{\theta }_U,h_U,M_U)=p\left(y\right|{\theta }_R,h_R,M_R)$, i.e., competing models have the same sampling density and

(vii)

$p\left(h\right|{\theta }_A,{\theta }_R,M_U)=p\left(h\right|{\theta }_R,M_R)$, i.e., in both models, the latent variables have the same prior distribution, instead of (vi) $p\left({\theta }_R\right|M_U)=p\left({\theta }_R\right|M_R)$.

Then, it is easy to show that in the first case

$$\frac{p\left(y\right|M_R)}{p\left(y\right|M_U)}=\frac{1}{Pr\left(D_R\right|y,M_R)Pr\left(D_A\right|M_U)}\times$$

$$\times \underset{{\Theta }_A\times {\Theta }_R\times H}{\int }{I}_{D_A\times D_R}({\theta }_A,{\theta }_R,h)\frac{p\left({\theta }_R\right|M_R)}{p\left({\theta }_R\right|M_U)}p({\theta }_A,{\theta }_R,h|y,M_U)d{\theta }_{A}d{\theta }_{R}dh,$$

(25)

and in the second case

$$\begin{array}{cc}\hfill \frac{p\left(y\right|M_U)}{p\left(y\right|M_R)}& =\frac{1}{Pr\left(D\right|y,M_U)}\times \hfill \\ \hfill & \times \underset{{\Theta }_A\times {\Theta }_R\times H}{\int }{I}_D({\theta }_A,{\theta }_R,h)\frac{p\left({\theta }_R\right|M_U)}{p\left({\theta }_R\right|M_R)}p({\theta }_R,h|y,M_R)p\left({\theta }_A\right|M_U)d{\theta }_{A}d{\theta }_{R}dh.\hfill \end{array}$$

(26)

Basing on identities (25) and (26), the following estimators of the Bayes factors can be formulated:

$$\begin{array}{c}\hfill \widehat{B}{F}_{R,U,D,\theta }==\frac{1}{\widehat{P}r\left(D_R\right|y,M_R)\widehat{P}r\left(D_A\right|M_U)}\frac{1}{k}\sum _{q=1}^k\frac{p\left({\theta }_R^{\left(q\right)}\right|M_R)}{p\left({\theta }_R^{\left(q\right)}\right|M_U)}{I}_{D_A\times D_R}({\theta }_A^{\left(q\right)},{\theta }_R^{\left(q\right)},h^{\left(q\right)}),\end{array}$$

(27)

where ${\left\{({\theta }_A^{\left(q\right)}{}^{\prime },{\theta }_R^{\left(q\right)}{}^{\prime },h^{\left(q\right)}{}^{\prime })\right\}}_{q=1}^k$ are drawn from $p({\theta }_A,{\theta }_R,h|y,M_U)$,

$$\widehat{B}{F}_{U,R,D,\theta }=\frac{1}{\widehat{P}r\left(D\right|y,M_U)}\frac{1}{k}\sum _{q=1}^k\frac{p\left({\theta }_R^{\left(q\right)}\right|M_U)}{p\left({\theta }_R^{\left(q\right)}\right|M_R)}{I}_D({\theta }_A^{\left(q\right)},{\theta }_R^{\left(q\right)},h^{\left(q\right)}),$$

(28)

where ${\left\{({\theta }_R^{\left(q\right)}{}^{\prime },h^{\left(q\right)}{}^{\prime })\right\}}_{q=1}^k$, ${\left\{{\theta }_A^{\left(q\right)}\right\}}_{q=1}^k$ are drawn from $p({\theta }_R,h|y,M_R)$ and $p\left({\theta }_A\right|M_U)$, respectively.

2.3. Different Conditional Sampling Distributions of the Vector of Observations

Now we assume that conditions (i)–(iv), (vi), and (vii) hold. Thus, the conditional distribution of the observable vector can be different, but the prior distributions for vector of latent variables and common parameters are the same in both competing models. Then, it is easy to show that

$$\begin{array}{cc}\hfill \frac{p\left(y\right|M_R)}{p\left(y\right|M_U)}& =\frac{1}{Pr\left(D_R\right|y,M_R)Pr\left(D_A\right|M_U)}\times \hfill \\ \hfill & \times \underset{{\Theta }_A\times {\Theta }_R\times H}{\int }{I}_{D_A\times D_R}({\theta }_A,{\theta }_R,h)\frac{p\left(y\right|{\theta }_R,h,M_R)}{p\left(y\right|{\theta }_A,h,M_U)}p({\theta }_A,{\theta }_R,h|y,M_U)d{\theta }_{A}d{\theta }_{R}dh,\hfill \end{array}$$

(29)

and for the reciprocal of this Bayes factor we have

$$\begin{array}{cc}\hfill \frac{p\left(y\right|M_U)}{p\left(y\right|M_R)}& =\frac{1}{Pr\left(D\right|y,M_U)}\times \hfill \\ \hfill & \times \underset{{\Theta }_A\times {\Theta }_R\times H}{\int }{I}_D({\theta }_A,{\theta }_R,h)\frac{p\left(y\right|{\theta }_R,h,M_U)}{p\left(y\right|{\theta }_R,h,M_R)}p({\theta }_R,h|y,M_R)p\left({\theta }_A\right|M_U)d{\theta }_{A}d{\theta }_{R}dh.\hfill \end{array}$$

(30)

Similarly to the above, basing on identities (29) and (30), the following estimators of the Bayes factors can be formulated:

$$\begin{array}{cc}\hfill {\widehat{BF}}_{R,U,D,y}& =\frac{1}{\widehat{P}r\left(D_R\right|y,M_R)\widehat{P}r\left(D_A\right|M_U)}\times \hfill \\ \hfill & \times \frac{1}{k}\sum _{q=1}^k\frac{p\left(y\right|{\theta }_R^{\left(q\right)},h^{\left(q\right)},M_R)}{p\left(y\right|{\theta }_A^{\left(q\right)},h^{\left(q\right)},M_U)}{I}_{D_A\times D_R}({\theta }_A^{\left(q\right)},{\theta }_R^{\left(q\right)},h^{\left(q\right)}),\hfill \end{array}$$

(31)

where ${\left\{({\theta }_A^{\left(q\right)}{}^{\prime },{\theta }_R^{\left(q\right)}{}^{\prime },h^{\left(q\right)}{}^{\prime })\right\}}_{q=1}^k$ are drawn from $p({\theta }_A,{\theta }_R,h|y,M_U)$, and

$${\widehat{BF}}_{U,R,D,y}=\frac{1}{\widehat{P}r\left(D\right|y,M_U)}\frac{1}{k}\sum _{q=1}^k\frac{p\left(y\right|{\theta }_R^{\left(q\right)},h^{\left(q\right)},M_U)}{p\left(y\right|{\theta }_R^{\left(q\right)},h^{\left(q\right)},M_R)}{I}_D({\theta }_A^{\left(q\right)},{\theta }_R^{\left(q\right)},h^{\left(q\right)}),$$

(32)

where ${\left\{({\theta }_R^{\left(q\right)}{}^{\prime },h^{\left(q\right)}{}^{\prime })\right\}}_{q=1}^k$, ${\left\{{\theta }_A^{\left(q\right)}\right\}}_{q=1}^k$ are drawn from $p({\theta }_R,h|y,M_R)$ and $p\left({\theta }_A\right|M_U)$, respectively.

Note that the estimators (19), (22), (27), (28), (31), and (32) are based on the arithmetic mean of the ratio of densities times the indicator function of an arbitrary subset of the space of model parameters and latent variables. Additionally, this arithmetic mean is corrected by the reciprocal of the posterior probability of the subset.

3. Simulation Study

In this section, we present a simple simulation study for models with latent variables in which, after having integrated out latent variables analytically, the true values of the marginal likelihoods can easily be assessed using, e.g., the corrected arithmetic mean estimator (CAME [25]). These assessments will be applied to obtain estimates of Bayes factors used as benchmark values for evaluation of the performance of the new method proposed.

Let us consider the Student t model

$$\begin{array}{cc}\hfill & y_t|\mu ,v\sim iit(\mu ,1,v),t=1,2,...,T,\hfill \\ \hfill & \mu \sim N({\mu }_0,{\sigma }_0^2),v\sim Exp\left({\lambda }_0\right)\hfill \end{array}$$

(33)

where $t(\mu ,1,v)$ denotes the Student t distribution with mode $\mu$, precision 1, and v degrees of freedom. $N({\mu }_0,{\sigma }_0^2)$ denotes the Normal distribution with mean ${\mu }_0$ and variance ${\sigma }_0^2$, in turn, $Exp\left({\lambda }_0\right)$ stands for the Exponential distribution with mean ${\lambda }_0^{-1}$ and variance ${\lambda }_0^{-2}$. The symbol iid stands for independent and identically distributed. Thus, the random variables $y_1,\dots ,y_T$ are independent and have the same Student t distribution.

The Student t distribution can be expressed as a scale mixture of Gaussian distributions by introducing a random variable $h_t$ that is inverse-gamma distributed. The model can be written as

$$\begin{array}{cc}\hfill & y_t|{\omega }_t,\mu ,v\sim iiN(\mu ,h_t^{-1}),\hfill \\ \hfill & h_t|v\sim iiG(v/2,v/2),\hfill \\ \hfill & \mu \sim N({\mu }_0,{\sigma }_0^2),v\sim Exp\left({\lambda }_0\right),\hfill \end{array}$$

(34)

where $G(v/2,v/2)$ denotes the gamma distribution with shape $v/2$ and rate $v/2$.

To simulate datasets we generated samples of size $T=500,1000,2000$, data points from model (33) with $\mu =0,0.25,0.5,v=8$, and ${\mu }_0=0;{\sigma }_0^2=1;{\lambda }_0=0.1$.

In turn, to simulate from the posterior distribution the Gibbs algorithm was used and run for 20,000 iterations. The conditional posterior distribution for $\mu$ is Normal, i.e.,

$$\mu |h,v,y\sim N_K({\mu }_1,{\sigma }_1^2),$$

(35)

where ${\sigma }_1^2={({\displaystyle \sum _{t=1}^T}{h}_t+{\sigma }_0^{-2})}^{-1},{\mu }_1={\sigma }_1^2({\displaystyle \sum _{t=1}^T}{y}_t{h}_t+{\mu }_0{\sigma }_0^{-2})$, while the conditional posterior distribution for $h_t$ is Gamma:

$$h_t|\mu ,v,y\sim G((v+1)/2,[{(y_t-\mu )}^2+v]/2),t=1,...,T.$$

(36)

An easy computation shows that the conditional posterior distribution of v is nonstandard:

$$v|h,\mu ,y\sim p\left(v\right|h,y,\mu )\propto {\left(\frac{v}{2}\right)}^{\frac{Tv}{2}}\Gamma {\left(\frac{v}{2}\right)}^{-T}{e}^{-v\kappa /2},$$

(37)

where $\kappa ={\displaystyle \sum _{t=1}^T}{h}_t-{\displaystyle \sum _{t=1}^T}ln{h}_t+2{\lambda }_0$.

However, reliable numerical methods for generating from this distribution do exist. We use one of them, the rejection sampling proposed by [27].

In the model under consideration, ${\theta }_U={(\mu ,v)}^{\prime }$, ${\theta }_R=v$, and $h={(h_1,\dots ,h_T)}^{\prime }$; thus, the number of latent variables is equal to the number of observations: 500, 1000, and 2000, respectively. The direct computation of the marginal likelihood for model (34) is intensive and unstable due to the presence of latent variables. However, by integrating out the latent variables from the joint distribution of parameters, latent variables, and data, $p\left(y\right|h,\mu ,v\left)p\right(h,\mu ,v)$, the conditional density of the data given parameters only can be written in the following form:

$$p\left(y\right|\mu ,v)=\prod _{t=1}^T\frac{\Gamma \left(\frac{v+1}{2}\right)}{\Gamma \left(\frac{v}{2}\right)\sqrt{v\pi }}{\left(1+\frac{{(y_t-\mu )}^2}{v}\right)}^{-\frac{v+1}{2}}.$$

(38)

Unfortunately, the marginal density of the data (i.e., $p\left(y\right))$ cannot be presented in closed form. However, due to the lack of latent variables as well as thanks to the small number of parameters, the marginal data density value can easily and precisely be evaluated by using the corrected arithmetic mean estimator (CAME [25]). Estimates obtained by the use of the CAME are treated as benchmark values.

In order to show performance of our new estimators of the Bayes factor, we assume that the subset D is an intersection of the space of parameters and latent variables, obtained from the condition that the ratio of conditional density functions $p\left(y\right|{\theta }_U,h,M_U)$ and $p\left(y\right|{\theta }_R,h,M_R)$ is between the lowest and the highest values of the ratio evaluated on the basis of pseudo-random sample ${\{{\theta }^{\left(q\right)},h^{\left(q\right)}\}}_{q=1}^k$ in both models, and the hyper-cuboid limited by the range of the sampler output, i.e.,

$$D=(A\times B)\cap C,$$

where

$$A=\underset{i}{\otimes }[\underset{m\in \{M_U,M_R\}}{max}\underset{q_m}{min}\left\{{\theta }_i^{\left(q_m\right)}\right\},\underset{m\in \{M_U,M_R\}}{min}\underset{q_m}{max}\left\{{\theta }_i^{\left(q_m\right)}\right\}],$$

$$B=\underset{t}{\otimes }[\underset{m\in \{M_U,M_R\}}{max}\underset{q_m}{min}\left\{h_t^{\left(q_m\right)}\right\},\underset{m\in \{M_U,M_R\}}{min}\underset{q_m}{max}\left\{h_t^{\left(q_m\right)}\right\}],$$

$$C=\{({\theta }_U^{\prime },h^{\prime }):L\le \frac{p\left(y\right|{\theta }_U,h,M_U)}{p\left(y\right|{\theta }_R,h,M_R)}\le Q\},$$

$$L=\underset{m\in \{M_U,M_R\}}{max}\underset{q_m}{min}\left\{\frac{p\left(y\right|{\theta }_U^{\left(q_m\right)},h^{\left(q_m\right)},M_U)}{p\left(y\right|{\theta }_R^{\left(q_m\right)},h^{\left(q_m\right)},M_R)}\right\},$$

$$Q=\underset{m\in \{M_U,M_R\}}{min}\underset{q_m}{max}\left\{\frac{p\left(y\right|{\theta }_U^{\left(q_m\right)},h^{\left(q_m\right)},M_U)}{p\left(y\right|{\theta }_R^{\left(q_m\right)},h^{\left(q_m\right)},M_R)}\right\}.$$

In the restricted models ($M_R)$, it is assumed that $\mu =0$, whereas in the unrestricted model ($M_U$) $\mu \ne 0$. The ratio of density functions of the conditional distributions of the observable vector is as follows:

$$\frac{p\left(y\right|{\theta }_U,h,M_U)}{p\left(y\right|{\theta }_R,h,M_R)}=e^{-\frac{1}{2}\sum _{t=1}^T[{(y_t-\mu )}^2-y_t^2]h_t}=e^{-\frac{1}{2}{\mu }^2\sum _{t=1}^T{h}_t+\mu \sum _{t=1}^T{y}_t{h}_t}.$$

(39)

In

Table 1. Mean (M), standard deviation (SD), average error (AE; corrected arithmetic mean estimator (CAME)—estimated), and root mean square error (RMSE) in the Student t model. The number of observations and of latent variables is equal to $T=500$.

Estimator and Its Sample Characteristics	for $\mathit{\mu }=0$		for $\mathit{\mu }=1/2$
	$log{\mathit{BF}}_{\mathit{U},\mathit{R}}$	$log{\mathit{BF}}_{\mathit{R},\mathit{U}}$	$log{\mathit{BF}}_{\mathit{U},\mathit{R}}$	$log{\mathit{BF}}_{\mathit{R},\mathit{U}}$
CAME	−1.256	1.256	21.191	−21.191
$log{\widehat{BF}}_{U,R}$: M	−1.254	–	19.826	–
SD	0.012	–	0.748	–
AE	−0.001	–	−1.340	–
RMSE	0.011	–	1.631	–
$log{\widehat{BF}}_{U,R,D}$: M	−1.255	–	21.02	–
SD	0.013	–	0.357	–
AE	−0.003	–	−0.209	–
RMSE	0.013	–	0.589	–
$log{\widehat{BF}}_{R,U}$: M	–	0.498	–	−22.981
SD	–	0.165	–	0.505
AE	–	−0.776	–	−1.870
RMSE	–	0.783	–	1.935
$log{\widehat{BF}}_{R,U,D}$: M	–	1.249	–	−21.238
SD	–	0.032	–	0.284
AE	–	−0.010	–	−0.100
RMSE	–	0.036	–	0.246

Notes: Results obtained for T = 500 observations from a Student t distribution with mean μ, precision equal to 1 and 8 degrees of freedom. Estimations of the Bayes factors were repeated 100 times. The Bayes factors were estimated with Monte Carlo sampler based on 20,000 iterations. The log BF_R,U denotes the decimal logarithm of the Bayes factor in favor of model M_R against model M_U. The log BF_U,R denotes the decimal logarithm of the Bayes factor in favor of model model M_U against model M_R. ${\widehat{BF}}_{U,R}$ denotes the estimator of the Bayes factor calculated using simulation from untruncated distributions, whereas ${\widehat{BF}}_{U,R,D}$ denotes the estimator of the Bayes factor based on simulation from truncated distributions to the set D (the subscript y is omitted for convenience).

,

Table 2. Mean (M), standard deviation (SD), average error (AE; CAME—estimated), and root mean square error (RMSE) in the Student t model. The number of observations and of latent variables is equal to $T=1000$.

Estimator and Its Sample Characteristics	for $\mathit{\mu }=0$		for $\mathit{\mu }=1/2$
	$log{\mathit{BF}}_{\mathit{U},\mathit{R}}$	$log{\mathit{BF}}_{\mathit{R},\mathit{U}}$	$log{\mathit{BF}}_{\mathit{U},\mathit{R}}$	$log{\mathit{BF}}_{\mathit{R},\mathit{U}}$
CAME	−1.309	1.309	37.266	−37.266
$log{\widehat{BF}}_{U,R}$: M	−1.305	–	34.806	–
SD	0.016	–	0.738	–
AE	0.004	–	−2.485	–
RMSE	0.020	–	2.612	–
$log{\widehat{BF}}_{U,R,D}$: M	−1.306	–	37.093	–
SD	0.017	–	0.423	–
AE	0.002	–	−0.151	–
RMSE	0.018	–	0.403	–
$log{\widehat{BF}}_{R,U}$: M	–	0.415	–	−40.079
SD	–	0.201	–	0.598
AE	–	−0.813	–	−2.626
RMSE	–	0.850	–	2.753
$log{\widehat{BF}}_{R,U,D}$: M	–	1.304	–	−37.341
SD	–	0.032	–	0.258
AE	–	0.0003	–	−0.137
RMSE	–	0.028	–	0.265

Notes: Results obtained for T = 1000 observations from a Student t distribution with mean μ, precision equal to 1 and 8 degrees of freedom. Meaning of symbols used the same as in Table 1.

and

Table 3. Mean (M), standard deviation (SD), average error (AE; CAME—estimated), and root mean square error (RMSE) in the Student t model. The number of observations and of latent variables is equal to $T=2000$.

Estimator and Its Sample Characteristics	for $\mathit{\mu }=0$		for $\mathit{\mu }=1/4$
	$log{\mathit{BF}}_{\mathit{U},\mathit{R}}$	$log{\mathit{BF}}_{\mathit{R},\mathit{U}}$	$log{\mathit{BF}}_{\mathit{U},\mathit{R}}$	$log{\mathit{BF}}_{\mathit{R},\mathit{U}}$
CAME	−1.500	1.500	23.837	−23.837
$log{\widehat{BF}}_{U,R}$: M	−1.482	–	22.683	–
SD	0.017	–	0.53	–
AE	0.021	–	−1.204	–
RMSE	0.030	–	1.321	–
$log{\widehat{BF}}_{U,R,D}$: M	−1.480	–	23.638	–
SD	0.016	–	0.371	–
AE	0.024	–	−0.134	–
RMSE	0.031	–	0.325	–
$log{\widehat{BF}}_{R,U}$: M	–	0.451	–	−25.521
SD	–	0.108	–	0.517
AE	–	−1.070	–	−1.715
RMSE	–	1.072	–	1.755
$log{\widehat{BF}}_{R,U,D}$: M	–	1.497	–	−23.874
SD	–	0.039	–	0.226
AE	–	−0.010	–	−0.051
RMSE	–	0.039	–	0.283

Notes: Results obtained for T = 2000 observations from a Student t distribution with mean μ, precision equal to 1 and 8 degrees of freedom. Meaning of symbols used the same as in Table 1.

, we present results obtained by using newly proposed estimators of Bayes factors (i.e., the corrected arithmetic means of the ratio of the densities of the conditional distributions of the observable vector, ${\widehat{BF}}_{U,R,D,y}$ and ${\widehat{BF}}_{R,U,D,y})$ and uncorrected arithmetic means of the ratios ${\widehat{BF}}_{U,R,y}$ and ${\widehat{BF}}_{R,U,y}$. Table 1, Table 2 and Table 3 present means, standard deviations, root mean square errors and average errors (relative to the CAM estimates) of the decimal logarithm of estimates obtained in models under consideration. As mentioned earlier, closed-form expressions for the marginal likelihoods do not exist. Therefore, the root mean square errors and average errors are determined relative to the CAM estimates obtained for each marginal likelihood separately.

To estimate $Pr\left(D\right|y,M_U)$, $Pr\left(D_R\right|y,M_R)$, and $Pr\left(D_A\right|M_U)$, we use simulation from appropriate posterior or prior distributions, e.g.,

$$\widehat{P}r\left(D\right|y,M_U)=\frac{1}{k}\sum _{q=1}^k{I}_D({\theta }_A^{\left(q\right)},{\theta }_R^{\left(q\right)},h^{\left(q\right)}),$$

(40)

where ${\left\{({\theta }_A^{\left(q\right)}{}^{\prime },{\theta }_R^{\left(q\right)}{}^{\prime },h^{\left(q\right)}{}^{\prime })\right\}}_{q=1}^k$ are drawn from the posterior distribution $p({\theta }_A,{\theta }_R,h|y,M_U)$. The remaining probabilities are calculated in a similar manner. Furthermore, we consider uncorrected arithmetic means of the ratios ${\widehat{BF}}_{U,R,y}$ and ${\widehat{BF}}_{R,U,y}$ to investigate sampling properties of these estimators.

Figure 1 indicates that even in such simple models performance of the uncorrected estimators are not satisfactory. The estimator ${\widehat{BF}}_{R,U,y}$ is downwardly “pseudo-biased” (respective average errors are negative). On the other hand, our simulation study demonstrates that the performance of the estimators ${\widehat{BF}}_{U,R,y}$ and ${\widehat{BF}}_{R,U,y}$ can be improved by trimming the posterior simulation support to a given subset of the space of latent variables and parameters, and next making the correction of the arithmetic mean of the ratios using the posterior probabilities of the subset. As can be seen from Table 1, Table 2 and Table 3 and Figure 1, Figure 2 and Figure 3, corrected estimators of the Bayes factors perform very well in comparison to uncorrected ones. The estimators ${\widehat{BF}}_{U,R,D,y}$ and ${\widehat{BF}}_{R,U,D,y}$ produce estimates which are very close to those obtained on the basis of the CAME, while the estimators ${\widehat{BF}}_{U,R,y}$ and ${\widehat{BF}}_{R,U,y}$, based on uncorrected arithmetic mean, provide biased estimates.

Finally, note that by using both estimators for the ratio of densities and their reciprocals, i.e., ${\widehat{BF}}_{U,R,D,y}$ and ${\widehat{BF}}_{R,U,D,y}$, one can check the accuracy of assessments and of adequate (from numerical point of view) selection of the subset D. This is because, roughly speaking, the equality $B{F}_{U,R}=\frac{1}{B{F}_{R,U}}$ implies that it is natural to use ${\widehat{BF}}_{U,R,D,y}$ and $\frac{1}{{\widehat{BF}}_{R,U,D,y}}$ to estimate $B{F}_{U,R}$. Different estimates for $B{F}_{U,R}$ can indicate numerical inadequacy of selection of the subset D and/or too small simulation sample.

4. Empirical Illustration: Formal Bayesian Comparison of Hybrid MSV-MGARCH Models

In this section, the new method will be applied in order to formally compare the hybrid Multivariate Stochastic Volatility–Multivariate Generalized Autoregressive Conditional Heteroskedasticity (MSV-MGARCH) models and thus to show that it can be used in practice. The hybrid MSV-MGARCH models were proposed in [28,29,30] for modeling financial time series. These hybrid models are characterized by relatively simple model structures that exploit advantages of both model classes: flexibility of the Multivariate Stochastic Volatility (MSV) class, where volatility is modeled by latent stochastic processes, and relative simplicity of the Multivariate GARCH (MGARCH) class. The simplest specification of MSV-MGARCH model (called LN-MSF-SBEKK) is based on one latent process (Multiplicative Stochastic Factor (MSF)) and on the scalar BEKK [31] covariance structure. The LN-MSF-SBEKK structure is obtained by multiplying the SBEKK conditional covariance matrix $H_t$ by a scalar random variable $h_t$ such that $\{ln{h}_t\}$ is a Gaussian AR(1) latent process with autoregression parameter $\varphi$. The hybrid LN-MSF-SBEKK specification has been recognized in the literature [32,33,34,35] and proved to be useful in multivariate modeling of financial time series as well as of macroeconomic data [36,37,38,39,40,41]. The drawback of the LN-MSF-MGARCH process is that it cannot be treated as a direct extension of the MGARCH process with the Student t conditional distribution. When $\varphi =0$, the LN-MSF-SBEKK process reduces itself to the SBEKK process with the conditional distribution being a continuous mixture of multivariate normal distributions with covariance matrices $H_t{h}_t$ and $h_t$ log-normally distributed. However, the multivariate Student t distribution can be expressed as a scale mixture of normal distributions with the inverted gamma as a mixing distribution. This fact was exploited in [42,43], where the IG-MSF-SBEKK specification was proposed as a natural hybrid extension of the SBEKK process with the Student t conditional distribution (t-SBEKK). In the new specification, the latent process $\{ln{h}_t\}$ remains an autoregressive process of order one, but it is non-Gaussian. For $\varphi =0$ the latent variables $h_t$ (where $t\in \mathbb{Z}$) are independent and have inverted gamma (IG) distribution. Unfortunately, for $\varphi \ne 0$ the unconditional distribution of the latent variables $h_t$ is unknown. To summarize, the IG-MSF-SBEKK model is obtained by multiplying the SBEKK conditional covariance matrix $H_t$ by a scalar random variable $h_t$ coming from such latent process (with autoregression parameter $\varphi$) that, for $\varphi =0$, $h_t$ has an inverted gamma distribution. Thus, $\varphi =0$ leads to the t-SBEKK specification, in which the conditional distribution is represented as a continuous mixture of multivariate normal distributions with covariance matrices $H_t{h}_t$, where $h_t$ is inverted gamma distributed. If $\varphi \ne 0$, the latent variables $h_t$ ($t\in \mathbb{Z}$) are dependent, so (in comparison to the t-SBEKK model) in the IG-MSF-SBEKK model there exists an additional source of dependence.

Let us consider a two-dimensional autoregressive process $r_t=(r_{1,t},r_{2,t})$ defined by the equation

$$r_t={\delta }_0+r_{t-1}\Delta +{\epsilon }_t,t=1,\cdots ,T,$$

(41)

where ${\delta }_0$ and $\Delta$ are, respectively, $2\times 1$ and $2\times 2$ matrix parameters, and T is the length of the observed time series. The hybrid MSF-MGARCH specification for the disturbance term ${\epsilon }_t$ is defined by the following equality:

$${\epsilon }_t={\zeta }_t{H}_t^{1/2}{h}_t^{1/2},$$

(42)

where

$$H_t=(1-{\beta }_1-{\beta }_2)A+{\beta }_1\left({\epsilon }_{t-1}{}^{\prime }{\epsilon }_{t-1}\right)+{\beta }_2{H}_{t-1},$$

(43)

$$ln{g}_t=\varphi ln{g}_{t-1}+ln{\gamma }_t,$$

(44)

$\left\{{\zeta }_t\right\}$ is a Gaussian white noise sequence with mean vector zero and unit covariance matrix; $\left\{{\gamma }_t\right\}$ is a sequence of independent positive random variables; ${\gamma }_t$ is inverted gamma distributed with mean $\frac{v}{v-2}$ for $v>2$, i.e., $\left\{{\gamma }_t\right\}\sim iiIG(v/2,v/2)$, ${\zeta }_t\perp {\gamma }_s$ for all $t,s\in \{1,\cdots ,T\}$, $0<\left|\varphi \right|<1$; ${\beta }_1$ and ${\beta }_2$ are positive scalar parameters such that ${\beta }_1+{\beta }_2<1$; and A is a free symmetric positive definite matrix of order 2. Under (41) and (42), the conditional distribution of $r_t$ (given the past of $r_t$ and the current latent variable $h_t)$ is Normal with mean vector ${\mu }_t={\delta }_0+r_{t-1}\Delta$ and covariance matrix ${\Sigma }_t=H_t{h}_t$. For $\varphi =0$ (this value is excluded in the definition of the hybrid models under consideration) $h_t={\gamma }_t$, so the distribution of $h_t$ is an inverted gamma. In this case, $r_t$ in (41) is, given its past, an IG mixture of two-variate normal distributions $N({\mu }_t,h_t{H}_t)$, i.e., it has the two-variate Student t distribution.

The assumptions presented so far determine the conditional distribution of the observations and latent variables given the parameters. Thus, it remains to formulate the prior distributions of parameters. We use the same prior distributions as in [42,43]. Six elements of $\delta =\left({\delta }_{0}vec{(\Delta )}^{\prime }\right)$ are assumed to be a priori independent of other parameters, with the $N(0,I_6)$ prior. Matrix A has an inverted Wishart prior distribution such that $A^{-1}$ has the Wishart prior distribution with mean $I_2$; the elements of $\beta =({\beta }_1,{\beta }_2{)}^{\prime }$ are jointly uniformly distributed over the unit simplex. As regards the initial conditions for $H_t$, we take $H_0=h_0{I}_2$ and treat $h_0>0$ as an additional parameter, a priori exponentially distributed with mean 1; $\varphi$ has the uniform distribution over (−1, 1), and for v we assume the gamma distribution with mean ${\lambda }_a/{\lambda }_v$, truncated to $(2,+\infty )$. We assume that ${\lambda }_v=0.1$ with two cases: ${\lambda }_a=3$ and ${\lambda }_a=1$ (note that ${\lambda }_a=1$ leads to exponential distribution for v).

Furthermore, we use the same bivariate data sets as those modeled in [30,42,43]. The first data set consists of the official daily exchange rates of the National Bank of Poland (NBP fixing rates) for the US dollar and German mark in the period 1 February 1996–31 December 2001. The length of the modeled time series of their daily growth rates (logarithmic return rates) is 1482. The second data set consists of the daily quotations of the main index of the Warsaw Stock Exchange (WIG) and the S&P500 index of NYSE—1727 logarithmic returns are modeled from the period 8 January 1999–1 February 2006.

In order to obtain pseudo-random sample from the posterior distribution of parameters and latent variables, we use MCMC simulation techniques described in [42,44] and implemented within the GAUSS programming environment.

Using fully Bayesian methods, we want to answer the question whether the IG-MSF-SBEKK model can be reduced to the t-SBEKK case. Thus, we consider the hypothesis that a scalar parameter $\varphi =0$ (the t-SBEKK model, $M_R)$ and the alternative hypothesis that $\varphi \ne 0$ (the IG-MSF-SBEKK model, $M_U)$. For the exchange rate data, the posterior probability that $\varphi <0$ is approximately $0.001$ only and $\varphi =0$ is included in the highest posterior density (HPD) interval of probability content as high as $0.996$. Thus, Lindley’s procedure leads to the conclusion that the t-SBEKK is inadequate. But for the stock data, the posterior probability that $\varphi <0$ is $0.017$ for ${\lambda }_a=3$ and $0.054$ for ${\lambda }_a=1$, and $\varphi =0$ is included in the HPD interval of probability content $0.87$ or $0.80$, depending on the prior hyperparameter ${\lambda }_a$. In the case of the stock data, Lindley’s testing procedure yields results that the t-SBEKK model cannot be rejected by the data.

Now, our new estimators of the Bayes factors will be used. We start with the assumption that the subset D is an intersection of the subspace of parameters and latent variables, obtained by the condition that the ratio of conditional density functions $p\left(h\right|{\theta }_U,M_U)$ and $p\left(h\right|{\theta }_R,M_R)$ is between the lowest and the highest values of the ratio evaluated at pseudo-random sample ${\{{\theta }^{\left(q\right)},h^{\left(q\right)}\}}_{q=1}^k$ in both models, and the hyper-cuboid limited by the range of the sampler output, i.e., $D=(A\times B)\cap C$, where

$$A=\underset{i}{\otimes }[\underset{m\in \{M_U,M_R\}}{max}\underset{q_m}{min}\left\{{\theta }_i^{\left(q_m\right)}\right\},\underset{m\in \{M_U,M_R\}}{min}\underset{q_m}{max}\left\{{\theta }_i^{\left(q_m\right)}\right\}],$$

$$B=\underset{t}{\otimes }[\underset{m\in \{M_U,M_R\}}{max}\underset{q_m}{min}\left\{h_t^{\left(q_m\right)}\right\},\underset{m\in \{M_U,M_R\}}{min}\underset{q_m}{max}\left\{h_t^{\left(q_m\right)}\right\}],$$

$$C=\{({\theta }_U^{\prime },h^{\prime }):L\le \frac{p\left(h\right|{\theta }_U,M_U)}{p\left(h\right|{\theta }_R,M_R)}\le Q\},$$

$$L=\underset{m\in \{M_U,M_R\}}{max}\underset{q_m}{min}\left\{\frac{p\left(h^{\left(q_m\right)}\right|{\theta }_U^{\left(q_m\right)},M_U)}{p\left(h^{\left(q_m\right)}\right|{\theta }_R^{\left(q_m\right)},M_R)}\right\},$$

$$Q=\underset{m\in \{M_U,M_R\}}{min}\underset{q_m}{max}\left\{\frac{p\left(h^{\left(q_m\right)}\right|{\theta }_U^{\left(q_m\right)},M_U)}{p\left(h^{\left(q_m\right)}\right|{\theta }_R^{\left(q_m\right)},M_R)}\right\}.$$

The ratio of the densities of the conditional distributions of the vector of latent variables is as follows:

$$\frac{p\left(h\right|{\theta }_U,M_U)}{p\left(h\right|{\theta }_R,M_R)}=\prod _{t=1}^T{h}_{t-1}^{\frac{v}{2}\varphi }{e}^{\frac{v}{2h_t}(1-h_{t-1}^{\varphi })}.$$

(45)

In

Table 4. Estimates of Bayes factors for the IG-MSF-SBEKK ($M_U$) and t-SBEKK ($M_R$) models.

Estimator	Exchange Rates, $\mathit{T}=1482$		Stock Indices, $\mathit{T}=1742$
	${\mathbf{\lambda }}_{\mathbf{a}}=\mathbf{1}$	${\mathbf{\lambda }}_{\mathbf{a}}=\mathbf{3}$	${\mathbf{\lambda }}_{\mathbf{a}}=\mathbf{1}$	${\mathbf{\lambda }}_{\mathbf{a}}=\mathbf{3}$
$log{\widehat{BF}}_{U,R,D}$	1.147	1.191	−0.110	0.400
(NSE)	(0.017)	(0.024)	(0.012)	(0.014)
$log{\widehat{BF}}_{R,U,D}$	−1.134	−1.182	0.092	−0.388
(NSE)	(0.009)	(0.019)	(0.011)	(0.024)
log SD ratio	−1.159	−1.167	0.118	−0.396

Notes: The Bayes factors were estimated with Monte Carlo sampler based on 2,000,000 iterations. The log ${\widehat{BF}}_{R,U,D}$ denotes the decimal logarithm of the estimator of the Bayes factor in favor of model M_R against model M_U. NSE denotes the numerical standard error and the log SD ratio stands for the decimal logarithm of the Savage–Dickey density ratio. Scale for the strength of evidence against M_R [45]: 0 < log BF_U,R ≤ 1/2—negligible; 1/2 < log BF_U,R ≤ 1—mild; 1 < log BF_U,R ≤ 2—strong; 2 < log BF_U,R—very strong.

, we present results obtained by using newly proposed estimators of Bayes factors—the corrected arithmetic means of the ratio of the densities of the conditional distributions of latent variables, $\widehat{B}{F}_{U,R,D}$ and ${\widehat{BF}}_{R,U,D}$; the subscript h is omitted for convenience. Results obtained on the basis of our new method confirm that in the case of exchange rates, the t-SBEKK model is strongly rejected by the data, whereas it seems that the t-SBEKK specification can be adequate for stock indices. Under equal prior model probabilities, the IG-MSF-SBEKK model for exchange rates is about 14–15 times more probable a posteriori than the t-SBEKK model, indicating a strong (but not very strong) evidence against the t-SBEKK model. In turn, for stock indices the decimal logarithm of the Bayes factor of t-SBEKK relative to IG-MSF-SBEKK depends on prior distribution of the number of degrees of freedom, v. The Bayes factor in favor of t-SBEKK is equal to about $1.25$ and $0.4$ in models with ${\lambda }_a=1$ and ${\lambda }_a=3$, respectively. Of course, according to scale of Kass and Raftery (1995) these results indicate evidence for a given model that is negligible: “not worth more than a bare of mention”.

For the sake of comparison, in the last row of the Table 4 the estimates of the Savage–Dickey ratio is presented. The marginal posterior density of $\varphi$ is evaluated at $\varphi =0$ on the basis of the histogram of the sampled values of $\varphi$ from the posterior distribution. The main conclusions pertaining to the validity of the reduction of the IG-MSF-SBEKK model to the t-SBEKK one are very similar.

Additionally, it is very important to stress that our method makes it possible to compare IG-MSF-SBEKK models whose prior distributions are different (e.g., in respect of the number of degrees of freedom). As can be seen from

Table 5. Estimates of Bayes factors for the IG-MSF-SBEKK models with different prior distributions of the number of degrees of freedom.

	Exchange Rates, $\mathit{T}=1482$		Stock Indices, $\mathit{T}=1742$
Estimator	$\mathit{\varphi }\ne 0$	$\mathit{\varphi }=0$	$\mathit{\varphi }\ne 0$	$\mathit{\varphi }=0$
	(IG-MSF-SBEKK)	(t-SBEKK)	(IG-MSF-SBEKK)	(t-SBEKK)
${\widehat{BF}}_{{\lambda }_a=1,{\lambda }_a=3,D}$	8.123	9.517	0.524	1.487
$log{\widehat{BF}}_{{\lambda }_a=1,{\lambda }_a=3,D}$	0.910	0.979	−0.281	0.172
${\widehat{BF}}_{{\lambda }_a=3,{\lambda }_a=1,D}$	0.123	0.105	1.850	0.622
$log{\widehat{BF}}_{{\lambda }_a=3,{\lambda }_a=1,D}$	−0.909	−0.978	0.269	−0.206

Notes: The Bayes factors were estimated with Monte Carlo sampler based on 2,000,000 iterations. The log ${\widehat{BF}}_{{\lambda }_a=1,{\lambda }_a=3,D}$ denotes the decimal logarithm of the estimator of the Bayes factor in favor of the IG-MSF-SBEKK model with ${\lambda }_a=1$ against the IG-MSF-SBEKK model with ${\lambda }_a=3$.

, for exchange rates the IG-MSF-SBEKK and t-SBEKK models with exponential distribution for v are more probable a posteriori than those with the gamma distribution for v with the hyperparameter ${\lambda }_a=3$. In contrary, for stock indices the IG-MSF-SBEKK model with the hyperparameter ${\lambda }_a=3$ is more probable a posteriori than the same model with ${\lambda }_a=1$, but this improvement is negligible.

5. Discussion

In this paper, a new method of the estimation of the Bayes factor is proposed. The idea of proposed estimators is based on correction of the arithmetic mean estimator of the ratio of conditional distributions by trimming the posterior sample to a certain subset of the space of parameters and latent variables, D, and correcting the arithmetic mean by the posterior probabilities of this subset. The new method makes it possible to compare a finite number of different models with a large number of parameters and/or latent variables in respect to the goodness of fit measured by their posterior probabilities. Our simulation study and empirical illustration show that the method performs well.

In this paper, the question of an adequate selection of the subset D used in the method proposed is not addressed. The choice of the subset which minimizes the variance of the estimator remains an open question. The simulation study and empirical example considered in the paper indicate that the choice of D as an intersection of the parameter space (with additional restrictions) and the hyper-cuboid limited by the range of the posterior sampler outputs leads to acceptable results.