Intrinsic Bayesian estimation of linear time series models

Abstract

Intrinsic loss functions (such as the Kullback–Leibler divergence, i.e. the entropy loss) have been used extensively in place of conventional loss functions for independent samples. But applications in serially correlated samples are scant. In the present study, we examine Bayes estimator of Linear Time Series (LTS) model under the entropy loss. We derive the Bayes estimator and show that it involves a frequentist expectation of regressors. We propose a Markov Chain Monte Carlo procedure that jointly simulates the posteriors of the LTS parameters with frequentist expectation of regressors. We conduct Bayesian estimation of an LTS model for seasonal effects in some U.S. macroeconomic variables.

Keywords:

• entropy loss
• latent parameters
• reference prior

AMS 1991 Subject Classifications:

• Primary 62F15
• Secondary 62A15
• 62C10
• 62F25

1. Introduction

To analyse dynamics of multi-variate economic systems, researchers frequently employ Linear Time Series (LTS) models (see for example, Sims, 1980 and the ensuing literature). Bayesian inference of such models often requires point estimate of parameters because reporting the entire posterior distribution is made difficult by a prohibitively large number of parameters. A critical aspect in Bayesian estimation is the choice of loss function.

In this study, we derive Bayes estimator of LTS models based on the intrinsic loss. We illustrate a computational problem arising from serial correlation in the models when applying the intrinsic loss and our solution to the problem.

A loss function $L(\theta ,\hat{\theta })$ measures the distance between the parameter θ and its estimate $\hat{\theta }$. Such a metric is often specified for convenience given the problem at hand instead of grounding on a general principle. Bernardo and Juárez (2003) noted that for inferential purposes, what matters most is not the distance between θ and $\hat{\theta }$, instead it is the intrinsic loss – the distance between the probability model $f(x\mid \hat{\theta })$ (corresponding to the estimate $\hat{\theta }$) and $f(x\mid \theta )$ (corresponding to the actual parameter θ). Robert (1994, 1996) proposed using the logarithmic divergence (also known as the Kullback–Leibler divergence or the entropy loss) as the intrinsic loss. The intrinsic loss has a number of desirable properties not generally possessed by conventional loss functions. For example, it is invariant to transformation of data x or parameter θ and has the additive property that the loss of the sum of two independent data sets is the sum over the two losses corresponding to each data set.

The intrinsic loss has been used for independent samples for Bayesian estimation. It has also be used in various contexts for time series data. For instance, Kitamura and Stutzer (1997) used the Kullback–Leibler distance to derive a frequentist estimator for nonlinear models. Solo et al. (2001) used the Kullback–Leibler distance for evaluation of signal processing model. Robertson et al. (2005) used the entropy divergence for evaluation of forecasting density. Fernandez-Villaverde and Rubio-Ramirez (2004) used the Kullback–Leibler distance to evaluate dynamic equilibrium models in economics. However, employing the intrinsic loss for Bayesian estimation of time series models leads to technical challenges.

To illustrate the difference of the intrinsic loss in independent models and serially correlated models, consider the following examples. First, suppose $\mathit{y}=\{y_1,\dots ,y_T\}$, where $y_t$ ($t=1,\dots ,T$) are independently identically distributed (iid) $N(\rho ,1)$ and we are interested in estimating the mean parameter ρ under the entropy loss $\begin{array}{rl}\kappa (\hat{\rho }\mid \rho )& =\int \log \left\{\frac{f(\mathit{y}\mid \rho )}{f(\mathit{y}\mid \hat{\rho })}\right\}f(\mathit{y}\mid \rho )\mathrm{d}\mathit{y}\\ & ={\text{E}}_{\mathit{y}\mid \rho }\log \left\{\frac{f(\mathit{y}\mid \rho )}{f(\mathit{y}\mid \hat{\rho })}\right\}.\end{array}$ By the assumption on the model $f(\mathit{y}\mid \rho )\propto \exp \{-\frac{1}{2}\times \sum _{t=1}^T(y_t-\rho {)}^2\}$. It is easy to verify that $\kappa (\hat{\rho }\mid \rho )=\left(T/2\right)(\rho -\hat{\rho }{)}^2$. In this case, the intrinsic loss coincides with the commonly used quadratic loss, which implies that the Bayes estimator of ρ is the posterior mean. Now consider an AR(1) model: $y_t=\rho y_{t-1}+{ϵ}_t,$ for $t=1,\dots ,T,$ where ${ϵ}_t$ is iid $N(0,1)$, and ρ is the only unknown parameter.

The entropy loss is still $\begin{array}{rl}\kappa (\hat{\rho }\mid \rho )& =\int \log \left\{\frac{f(\mathit{y}\mid \rho )}{f(\mathit{y}\mid \hat{\rho })}\right\}f(\mathit{y}\mid \rho )\mathrm{d}\mathit{y}\\ & ={\text{E}}_{\mathit{y}\mid \rho }\log \left\{\frac{f(\mathit{y}\mid \rho )}{f(\mathit{y}\mid \hat{\rho })}\right\},\end{array}$ but now f is the density of the AR variable. Substitute in the distribution of data $\begin{array}{r}\kappa (\hat{\rho }\mid \rho )=\frac{(\rho -\hat{\rho }{)}^2}{2}{\text{E}}_{\mathit{y}\mid \rho }\sum _{t=1}^T{y}_{t-1}^2=(\rho -\hat{\rho }{)}^2\delta (y_0,\rho ),\end{array}$ where $\begin{array}{rl}\delta (y_0,\rho )& =\frac{1}{2}{\text{E}}_{\mathit{y}\mid \rho }\sum _{t=1}^T{y}_{t-1}^2\\ & =\frac{1}{2}\sum _{t=1}^T\left\{{\rho }^{2(t-1)}{y}_0^2+1+{\rho }^2+\cdots +{\rho }^{2(t-2)}\right\}\\ & =\frac{1}{2}\left\{\frac{1-{\rho }^{2T}}{1-{\rho }^2}{y}_0^2+\frac{1}{1-{\rho }^2}\left(T-\frac{1-{\rho }^{2T}}{1-{\rho }^2}\right)\right\}.\end{array}$ It is obvious that $\delta (y_0,\rho )$ is an increasing function of ${\rho }^2$ and is nonnegative for any ρ. A Bayes estimator (which is called a generalised Bayes estimator if the prior is improper) minimises the Bayesian posterior expected loss. If the entropy loss is employed, the Bayes estimator for ρ with a given initial condition $y_0$ is ${\hat{\rho }}_E=\text{argmin}{\text{E}}_{\rho \mid \mathit{y}}\left[\delta \right(y_0,\rho \left)\right(\rho -\hat{\rho }{)}^2]={\text{E}}_{\rho \mid \mathit{y}}\{\delta (y_0,\rho )\rho \left\}/{\text{E}}_{\rho \mid \mathit{y}}\right\{\delta (y_0,\rho )\}.$ Note that if ρ is positive, ρ and $\delta (y_0,\rho )$ are positively correlated. It follows that the Bayes estimator under the entropy loss for a positive ρ is larger than the posterior mean. It is well known that the MLE (and the posterior mean under constant prior) of ${\hat{\rho }}_M$ is biased downward, especially when the true parameter is close to unity (see MacKinnon & Smith, 1998) Note that under the constant prior $\sum _{t=1}^T{y}_{t-1}^2$ is the posterior precision (i.e. the inverse of posterior variance) for ρ. Hence the weight for the square of estimation error in the intrinsic loss function, $\delta (y_0,\rho )={\text{E}}_{\mathit{y}\mid \rho }\sum _{t=1}^T{y}_{t-1}^2$, is larger in the region of ρ where the posterior precision is high. It is in the spirit of Zellner's (1978, 1998) ‘precision of the estimation’ loss. This is in contrast of quadratic loss that imposes the same weight on all regions of ρ.

Now we turn to the model of interest. The LTS of a p dimensional column endogenous variable ${\mathit{y}}_t$ and a q dimensional column exogenous (predetermined) variable ${\mathit{x}}_{0t}$ $(t=1,\dots ,T)$ has the form: (1) $\begin{array}{r}{\mathit{y}}_t^{\prime }={\mathit{x}}_{0t}^{\prime }{\mathit{B}}_0+\sum _{j=1}^L{\mathit{y}}_{t-j}^{\prime }{\mathit{B}}_j+{\mathit{ϵ}}_t,\end{array}$(1) where L is a known positive integer, ${\mathit{B}}_0$ is a $q\times p$ unknown matrix, ${\mathit{B}}_j$ is an unknown $p\times p$ matrix, ${\mathit{ϵ}}_1,\dots ,{\mathit{ϵ}}_T$ are iid $N_p(\mathbf{0},\mathbf{\Sigma })$ errors, and $\mathbf{\Sigma }$ is an unknown $p\times p$ positive definite matrix.

A special case of the above LTS model is all of the lag coefficients ${\mathit{B}}_1,\dots ,{\mathit{B}}_L$ are zero (i.e. all regressors (with q>p) are exogenous variables.) The exogenous variables may be functions of time. For example, in modelling of climate temperature or holiday consumer spending, seasonal dummies may be introduced in the model. In economic applications, these exogenous variables may also be variables of government policies. Another special case is ${\mathit{x}}_{0t}$ is a $1\times p$ constant vector with elements of unity. The regressors only include lags of the variable ${\mathit{y}}_t$. The LTS model becomes a Vector AutoRegression (VAR), which is commonly used for modelling of macroeconomic time series.

We can rewrite Equation (1(1) $\begin{array}{r}{\mathit{y}}_t^{\prime }={\mathit{x}}_{0t}^{\prime }{\mathit{B}}_0+\sum _{j=1}^L{\mathit{y}}_{t-j}^{\prime }{\mathit{B}}_j+{\mathit{ϵ}}_t,\end{array}$(1) ) in the familiar matrix form (2) $\begin{array}{r}\mathit{Y}=\mathit{X}\mathbf{\Phi }+\mathit{ϵ},\end{array}$(2) where $\begin{array}{rl}\mathit{X}& =\left(\begin{array}{c}{\mathit{x}}_1^{\prime }\\ ⋮\\ {\mathit{x}}_T^{\prime }\end{array}\right)=({\mathit{X}}_0,{\mathit{X}}_1);{\mathit{X}}_0=\left(\begin{array}{c}{\mathit{x}}_{01}^{\prime }\\ ⋮\\ {\mathit{x}}_{0T}^{\prime }\end{array}\right),\\ {\mathit{X}}_1& =\left(\begin{array}{ccc}{\mathit{y}}_0^{\prime }& \cdots & {\mathit{y}}_{1-L}^{\prime }\\ ⋮& ⋮& ⋮\\ {\mathit{y}}_{T-1}^{\prime }& \cdots & {\mathit{y}}_{T-L}^{\prime }\end{array}\right);\\ \mathbf{\Phi }& =\left(\begin{array}{c}{\mathit{B}}_0\\ {\mathbf{\Phi }}_1\end{array}\right)=\left(\begin{array}{c}{\mathit{B}}_0\\ {\mathit{B}}_1\\ ⋮\\ {\mathit{B}}_L\end{array}\right);\mathit{Y}=\left(\begin{array}{c}{\mathit{y}}_1^{\prime }\\ ⋮\\ {\mathit{y}}_T^{\prime }\end{array}\right),\mathit{ϵ}=\left(\begin{array}{c}{\mathit{ϵ}}_1^{\prime }\\ ⋮\\ {\mathit{ϵ}}_T^{\prime }\end{array}\right).\end{array}$ Here ${\mathit{X}}_0$ and ${\mathit{X}}_1$ are $T\times q$ and $T\times Lp$; the former does not depend on parameters $\mathbf{\Sigma }$ and $\mathbf{\Phi }$, but the latter does. $\mathit{Y}$ and $\mathit{ϵ}$ are $T\times p$ matrices, $\mathbf{\Phi }$ is a $(q+Lp)\times p$ matrix of unknown parameters, ${\mathit{x}}_t^{\prime }$ is a $1\times (q+Lp)$ row vector, and $\mathit{X}$ is a $T\times (q+Lp)$ matrix of observations. The likelihood function of $(\mathbf{\Phi },\mathbf{\Sigma })$ based on $\mathit{Y}$ is then (3) $\begin{array}{rl}& f(\mathit{Y}\mid \mathbf{\Phi },\mathbf{\Sigma })\\ & \propto \frac{1}{|\mathbf{\Sigma }{|}^{T/2}}\exp \left\{-\frac{1}{2}\sum _{t=1}^T({\mathit{y}}_t-{\mathit{x}}_t\mathbf{\Phi }){\mathbf{\Sigma }}^{-1}({\mathit{y}}_t-{\mathit{x}}_t\mathbf{\Phi }{)}^{\prime }\right\}\\ & =\frac{1}{|\mathbf{\Sigma }{|}^{T/2}}\mathrm{e}\mathrm{t}\mathrm{r}\left\{-\frac{1}{2}(\mathit{Y}-\mathit{X}\mathbf{\Phi }){\mathbf{\Sigma }}^{-1}(\mathit{Y}-\mathit{X}\mathbf{\Phi }{)}^{\prime }\right\}.\end{array}$(3) Here and hereafter $\mathrm{e}\mathrm{t}\mathrm{r}\left(\mathbf{A}\mathbf{\right)}$ is $\exp \left(\text{tr}\right(\mathbf{A}\mathbf{\left)}\mathbf{\right)}$ of a matrix $\mathbf{A}\mathbf{.}$

The present paper achieves two goals. The first one is derivation of the Bayes estimator of LTS model under the entropy loss. We show that the entropy loss on $(\mathbf{\Phi },\mathbf{\Sigma })$ is non-separable in $\mathbf{\Sigma }$ and $\mathbf{\Phi }$, which can be written as the sum of losses pertaining to the covariance matrix $\mathbf{\Sigma }$ and normalised estimation error of $\mathbf{\Phi }$. The form of the $\mathbf{\Phi }$-part of the entropy loss for LTS is $\mathrm{t}\mathrm{r}\left\{{\hat{\mathbf{\Sigma }}}^{-1}\right(\mathbf{\Phi }-\hat{\mathbf{\Phi }}{)}^{\prime }{\mathrm{E}}_{(\mathit{X}\mid \mathbf{\Phi },\mathbf{\Sigma })}\left({\mathit{X}}^{\prime }\mathit{X}\right)(\mathbf{\Phi }-\hat{\mathbf{\Phi }})\}$, where $\hat{\mathbf{\Sigma }}$ is the Bayes estimator of $\mathbf{\Sigma }$. Under the entropy loss, the Bayes estimator distinctly differs from the posterior mean and differs from that of the iid multivariate normal model. The part of the intrinsic loss function associated with the regression coefficients turns out to be related with a conventional loss function. For estimation of a matrix parameter such as $\mathbf{\Phi }$ in the simultaneous equations context, Zellner (1978, 1998) proposed a ‘precision of estimation’ loss that can also be written as $\mathrm{t}\mathrm{r}\left\{{\mathbf{\Sigma }}^{-1}\right(\mathbf{\Phi }-\hat{\mathbf{\Phi }}{)}^{\prime }\left({\mathit{X}}^{\prime }\mathit{X}\right)(\mathbf{\Phi }-\hat{\mathbf{\Phi }})\}$. However, in Zellner's simultaneous equations model, ${\mathit{X}}^{\prime }\mathit{X}$ is taken as given, but in LTS the predetermined variable ${\mathit{X}}^{\prime }\mathit{X}$ depends on parameters $(\mathbf{\Phi },\mathbf{\Sigma })$.

The second goal concerns numerical estimation of the intrinsic Bayes estimator via Markov Chain Monte Carlo (MCMC). We propose a general algorithm that generate regressors as latent parameters in simulation of posteriors of parameters of LTS models. Data augmentation in this study differs from that in Tanner and Wang's 1987 seminal paper in motivation and implementation. Tanner and Wang use data augmentation to alter the likelihood function for easier MCMC simulation of the posteriors. In this study, the likelihood function of the generated data is the same as the likelihood of the sample data. Here, data augmentation does not make it easier for posterior simulation. Instead, it makes it possible to compute frequentist moment ${\mathrm{E}}_{(\mathit{X}\mid \mathbf{\Phi },\mathbf{\Sigma })}\left({\mathit{X}}^{\prime }\mathit{X}\right)$ of the LTS variables. The frequentist moment, simulated jointly with parameters, is used to produce Bayes estimates under the entropy loss.

Besides the choice of loss function, the choice of prior also plays a pivotal role in Bayesian estimation. Jeffreys prior on $\mathbf{\Sigma }$ (see Zellner, 1971) is a noninformative prior for $\mathbf{\Sigma }$ that gives rise to conditional posteriors in well known distributions. Ni et al. (2007) conducted Bayesian estimation of VAR model under the entropy loss, using the Jefferys prior for $\mathbf{\Sigma }$. However, despite its popularity the Jeffreys prior is known for producing unsatisfactory results in multi-parameter settings. In this study we simulate the LTS model under a combination of normal prior on regression parameters and Yang and Berger (1994) reference prior on $\mathbf{\Sigma }$. The conditional posteriors of $\mathbf{\Sigma }$ are simulated using a Metropolis-Hastings algorithm. Our empirical application shows that despite the fact that LTS models involve a large number of parameters and a large number of latent variables, the data-augmentation algorithm is quite efficient.

In Section 2 of the paper, we derive the Bayes estimator of LTS models under the entropy loss function and discuss computation of the weighting matrix in the Bayes estimator. In Section 3, we present a general algorithm using generated data as latent parameters. In Section 4, we lay out the MCMC algorithm for computing ($\mathbf{\Phi }$, $\mathbf{\Sigma }$) in the LTS model. In Section 5, we first compare intrinsic Bayes estimator with other estimators in a numerical example and then estimate a LTS model using seasonally unadjusted macroeconomic data. In Section 6 we offer concluding remarks.

2. Entropy loss function for the iid multivariate and LTS models

2.1. Entropy loss function for the iid model

We first consider the entropy loss function (Robert, 1994, p. 74) for a multivariate normal distribution. Let $\mathit{Y}=({\mathit{x}}_1,{\mathit{x}}_2,\dots ,{\mathit{x}}_T)$ be a random sample from $N_p(\mathit{\mu },\mathbf{\Sigma })$. One can compute the entropy loss function as $\begin{array}{rl}& L(\tilde{\mathit{\mu }},\tilde{\mathbf{\Sigma }};\mathit{\mu },\mathbf{\Sigma })\\ & =\int \log \left\{\frac{p(\mathit{x}\mid \mathit{\mu },\mathbf{\Sigma })}{p(\mathit{x}\mid \tilde{\mathit{\mu }},\tilde{\mathbf{\Sigma }})}\right\}p(\mathit{x}\mid \mathit{\mu },\mathbf{\Sigma })\text{d}\mathit{x}\\ & =\frac{T}{2}\left\{\mathrm{t}\mathrm{r}\left(\mathbf{\Sigma }{\tilde{\mathbf{\Sigma }}}^{-1}\right)-\log |\mathbf{\Sigma }{\tilde{\mathbf{\Sigma }}}^{-1}|-p\right.\\ & +\left.(\tilde{\mathit{\mu }}-\mathit{\mu }{)}^{\prime }{\tilde{\mathbf{\Sigma }}}^{-1}(\tilde{\mathit{\mu }}-\mathit{\mu })\right\}.\end{array}$ Here $p(\mathit{x}\mid \mathit{\mu },\mathbf{\Sigma })$ is the density of $N_p(\mathit{\mu },\mathbf{\Sigma })$. Clearly, the loss function has two parts. One part is related to the means $\tilde{\mathit{\mu }}$ and $\mathit{\mu }$ (with $\tilde{\mathbf{\Sigma }}$ as the weighting matrix), and the other part is related to $\tilde{\mathbf{\Sigma }}$ and $\mathbf{\Sigma }.$ The following fact states that Bayes estimator for μ is the posterior mean of μ but that of $\mathbf{\Sigma }$ is larger than the posterior mean of $\mathbf{\Sigma }$.

Fact 2.1

Under the entropy loss L, the generalised Bayes estimator of $(\mathit{\mu },\mathbf{\Sigma })$ is $\begin{array}{rl}{\hat{\mathit{\mu }}}_{\mathrm{i}\mathrm{i}\mathrm{d}}& =\mathrm{E}\left\{\mathit{\mu }|\mathit{Y}\right\},\\ {\hat{\mathbf{\Sigma }}}_{\mathrm{i}\mathrm{i}\mathrm{d}}& =\mathrm{E}(\mathbf{\Sigma }\mid \mathit{Y})+\mathrm{v}\mathrm{a}\mathrm{r}(\mathit{\mu }\mid \mathit{Y}).\end{array}$

Note that $\mathit{Y}$ represents data, expectation $\mathrm{E}(.)$ and variance $\mathrm{v}\mathrm{a}\mathrm{r}(.)$ are with respect to the posterior distribution. The proof is in the appendix.

2.2. Entropy loss functions for LTS models

Recall that for the LTS model (2(2) $\begin{array}{r}\mathit{Y}=\mathit{X}\mathbf{\Phi }+\mathit{ϵ},\end{array}$(2) ), the likelihood function of $(\mathbf{\Phi },\mathbf{\Sigma })$ is of the form (3(3) $\begin{array}{rl}& f(\mathit{Y}\mid \mathbf{\Phi },\mathbf{\Sigma })\\ & \propto \frac{1}{|\mathbf{\Sigma }{|}^{T/2}}\exp \left\{-\frac{1}{2}\sum _{t=1}^T({\mathit{y}}_t-{\mathit{x}}_t\mathbf{\Phi }){\mathbf{\Sigma }}^{-1}({\mathit{y}}_t-{\mathit{x}}_t\mathbf{\Phi }{)}^{\prime }\right\}\\ & =\frac{1}{|\mathbf{\Sigma }{|}^{T/2}}\mathrm{e}\mathrm{t}\mathrm{r}\left\{-\frac{1}{2}(\mathit{Y}-\mathit{X}\mathbf{\Phi }){\mathbf{\Sigma }}^{-1}(\mathit{Y}-\mathit{X}\mathbf{\Phi }{)}^{\prime }\right\}.\end{array}$(3) ). The entropy loss for the LTS model is (4) $\begin{array}{r}L(\tilde{\mathbf{\Phi }},\tilde{\mathbf{\Sigma }};\mathbf{\Phi },\mathbf{\Sigma })={\mathrm{E}}_{(\mathit{Y}\mid \mathbf{\Phi },\mathbf{\Sigma })}\log \left\{\frac{f(\mathit{Y}\mid \mathbf{\Phi },\mathbf{\Sigma })}{f(\mathit{Y}\mid \tilde{\mathbf{\Phi }},\tilde{\mathbf{\Sigma }})}\right\},\end{array}$(4) where for computing the expectation on the right-hand side, $(\tilde{\mathbf{\Phi }},\tilde{\mathbf{\Sigma }})$ is not a function of $\mathit{Y}$. The entropy loss $L(\tilde{\mathbf{\Phi }},\tilde{\mathbf{\Sigma }};\mathbf{\Phi },\mathbf{\Sigma })$ can be decomposed into two parts. One part measures the loss associated with the covariance matrix only, while the second part measures the loss of coefficients $\mathbf{\Phi }$, but related to the covariance matrix $\mathbf{\Sigma }$ as well. Because ${\mathit{y}}_t-{\mathit{X}}_t\mathbf{\Phi }$ $(t=1,\dots ,T)$ are iid $N_p(\mathbf{0},\mathbf{\Sigma })$, we have $\begin{array}{rl}& {\mathrm{E}}_{(\mathit{X}\mid \mathbf{\Phi },\mathbf{\Sigma })}(\mathit{Y}-\mathit{X}\mathbf{\Phi })=\mathbf{0},\\ & {\mathrm{E}}_{(\mathit{X}\mid \mathbf{\Phi },\mathbf{\Sigma })}\left\{\right(\mathit{Y}-\mathit{X}\mathbf{\Phi }{)}^{\prime }\mathit{X}\mathbf{\Phi }\}=\mathbf{0},\\ & {\mathrm{E}}_{(\mathit{X}\mid \mathbf{\Phi },\mathbf{\Sigma })}\left\{\right(\mathit{Y}-\mathit{X}\mathbf{\Phi }{)}^{\prime }(\mathit{Y}-\mathit{X}\mathbf{\Phi })\}=T\mathbf{\Sigma }.\end{array}$ Then $\begin{array}{rl}& {\mathrm{E}}_{(\mathit{X}\mid \mathbf{\Phi },\mathbf{\Sigma })}\{\log \frac{|\mathbf{\Sigma }{|}^{-T/2}\mathrm{e}\mathrm{t}\mathrm{r}\left\{-\frac{1}{2}\right(\mathit{Y}-\mathit{X}\mathbf{\Phi }\left){\mathbf{\Sigma }}^{-1}\right(\mathit{Y}-\mathit{X}\mathbf{\Phi }{)}^{\prime }\}}{|\hat{\mathbf{\Sigma }}{|}^{-T/2}\mathrm{e}\mathrm{t}\mathrm{r}\left\{-\frac{1}{2}\right(\mathit{Y}-\mathit{X}\hat{\mathbf{\Phi }}\left){\hat{\mathbf{\Sigma }}}^{-1}\right(\mathit{Y}-\mathit{X}\hat{\mathbf{\Phi }}{)}^{\prime }\}}\}\\ & =\frac{T}{2}(\log |\hat{\mathbf{\Sigma }}{\mathbf{\Sigma }}^{-1}|-p)\\ & +\frac{1}{2}\mathrm{t}\mathrm{r}\left[{\mathrm{E}}_{(\mathit{X}\mid \mathbf{\Phi },\mathbf{\Sigma })}\left\{(\mathit{Y}-\mathit{X}\hat{\mathbf{\Phi }}){\hat{\mathbf{\Sigma }}}^{-1}(\mathit{Y}-\mathit{X}\hat{\mathbf{\Phi }}{)}^{\prime }\right\}\right]\\ & =\frac{T}{2}(\log |\hat{\mathbf{\Sigma }}{\mathbf{\Sigma }}^{-1}|-p)+\frac{T}{2}\mathrm{t}\mathrm{r}\left(\hat{\mathbf{\Sigma }}{\mathbf{\Sigma }}^{-1}\right)\\ & +\frac{1}{2}\mathrm{t}\mathrm{r}\left[{\mathrm{E}}_{(\mathit{X}\mid \mathbf{\Phi },\mathbf{\Sigma })}\left\{\mathit{X}(\mathbf{\Phi }-\hat{\mathbf{\Phi }}){\hat{\mathbf{\Sigma }}}^{-1}(\mathbf{\Phi }-\hat{\mathbf{\Phi }}{)}^{\prime }{\mathit{X}}^{\prime }\right\}\right]\\ & =\frac{T}{2}\left\{\mathrm{t}\mathrm{r}\left({\hat{\mathbf{\Sigma }}}^{-1}\mathbf{\Sigma }\right)-\log |{\hat{\mathbf{\Sigma }}}^{-1}\mathbf{\Sigma }|-p\right\}\\ & +\frac{1}{2}\mathrm{t}\mathrm{r}\left\{{\hat{\mathbf{\Sigma }}}^{-1}(\mathbf{\Phi }-\hat{\mathbf{\Phi }}{)}^{\prime }{\mathrm{E}}_{(\mathit{X}\mid \mathbf{\Phi },\mathbf{\Sigma })}({\mathit{X}}^{\prime }\mathit{X}\left)\right(\mathbf{\Phi }-\hat{\mathbf{\Phi }})\right\}.\end{array}$ The result of this derivation can be summarised by the following lemma.

Lemma 2.1

The entropy loss function for the LTS model is (5) $\begin{array}{rl}& \frac{T}{2}\{\mathrm{t}\mathrm{r}\left({\tilde{\mathbf{\Sigma }}}^{-1}\mathbf{\Sigma }\right)-\log |{\tilde{\mathbf{\Sigma }}}^{-1}\mathbf{\Sigma }|-p\}\\ & +\frac{1}{2}\mathrm{t}\mathrm{r}\{{\tilde{\mathbf{\Sigma }}}^{-1}(\mathbf{\Phi }-\tilde{\mathbf{\Phi }}{)}^{\prime }\mathit{W}(\mathbf{\Phi }-\tilde{\mathbf{\Phi }})\},\end{array}$(5) where (6) $\begin{array}{r}\mathit{W}=\frac{1}{T}{\mathrm{E}}_{(\mathit{X}\mid \mathbf{\Phi },\mathbf{\Sigma })}\left({\mathit{X}}^{\prime }\mathit{X}\right).\end{array}$(6)

This Lemma can be proved using algebra similar to that in the iid case. However, there is an important difference. For the LTS model, Bayes estimators involve matrix $\mathit{W}$, a frequentist expectation of ${\mathit{X}}^{\prime }\mathit{X}$ for given parameters $(\mathbf{\Phi },\mathbf{\Sigma })$. For the iid case, no such term is present. The next theorem gives the form of the Bayes estimators under the entropy loss.

Theorem 2.1

The generalised Bayes estimator of $(\mathbf{\Phi },\mathbf{\Sigma })$ under the entropy loss is (7) ${\hat{\mathbf{\Phi }}}_E=\{\mathrm{E}(\mathit{W}\mid \mathit{Y}){\}}^{-1}\mathrm{E}(\mathit{W}\mathbf{\Phi }\mid \mathit{Y}),$(7)

(8) $\begin{array}{r}{\hat{\mathbf{\Sigma }}}_E=\mathrm{E}(\mathbf{\Sigma }\mid \mathit{Y})+\mathrm{E}\{(\mathbf{\Phi }-{\hat{\mathbf{\Phi }}}_E{)}^{\prime }\mathit{W}(\mathbf{\Phi }-{\hat{\mathbf{\Phi }}}_E)\mid \mathit{Y}\}.\end{array}$(8)

The above theorem can be proved similarly as Fact 2.1.

Under the special case with no lag coefficients in the regression, we have $\mathit{W}={\mathrm{E}}_{(\mathit{X}\mid \mathbf{\Phi },\mathbf{\Sigma })}\left({\mathit{X}}^{\prime }\mathit{X}\right)={\mathit{X}}_0^{\prime }{\mathit{X}}_0$, which is not a function of $\mathbf{\Sigma }$ and $\mathbf{\Phi }$. It follows that the Bayes estimator of ${\hat{\mathbf{\Phi }}}_E$ is the posterior mean, as it is for the iid model. This observation is stated in the following remark.

Remark 1

If ${\mathit{B}}_j=0$ for j>0, then $\begin{array}{rl}{\hat{\mathbf{\Phi }}}_E& =\mathrm{E}\{\mathbf{\Phi }\mid \mathit{Y}\},\\ {\hat{\mathbf{\Sigma }}}_E& =\mathrm{E}(\mathbf{\Sigma }\mid \mathit{Y})+\frac{1}{T}{\mathit{X}}_0^{\prime }{\mathit{X}}_0\mathrm{V}\mathrm{a}\mathrm{r}(\mathbf{\Phi }\mid \mathit{Y}).\end{array}$

However, the Bayes estimator for the LTS model is generally different from the iid case. The Bayes estimator ${\hat{\mathbf{\Phi }}}_E$ for the LTS model is not the posterior mean. To compare the estimator ${\hat{\mathbf{\Phi }}}_E$ with the posterior mean, note that in general $\begin{array}{r}{\hat{\mathbf{\Phi }}}_E=\mathrm{E}(\mathbf{\Phi }\mid \mathit{Y})+\{\mathrm{E}(\mathit{W}\mid \mathit{Y}){\}}^{-1}\mathrm{C}\mathrm{o}\mathrm{v}(\mathit{W},\mathbf{\Phi }\mid \mathit{Y}).\end{array}$ Because $\mathit{W}={\mathrm{E}}_{(\mathit{X}\mid \mathbf{\Phi },\mathbf{\Sigma })}\left({\mathit{X}}^{\prime }\mathit{X}\right)$ and $\mathbf{\Phi }$ are likely to be positively correlated, the Bayes estimator of $\mathbf{\Phi }$ under the intrinsic loss is likely to be larger than the posterior mean. It is known that MLE and the posterior mean of $\mathbf{\Phi }$ under a diffuse prior is likely to have a downward bias when the true parameters are closed to random walk, a typical pattern of macroeconomic data. The form of Bayes estimator of $\mathbf{\Phi }$ based on the intrinsic loss is helpful in correcting the bias in the posterior mean.

The estimator in LTS model involves the frequentist expectation $\mathit{W}$. The $\mathit{W}$ matrix depends on specifications of the regressors $\mathit{X}$. If the regressors are specified as functions of lags of $\mathit{Y},$ the computation for $\mathit{W}$ matrix becomes nontrivial.

Using notation in Equation (1(1) $\begin{array}{r}{\mathit{y}}_t^{\prime }={\mathit{x}}_{0t}^{\prime }{\mathit{B}}_0+\sum _{j=1}^L{\mathit{y}}_{t-j}^{\prime }{\mathit{B}}_j+{\mathit{ϵ}}_t,\end{array}$(1) ), the frequentist expectation matrix $\mathit{W}$ can be written as $\begin{array}{r}{\mathrm{E}}_{(\mathit{X}\mid \mathbf{\Phi },\mathbf{\Sigma })}\left({\mathit{X}}^{\prime }\mathit{X}\right)=\left(\begin{array}{cc}{\mathit{X}}_0^{\prime }{\mathit{X}}_0& {\mathrm{E}}_{(\mathit{X}\mid \mathbf{\Phi },\mathbf{\Sigma })}\left({\mathit{X}}_0^{\prime }{\mathit{X}}_1\right)\\ {\mathrm{E}}_{(\mathit{X}\mid \mathbf{\Phi },\mathbf{\Sigma })}\left({\mathit{X}}_1^{\prime }{\mathit{X}}_0\right)& {\mathrm{E}}_{(\mathit{X}\mid \mathbf{\Phi },\mathbf{\Sigma })}\left({\mathit{X}}_1^{\prime }{\mathit{X}}_1\right)\end{array}\right).\end{array}$

For exogenous variables ${\mathit{X}}_0$, there is no need deriving general closed-form expression for the terms in the above matrix as functions of parameter $\mathbf{\Phi }$ and $\mathbf{\Sigma }$. On the other hand, due to the serial correlations of ${\mathit{y}}_t$, computation of ${\mathrm{E}}_{(\mathit{X}\mid \mathbf{\Phi },\mathbf{\Sigma })}\left({\mathit{X}}_1^{\prime }{\mathit{X}}_1\right)$ is not straightforward. In the presence of exogenous variables no analytical expression for $\mathit{W}$ is available. In the following, we discuss approaches to Bayesian estimation under the entropy loss for the general LTS model.

3. Approaches of computing the expectation ${\mathrm{E}}_{(\mathit{X}\mid \mathbf{\Phi },\mathbf{\Sigma })}\left({\mathit{X}}^{\prime }\mathit{X}\right)$

Theorem 2.1 shows that under the entropy loss the Bayes estimator for $(\mathbf{\Phi },\mathbf{\Sigma })$ involves the frequentist expectation $\mathit{W}={\mathrm{E}}_{(\mathit{X}\mid \mathbf{\Phi },\mathbf{\Sigma })}\left({\mathit{X}}^{\prime }\mathit{X}\right)$, and we need to compute the posterior moments $\mathrm{E}(\mathit{W}\mid \mathit{Y})$, $\mathrm{E}(\mathit{W}\mathbf{\Phi }\mid \mathit{Y})$, and $\mathrm{E}({\mathbf{\Phi }}^{\prime }\mathit{W}\mathbf{\Phi }\mid \mathit{Y})$.

The frequentist expectation ${\mathrm{E}}_{(\mathit{X}\mid \mathbf{\Phi },\mathbf{\Sigma })}\left({\mathit{X}}^{\prime }\mathit{X}\right)$ depends on $\mathbf{\Phi }$ and $\mathbf{\Sigma }$. For the LTS model ${\mathrm{E}}_{(\mathit{X}\mid \mathbf{\Phi },\mathbf{\Sigma })}\left({\mathit{X}}^{\prime }\mathit{X}\right)$ does not have an analytical form and needs to be computed numerically for a given ($\mathbf{\Phi },\mathbf{\Sigma }$). We use $\mathit{Y}$ and $\mathit{X}$ to denote observed data in the LTS model (2(2) $\begin{array}{r}\mathit{Y}=\mathit{X}\mathbf{\Phi }+\mathit{ϵ},\end{array}$(2) ). We generate ${\mathit{Y}}^{\ast }$ and ${\mathit{X}}^{\ast }$ from the same model in (2(2) $\begin{array}{r}\mathit{Y}=\mathit{X}\mathbf{\Phi }+\mathit{ϵ},\end{array}$(2) ) given parameters $\mathbf{\Phi }$ and $\mathbf{\Sigma }$, in order to compute ${\mathrm{E}}_{(\mathit{X}\mid \mathbf{\Phi },\mathbf{\Sigma })}\left({\mathit{X}}^{\prime }\mathit{X}\right)$. There is only one observed data set $\mathit{Y}$ and $\mathit{X}$ but there are many sets of generated ${\mathit{Y}}^{\ast }$ and ${\mathit{X}}^{\ast }$. Suppose $\mathbf{\Phi }$ and $\mathbf{\Sigma }$ need to be simulated by an MCMC algorithm, then ${\mathit{Y}}^{\ast }$ and ${\mathit{X}}^{\ast }$ need to be generated for each draw of $\mathbf{\Phi }$ and $\mathbf{\Sigma }$.

One approach to computing ${\mathrm{E}}_{(\mathit{X}\mid \mathbf{\Phi },\mathbf{\Sigma })}\left({\mathit{X}}^{\prime }\mathit{X}\right)$ is straight forward but time-consuming: for each ${\mathbf{\Phi }}_k$ and ${\mathbf{\Sigma }}_k$ drawn in the kth MCMC cycle we generate many sets of ${\mathit{X}}^{\ast }$ and use the average of ${{\mathit{X}}^{\ast }}^{\prime }{\mathit{X}}^{\ast }$ to approximate ${\mathrm{E}}_{(\mathit{X}\mid \mathbf{\Phi },\mathbf{\Sigma })}\left({\mathit{X}}^{\prime }\mathit{X}\right)$. While this approach is possible in theory its high computational cost renders it infeasible in practice. For practical purposes, we must take an alternative approach to compute Bayes estimates.

Fortunately, we have an alternative approach that does not require much additional computational cost beyond simulating $(\mathbf{\Phi },\mathbf{\Sigma })$. Suppose we simulate one set of data ${\mathit{X}}_k^{\ast }$ from the LTS model in each MCMC cycle with simulated parameters of $({\mathbf{\Phi }}_{k-1},{\mathbf{\Sigma }}_{k-1})$, and then simulate the parameters of the next MCMC cycle $({\mathbf{\Phi }}_k,{\mathbf{\Sigma }}_k)$ conditional on both the sample data $\mathit{X}$ and the simulated data ${\mathit{X}}_k^{\ast }$. We will demonstrate that the posterior moments such as $\mathrm{E}(\mathit{W}\mid \mathit{Y})$, $\mathrm{E}(\mathit{W}\mathbf{\Phi }\mid \mathit{Y})$, and $\mathrm{E}({\mathbf{\Phi }}^{\prime }\mathit{W}\mathbf{\Phi }\mid \mathit{Y})$ can be computed through simulated parameters $({\mathbf{\Phi }}_k,{\mathbf{\Sigma }}_k)$ and the jointly simulated data ${\mathit{X}}_k^{\ast }$ (for $k=1,\dots ,M$). The simulated data are in essence latent parameters. They are not the subject of our interests per se but are useful for simulation of parameter of interest (i.e. the frequentist expectation $\mathit{W}$). Data augmentation is not uncommon in Bayesian simulations, but as we noted in the introduction, this data-augmented simulation approach differs from its other uses in the econometrics and statistics literature. One question of practical importance remains though: The number of elements in simulated matrix ${\mathit{X}}_1^{\ast }$ has the dimension of $T\times Lp$, which can be quite large. Do we have to simulate very long Markov chains to assure the averages are good approximates of the posterior mean? Fortunately, our numerical results show that the answer to the question is “no”.

In the following we propose a general algorithm that formalises the data-augmentation idea discussed above.

3.1. A general algorithm using data as latent parameters

Suppose that observed data $\mathit{X}$ has the density $f(\mathit{x}\mid \mathit{\theta })$, where parameter vector $\mathit{\theta }$ is unknown. A prior $\pi \left(\mathit{\theta }\right)$ can be informative or noninformative. Let ${\mathbf{X}}^{\ast }$ be a random vector (or a matrix) with the density $f({\mathbf{x}}^{\ast }\mid \mathit{\theta })$. Let $h\left(\mathit{\theta }\right)$ be a function of the parameters $\mathit{\theta }$. We are interested in the posterior mean of the quantity ${\mathrm{E}}_{({\mathit{X}}^{\ast }\mid \mathit{\theta })}g({\mathbf{X}}^{\ast },h(\mathit{\theta }\left)\right)$ given the data $\mathit{X}$.

Our algorithm is based on the following fact: $\begin{array}{rl}& \mathrm{E}\left\{{\mathrm{E}}_{({\mathbf{X}}^{\ast }\mid \mathit{\theta })}g\right({\mathit{X}}^{\ast },h\left(\mathit{\theta }\right))\mid \mathit{X}\}\\ & =\frac{{\displaystyle \int \{\int g({\mathit{x}}^{\ast },h(\mathit{\theta }\left)\right)f\left({\mathit{x}}^{\ast }\mid \mathit{\theta }\right)\mathrm{d}{\mathit{x}}^{\ast }\}f(\mathit{X}\mid \mathit{\theta })\pi \left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }}}{{\displaystyle \int f(\mathit{X}\mid \mathit{\theta })\pi \left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }}}\\ & =\int \int g({\mathit{x}}^{\ast },h(\mathit{\theta }\left)\right)\pi ({\mathit{x}}^{\ast },\mathit{\theta }\mid \mathit{X})\mathrm{d}{\mathit{x}}^{\ast }\mathrm{d}\mathit{\theta },\end{array}$ where (9) $\begin{array}{r}\pi ({\mathit{x}}^{\ast },\mathit{\theta }\mid \mathit{X})=\frac{f({\mathit{x}}^{\ast }\mid \mathit{\theta })f(\mathit{X}\mid \mathit{\theta })\pi \left(\mathit{\theta }\right)}{{\displaystyle \int \int f({\tilde{\mathit{x}}}^{\ast }\mid \tilde{\mathit{\theta }})f(\mathit{X}\mid \tilde{\mathit{\theta }})\pi \left(\tilde{\mathit{\theta }}\right)\mathrm{d}{\tilde{\mathit{x}}}^{\ast }\mathrm{d}\tilde{\mathit{\theta }}}}.\end{array}$(9) If we have a random sample $({\mathit{X}}_k^{\ast },{\mathit{\theta }}_k),k=1,\dots ,M$, from the joint distribution of (9(9) $\begin{array}{r}\pi ({\mathit{x}}^{\ast },\mathit{\theta }\mid \mathit{X})=\frac{f({\mathit{x}}^{\ast }\mid \mathit{\theta })f(\mathit{X}\mid \mathit{\theta })\pi \left(\mathit{\theta }\right)}{{\displaystyle \int \int f({\tilde{\mathit{x}}}^{\ast }\mid \tilde{\mathit{\theta }})f(\mathit{X}\mid \tilde{\mathit{\theta }})\pi \left(\tilde{\mathit{\theta }}\right)\mathrm{d}{\tilde{\mathit{x}}}^{\ast }\mathrm{d}\tilde{\mathit{\theta }}}}.\end{array}$(9) ), we can estimate $\mathrm{E}\left({\mathrm{E}}_{({\mathit{X}}^{\ast }\mid \mathit{\theta })}\right\{g({\mathit{X}}^{\ast },\mathit{\theta })\left\}h\right(\mathit{\theta })\mid \mathit{X})$ by using the result $\begin{array}{rl}\hat{\mathrm{E}}\left[\mathrm{E}\right\{g({\mathit{X}}^{\ast },h(\mathit{\theta }\left)\right)\mid \mathit{\theta }\}\mid \mathit{X}]& ={\hat{\mathrm{E}}}_{({\mathit{X}}^{\ast },\mathit{\theta }\mid \mathit{X})}\left\{g\right({\mathit{X}}^{\ast },h\left(\mathit{\theta }\right)\left)\right\}\\ & =\frac{1}{M}\sum _{k=1}^{M}g\left({\mathit{X}}_k^{\ast }h\right({\mathit{\theta }}_k\left)\right).\end{array}$

The problem becomes to generate observations from the joint distribution of $({\mathit{X}}^{\ast },\mathit{\theta })$ given the data $\mathit{X}$. For this task the following MCMC method can be used.

Suppose that at the beginning of cycle k we have $({\mathit{X}}_{k-1}^{\ast },{\mathit{\theta }}_{k-1})$.

Simulating full conditional posterior: We sample from $\begin{array}{r}\pi (\mathit{\theta }\mid {\mathit{X}}^{\ast },\mathit{X})\propto f({\mathit{X}}^{\ast }\mid \mathit{\theta })f(\mathit{X}\mid \mathit{\theta })\pi \left(\mathit{\theta }\right).\end{array}$

Step 1. Simulate ${\mathit{X}}_k^{\ast }\sim f({\mathit{x}}^{\ast }\mid {\mathit{\theta }}_{k-1}).$

Step 2. Simulate ${\mathit{\theta }}_k\sim \pi (\mathit{\theta }\mid {\mathit{X}}_k^{\ast },\mathit{X})\propto f({\mathit{X}}_k^{\ast }\mid \mathit{\theta })f(\mathit{X}\mid \mathit{\theta })\pi \left(\mathit{\theta }\right).$

4. Bayesian estimation of $(\mathbf{\Phi },\mathbf{\Sigma })$ in LTS models

4.1. Priors

The Bayes estimator of LTS depends on the prior of $(\mathbf{\Phi },\mathbf{\Sigma })$ We assume prior independence so the prior for $\pi (\mathbf{\Phi },\mathbf{\Sigma })$ is $\pi \left(\mathit{\phi }\right)\pi \left(\mathbf{\Sigma }\right)$, the product of priors for $\mathbf{\Phi }$ and $\mathbf{\Sigma }$.

For estimation of regression coefficient $\mathbf{\Phi }$, a popular informative prior of $\mathit{\phi }=\mathrm{v}\mathrm{e}\mathrm{c}\left(\mathbf{\Phi }\right)$ is the normal distribution, $N({\mathit{\phi }}_0,{\mathit{M}}_0)$, with hyperparameters ${\mathit{\phi }}_0$ and ${\mathit{M}}_0$: (10) $\begin{array}{r}{\pi }_N\left(\mathit{\phi }\right)\propto |{\mathit{M}}_0{|}^{-1/2}\exp \left\{-\frac{1}{2}(\mathit{\phi }-{\mathit{\phi }}_0{)}^{\prime }{\mathit{M}}_0^{-1}(\mathit{\phi }-{\mathit{\phi }}_0)\right\}.\end{array}$(10)

A popular class of non-informative prior on $\mathbf{\Sigma }$ is ${\pi }_b\left(\mathbf{\Sigma }\right)\propto 1/|\mathbf{\Sigma }{|}^{b/2}.$ If b = p + 1, ${\pi }_b\left(\mathbf{\Sigma }\right)$ becomes the Jeffreys prior (see Zellner, 1971)

Ni et al. (2007) examined intrinsic Bayes estimator under prior ${\pi }_b(\mathit{\phi },\mathbf{\Sigma })={\pi }_N\left(\mathit{\phi }\right){\pi }_b\left(\mathbf{\Sigma }\right).$

In the appendix we show posteriors $\pi (\mathbf{\Sigma }\mid {\mathit{X}}^{\ast },\mathit{X})$ and $\pi (\mathbf{\Phi }\mid {\mathit{X}}^{\ast },\mathit{X})$ can be obtained in analytical form.

As mentioned in the introduction, in multiple-parameter settings the Jeffreys prior often has undesirable properties. Bernardo (1979) proposed an approach of deriving a reference prior by breaking a single multiparameter problem into a consecutive series of problems with fewer numbers of parameters. For examples where the reference priors produce more desirable estimates than the Jeffreys priors, see Berger and Bernardo (1992) and Sun and Berger (1998), among others. In estimating the variance-covariance matrix $\mathbf{\Sigma }$ based on an iid random sample from a normal population with known mean, Yang and Berger (1994) re-parameterised matrix $\mathbf{\Sigma }$ as ${\mathit{O}}^{\prime }\mathit{D}O$, where $\mathit{D}$ is a diagonal matrix the elements of which are the eigenvalues of $\mathbf{\Sigma }$ (in increasing or decreasing order) and $\mathit{O}$ is an orthogonal matrix. The following reference prior is derived by giving vectorised $\mathit{D}$ higher priority over vectorised $\mathit{O}$: (11) $\begin{array}{r}{\pi }_R\left(\mathbf{\Sigma }\right)\propto \frac{1}{|\mathbf{\Sigma }|\prod _{1\le i<j\le p}(d_i-d_j)},\end{array}$(11) where $d_1>d_2>\cdots >d_p$ are the eigenvalues of $\mathbf{\Sigma }.$

For numerical and empirical exercise in the is study we use the normal-reference prior ${\pi }_N\left(\mathit{\phi }\right){\pi }_R\left(\mathbf{\Sigma }\right)$. The conditional densities under the normal-reference prior are in the appendix. Ni and Sun (2003) proved that the posteriors of $(\mathbf{\Phi },\mathbf{\Sigma })$ are proper under the normal-reference prior. But the conditional posterior $\pi (\mathbf{\Sigma }\mid \mathbf{\Phi },{\mathit{X}}^{\ast },\mathit{X})$ does not have an analytical form and must be sampled numerically.

4.2. A simulation algorithm for LTS models under the normal-reference prior

We employ an MCMC method to sample from the posterior. In particular, we use the Gibbs sampling method (cf. Gelfand & Smith, 1990). The following algorithm simulates the posteriors of LTS parameters conditional on both the sample and generated data.

Suppose that at cycle k, we have $({\mathbf{\Phi }}_{k-1},{\mathbf{\Sigma }}_{k-1})$ (with an initial draw of $\mathbf{\Sigma }$ and $\mathbf{\Phi }$, e.g. the MLE.)

Algorithm MCMC:

Step 1. Generate ${\mathit{Y}}_k^{\ast }|{\mathbf{\Sigma }}_{k-1},{\mathbf{\Phi }}_{k-1}$.

Simulate ${\mathit{y}}_{k,t}^{{\ast }^{\prime }}\sim N({\mathit{x}}_{0t}^{\prime }{\mathit{B}}_0+\sum _{i=1}^L{\mathit{y}}_{k,t-i}^{{\ast }^{\prime }}{\mathit{B}}_{k-1,i},{\mathbf{\Sigma }}_{k-1}),$ for $t=1,\dots ,T$. Define $\begin{array}{r}{\mathit{Y}}_k^{\ast }=\left(\begin{array}{c}{\mathit{y}}_{k,1}^{{\ast }^{\prime }}\\ ⋮\\ {\mathit{y}}_{k,T}^{{\ast }^{\prime }}\end{array}\right)\text{and}{\mathit{X}}_k^{\ast }=\left(\begin{array}{cccc}{\mathit{x}}_{01}^{\prime }& {\mathit{y}}_{k,-1}^{{\ast }^{\prime }}& \cdots & {\mathit{y}}_{k,-L}^{{\ast }^{\prime }}\\ {\mathit{x}}_{02}^{\prime }& {\mathit{y}}_{k,0}^{{\ast }^{\prime }}& \cdots & {\mathit{y}}_{k,1-L}^{{\ast }^{\prime }}\\ ⋮& ⋮& \ddots & ⋮\\ {\mathit{x}}_{0T}^{\prime }& {\mathit{y}}_{k,T-1}^{{\ast }^{\prime }}& \cdots & {\mathit{y}}_{k,T-L}^{{\ast }^{\prime }}\end{array}\right).\end{array}$

Step 2. Generate ${\mathbf{\Phi }}_k|{\mathbf{\Sigma }}_{k-1},{\mathit{Y}}_k^{\ast },\mathit{Y}$.

Simulate ${\mathit{\phi }}_k=\mathrm{v}\mathrm{e}\mathrm{c}\left({\mathbf{\Phi }}_k\right)\sim N({\mathit{\mu }}_k,{\mathit{V}}_k),$ where (12) $\begin{array}{rl}{\mathit{\mu }}_k& ={\mathit{V}}_k^{\ast }\left[\right\{{\mathbf{\Sigma }}_{k-1}^{-1}\otimes ({\mathit{X}}^{\prime }\mathit{X}\\ & +{{\mathit{X}}_k^{\ast }}^{\prime }{\mathit{X}}_k^{\ast }\left)\right\}{\hat{\mathit{\phi }}}_{\mathrm{M}\mathrm{L}\mathrm{E}k}^{\ast }+{\mathit{M}}_0^{-1}{\mathit{\phi }}_0];\end{array}$(12) (13) $\begin{array}{rl}{\hat{\mathit{\phi }}}_{\mathrm{M}\mathrm{L}\mathrm{E}k}^{\ast }& =\mathrm{v}\mathrm{e}\mathrm{c}\{({\mathit{X}}^{\prime }\mathit{X}+{{\mathit{X}}_k^{\ast }}^{\prime }{\mathit{X}}_k^{\ast }{)}^{-1}({\mathit{X}}^{\prime }\mathit{Y}+{{\mathit{X}}_k^{\ast }}^{\prime }{\mathit{Y}}_k^{\ast })\}.\end{array}$(13) (14) $\begin{array}{rl}{\mathit{V}}_k^{\ast }& =\{{\mathit{M}}_0^{-1}+{\mathbf{\Sigma }}_{k-1}^{-1}\otimes ({\mathit{X}}^{\prime }\mathit{X}+{{\mathit{X}}_k^{\ast }}^{\prime }{\mathit{X}}_k^{\ast }){\}}^{-1}.\end{array}$(14)

Steps 3 to 6 generate ${\mathbf{\Sigma }}_k|{\mathbf{\Sigma }}_{k-1},{\mathbf{\Phi }}_k,{\mathit{Y}}_k^{\ast },\mathit{Y}$.

Step 3: Calculate ${\mathit{S}}_k=S\left({\mathbf{\Phi }}_k\right)+S^{\ast }\left({\mathbf{\Phi }}_k\right)=(\mathit{Y}-\mathit{X}\times {\mathbf{\Phi }}_k{)}^{\prime }(\mathit{Y}-\mathit{X}{\mathbf{\Phi }}_k)+({\mathit{Y}}^{\ast }-{\mathit{X}}^{\ast }{\mathbf{\Phi }}_k{)}^{\prime }({\mathit{Y}}^{\ast }-{\mathit{X}}^{\ast }{\mathbf{\Phi }}_k).$ Decompose ${\mathbf{\Sigma }}_{k-1}=\mathit{O}\mathit{D}{\mathit{O}}^{\prime },$ where $\mathit{O}$ is an orthogonal matrix, $\mathit{D}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(d_1,\dots ,d_p)$ and $d_1>d_2\cdots >d_p.$ Let $d_i^{\mathrm{\#}}=\log \left(d_i\right),$ $\begin{array}{r}{\mathit{D}}^{\mathrm{\#}}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(d_1^{\mathrm{\#}},\dots ,d_p^{\mathrm{\#}})\mathrm{a}\mathrm{n}\mathrm{d}{\mathbf{\Sigma }}_{k-1}^{\mathrm{\#}}=\mathit{O}{\mathit{D}}^{\mathrm{\#}}{\mathit{O}}^{\prime }.\end{array}$

Step 4: Select a random symmetric $p\times p$ matrix $\mathit{V}$, with elements, $v_{ij}=z_{ij}/\sqrt{\sum _{l\le m}{z}_{lm}^2},$ where $z_{ij}\sim N(0,1)$ ($1\le i\le j\le p$, the other elements of $\mathit{V}$ are defined by symmetry).

Step 5: Generate $\lambda \sim N(0,1)$ and set $\mathbf{\Psi }={\mathbf{\Sigma }}_{k-1}^{\mathrm{\#}}+\lambda \mathit{V}.$ Decompose $\mathbf{\Psi }=\mathit{Q}{\mathit{C}}^{\mathrm{\#}}{\mathit{Q}}^{\prime },$ where $\mathit{Q}$ is an orthogonal matrix, ${\mathit{C}}^{\mathrm{\#}}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(c_1^{\mathrm{\#}},\dots ,c_p^{\mathrm{\#}})$ and $c_1^{\mathrm{\#}}>c_2^{\mathrm{\#}}\cdots >c_p^{\mathrm{\#}}.$ Compute $\begin{array}{rl}{\text{β}}_k& =T\sum _{i=1}^p(d_i^{\mathrm{\#}}-c_i^{\mathrm{\#}})\\ & +\frac{1}{2}\mathrm{t}\mathrm{r}\left[\right\{(\mathrm{e}\mathrm{x}\mathrm{p}{\mathbf{\Sigma }}_{k-1}^{\mathrm{\#}}{)}^{-1}-(\exp \mathbf{\Psi }{)}^{-1}\left\}{\mathit{S}}_k\right]\\ & +\sum _{i<j}\log (d_i^{\mathrm{\#}}-d_j^{\mathrm{\#}})-\sum _{i<j}\log (c_i^{\mathrm{\#}}-c_j^{\mathrm{\#}}).\end{array}$

Step 6: Define $\mathit{C}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(\exp (c_1^{\mathrm{\#}}),\dots ,\exp (c_p^{\mathrm{\#}}\left)\right)$ and $\tilde{\mathbf{\Sigma }}=\mathit{Q}\mathit{C}{\mathit{Q}}^{\prime }$. Simulate $u\sim \mathrm{u}\mathrm{n}\mathrm{i}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{m}(0,1)$ and let $\begin{array}{r}{\mathbf{\Sigma }}_k=\left\{\begin{array}{ll}\tilde{\mathbf{\Sigma }},& \mathrm{i}\mathrm{f}u\le min\{1,\exp ({\text{β}}_k\left)\right\},\\ {\mathbf{\Sigma }}_{k-1},& \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}.\end{array}\right.\end{array}$

Note the acceptance probability $\exp \left({\text{β}}_k\right)=\pi (\mathbf{\Psi }|{\mathbf{\Phi }}_k,{\mathit{Y}}_k^{\ast },\mathit{Y})/\pi ({\mathbf{\Sigma }}_{k-1}|{\mathbf{\Phi }}_k,{\mathit{Y}}_k^{\ast },\mathit{Y})$, where the conditional posterior $\pi (\mathbf{\Sigma }|\mathbf{\Phi },{\mathit{Y}}^{\ast },\mathit{Y})$ is given in (A14(A14) $\begin{array}{rl}& \pi ({\mathbf{\Sigma }}^{\mathrm{\#}}|\mathbf{\Phi },\mathit{Y},{\mathit{Y}}^{\ast })=\pi \left({\mathbf{\Sigma }}^{\mathrm{\#}}|\mathit{S}\right(\mathbf{\Phi }),{\mathit{S}}^{\ast }(\mathbf{\Phi }\left)\right)\\ & \propto \frac{\mathrm{e}\mathrm{t}\mathrm{r}[-T|{\mathbf{\Sigma }}^{\mathrm{\#}}|-\frac{1}{2}(\mathrm{e}\mathrm{x}\mathrm{p}{\mathbf{\Sigma }}^{\ast }{)}^{-1}\left\{\mathit{S}\right(\mathbf{\Phi })+{\mathit{S}}^{\ast }(\mathbf{\Phi }\left)\right\}]}{\prod _{i<j}(d_i^{\mathrm{\#}}-d_j^{\mathrm{\#}})},\end{array}$(A14) ). To accelerate the convergence, we repeat Step 6 up to five times until a new candidate is accepted.

4.3. Computing the posterior average loss

From Lemma 2.1, given the estimate $(\hat{\mathbf{\Phi }},\hat{\mathbf{\Sigma }})$, which is computed for a given data sample $\mathit{Y}$, we write the posterior average loss ${\mathrm{E}}_{\left(\right(\mathbf{\Phi },\mathbf{\Sigma }\left)|\mathbf{Y}\mathbf{\right)}}L(\hat{\mathbf{\Phi }},\hat{\mathbf{\Sigma }},\mathbf{\Phi },\mathbf{\Sigma })$ as ${\mathrm{E}}_{(\mathbf{\Phi },\mathbf{\Sigma })|\mathit{Y}}{L}_1(\hat{\mathbf{\Sigma }};\mathbf{\Sigma })+{\mathrm{E}}_{(\mathbf{\Phi },\mathbf{\Sigma })|\mathit{Y}}{L}_2(\hat{\mathbf{\Sigma }},\hat{\mathbf{\Phi }},\mathbf{\Sigma },\mathbf{\Phi })$. We decompose the ${\mathbf{\Sigma }}_k$ in the kth MCMC cycle as ${\mathbf{\Sigma }}_k=\mathit{Q}{\mathit{D}}_k{\mathit{Q}}^{\prime }$, where ${\mathit{D}}_k$ is the diagonal matrix that consists of eigenvalues of ${\mathbf{\Sigma }}_k$: ${\mathit{D}}_k=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mathit{d}}_{k1},{\mathit{d}}_{k2},\dots ,{\mathit{d}}_{kp})$, and $\mathit{Q}$ is an orthogonal matrix with $\mathit{Q}{\mathit{Q}}^{\prime }=\mathit{I}$.

The posterior average loss under the intrinsic loss can be computed using the posterior draws of $({\mathbf{\Phi }}_k,{\mathbf{\Sigma }}_k,{\mathit{X}}_k^{\ast })$ generated by the MCMC procedure $(k=1,2,\dots ,M)$, with (15) $\begin{array}{rl}& {\hat{\mathrm{E}}}_{\left(\right(\mathbf{\Phi },\mathbf{\Sigma }\left)|\mathbf{Y}\mathbf{\right)}}{L}_1(\hat{\mathbf{\Sigma }};\mathbf{\Sigma })\\ & ={\hat{\mathrm{E}}}_{\left(\right(\mathbf{\Phi },\mathbf{\Sigma }\left)|\mathbf{Y}\mathbf{\right)}}\{\frac{T}{2}\left\{\mathrm{t}\mathrm{r}\right({\hat{\mathbf{\Sigma }}}^{-1}\mathbf{\Sigma })-\mathrm{l}\mathrm{o}\mathrm{g}|{\hat{\mathbf{\Sigma }}}^{-1}\mathbf{\Sigma }|-p\}\\ & =\frac{T}{2}\{\mathrm{t}\mathrm{r}({\hat{\mathbf{\Sigma }}}^{-1}\frac{1}{M}\sum _{k=1}^M{\mathbf{\Sigma }}_k)\\ & +\log |\hat{\mathbf{\Sigma }}|-p-\frac{1}{M}\sum _{k=1}^M\sum _{i=1}^p\log |{\mathit{d}}_{ki}|\};\end{array}$(15) (16) $\begin{array}{rl}& {\hat{\mathrm{E}}}_{\left(\right(\mathbf{\Phi },\mathbf{\Sigma }\left)|\mathbf{Y}\mathbf{\right)}}{L}_2(\hat{\mathbf{\Phi }},\hat{\mathbf{\Sigma }};\mathbf{\Phi },\mathbf{\Sigma })\\ & ={\hat{\mathrm{E}}}_{\left(\right(\mathbf{\Phi },\mathbf{\Sigma }\left)|\mathbf{Y}\mathbf{\right)}}\frac{1}{2}\mathrm{t}\mathrm{r}\left[{\hat{\mathbf{\Sigma }}}^{-1}\right\{(\mathbf{\Phi }-\hat{\mathbf{\Phi }}{)}^{\prime }\mathit{W}(\mathbf{\Phi }-\hat{\mathbf{\Phi }}\left)\right\}]\\ & =\frac{1}{2}\{\mathrm{t}\mathrm{r}({\hat{\mathbf{\Sigma }}}^{-1}\frac{1}{M}\sum _{k=1}^M{{\mathbf{\Phi }}_k}^{\prime }{{\mathit{X}}_k^{\ast }}^{\prime }{\mathit{X}}_k^{\ast }{\mathbf{\Phi }}_k)\\ & +\mathrm{t}\mathrm{r}(\hat{\mathbf{\Phi }}{\hat{\mathbf{\Sigma }}}^{-1}{\hat{\mathbf{\Phi }}}^{\prime }\frac{1}{M}\sum _{k=1}^M{{\mathit{X}}_k^{\ast }}^{\prime }{\mathit{X}}_k^{\ast })\}\\ & -\mathrm{t}\mathrm{r}({\hat{\mathbf{\Sigma }}}^{-1}{\hat{\mathbf{\Phi }}}^{\prime }\frac{1}{M}\sum _{k=1}^M{{\mathit{X}}_k^{\ast }}^{\prime }{\mathit{X}}_k^{\ast }{\mathbf{\Phi }}_k),\end{array}$(16) where $\mathit{W}={\mathrm{E}}_{(\mathit{X}\mid \mathbf{\Phi },\mathbf{\Sigma })}\left({\mathit{X}}^{\prime }\mathit{X}\right)$.

Note that all terms in the posterior entropy loss are functions of simulated $\mathbf{\Sigma }$, $\mathbf{\Phi }$, and ${\mathit{X}}^{\ast }$ over the MCMC cycles. The moments of the simulated parameters can be computed in the MCMC cycles, just as the posterior mean, without the need of storage of all of the simulated parameters.

5. A numerical example and an empirical study

5.1. A numerical example

In this section we first simulate data from an LTS model (17) $\begin{array}{r}{\mathit{y}}_t=\mathit{c}+({\mathit{y}}_{t-1},{\mathit{x}}_t)\mathit{B}+{\mathit{ϵ}}_t,\end{array}$(17) for $t=1,\dots ,T$. The dimension of the VAR variable ${\mathit{y}}_t$ is 5. The exogenous variable ${\mathit{x}}_t$ is a scalar representing seasonal cycles, with ${\mathit{x}}_1=0$, ${\mathit{x}}_2=1$, ${\mathit{x}}_3=0$, ${\mathit{x}}_4=-1$, and ${\mathit{x}}_t={\mathit{x}}_{t-4}$ for t>4.

Now we let the true parameters be $\begin{array}{rl}\mathbf{\Sigma }& =\left(\begin{array}{ccccc}1& 0& 0& 0& 0\\ 0& 2& 0& 0& 0\\ 0& 0& 3& 0& 0\\ 0& 0& 0& 4& 0\\ 0& 0& 0& 0& 5\end{array}\right),\\ \mathit{B}& =\left(\begin{array}{ccccc}0& 0& 0& 0& 0\\ 0.5& 0& 0& 0& 0\\ 0& 0.5& 0& 0& 0\\ 0& 0& 0.5& 0& 0\\ 0& 0& 0& 0.5& 0\\ 0& 0& 0& 0& 0.5\\ 0.1& 0.2& 0.3& 0.4& 0.5\end{array}\right).\end{array}$

The last row of matrix $\mathit{B}$ are the parameters of the seasonal dummies. The discussion in Section 2 shows with this parameter setting there is no closed-form expression of the frequentist expectation $\text{E}({\mathit{X}}^{\prime }\mathit{X}|\mathbf{\Sigma },\mathbf{\Phi })$ and we need to simulate ${\mathrm{E}}_{(\mathit{X}\mid \mathbf{\Phi },\mathbf{\Sigma })}\left({\mathit{X}}^{\prime }\mathit{X}\right)$ using the data-augmentation algorithms proposed in Section 4.

We generate one data sample (T) of 100 observations from the LTS model with the above parameters for $\mathbf{\Sigma }$ and $\mathbf{\Phi }$. The MLE of the parameters are

$\begin{array}{rl}{\hat{\mathbf{\Sigma }}}_{\mathrm{M}\mathrm{L}\mathrm{E}}& =\left(\begin{array}{ccccc}0.93201& 0.12358& -0.08462& 0.27246& -0.11782\\ 0.12358& 1.76018& 0.02217& -0.16508& 0.26475\\ -0.08462& 0.02217& 2.44980& -0.07874& -0.01940\\ \$0.27246\$& \$-\$0.16508& \$-\$0.07874& 3.44815& -0.80056\\ \$-\$0.11782& 0.26475& -0.01940& -0.80056& 4.03612\end{array}\right),\\ {\hat{\mathbf{\Phi }}}_{\mathrm{M}\mathrm{L}\mathrm{E}}& =\left(\begin{array}{ccccc}-0.05633& 0.14700& 0.01158& 0.16497& 0.26224\\ 0.51936& 0.24397& 0.21705& -0.1209& -0.25149\\ 0.0214& 0.41804& -0.03495& 0.06221& 0.28700\\ -0.04826& -0.03265& 0.50279& -0.03351& -0.10127\\ -0.06176& 0.01439& 0.14174& 0.37051& -0.01758\\ -0.0194& -0.0014& -0.08389& 0.01734& 0.39826\\ 0.15900& 0.22462& -0.21547& -0.01671& 1.06891\end{array}\right).\end{array}$

We conduct the simulation with a diffuse prior on $\mathbf{\Phi }$ and the Yang–Berger reference prior on $\mathbf{\Sigma }$. The length of MCMC cycles set at 100,000. The intrinsic Bayes estimates are

(18) $\begin{array}{rl}{\hat{\mathbf{\Sigma }}}_E& =\left(\begin{array}{ccccc}1.25871& 0.17644& -0.08826& 0.26711& -0.13247\\ 0.17644& 2.33691& 0.00497& -0.15278& 0.28849\\ -0.08826& 0.00497& 2.97993& -0.03929& -0.07036\\ 0.26711& -0.15278& -0.03929& 3.76133& -0.74181\\ -0.13247& 0.28849& -0.07036& -0.74181& 4.57260\end{array}\right),\end{array}$(18) (19) $\begin{array}{rl}{\hat{\mathbf{\Phi }}}_E& =\left(\begin{array}{ccccc}-0.05655& 0.09316& -0.00092& 0.11812& 0.16943\\ 0.46112& 0.18679& 0.13921& -0.09552& -0.13973\\ 0.02891& 0.37594& -0.01426& 0.05658& 0.18949\\ -0.04417& -0.03156& 0.45267& -0.0236& -0.10665\\ -0.05809& 0.02705& 0.11909& 0.33560& 0.02602\\ -0.02400& 0.00105& -0.08313& 0.02921& 0.34976\\ 0.11188& 0.08832& -0.08980& 0.04434& 0.34897\end{array}\right).\end{array}$(19)

The acceptance rate for the Metropolis step employed for sampling of $\mathbf{\Sigma }$ from the posterior conditional on other parameters and data is $24\mathrm{\%}$.

Now we compare the Bayes estimator with the MLE and posterior mean. The posterior mean obtained using Option 1 is

(20) $\begin{array}{rl}{\hat{\mathbf{\Sigma }}}_{\mathrm{M}\mathrm{e}\mathrm{a}\mathrm{n}}& =\left(\begin{array}{ccccc}\$1.17992\$& 0.16657& -0.08449& 0.25618& -0.12700\\ 0.16657& 2.20489& 0.00476& -0.14659& 0.27587\\ -0.08449& 0.00476& 2.81882& -0.03801& -0.06557\\ 0.25618& -0.14659& -0.03801& 3.57382& -0.71561\\ -0.12700& 0.27587& -0.06557& -0.71561& 4.35069\end{array}\right),\end{array}$(20) (21) $\begin{array}{rl}{\hat{\mathbf{\Phi }}}_{\mathrm{M}\mathrm{e}\mathrm{a}\mathrm{n}}& =\left(\begin{array}{ccccc}-0.0599& 0.09731& -0.00193& 0.12276& 0.17558\\ 0.43287& 0.18741& 0.14185& -0.09236& -0.13907\\ 0.02772& 0.35030& -0.01270& 0.05628& 0.19118\\ -0.04650& -0.03116& 0.42790& -0.02733& -0.10462\\ -0.05718& 0.02602& 0.12006& 0.31275& 0.02437\\ -0.02373& -0.00041& -0.08285& 0.02851& 0.32792\\ 0.11167& 0.08778& -0.08838& 0.04414& 0.34923\end{array}\right).\end{array}$(21)

It is known that the posterior mean of $\mathbf{\Sigma }$ minimises the expected posterior loss of $L_1$. But it does not minimise the intrinsic loss because the estimator of $\mathbf{\Sigma }$ also influences the weight of $\mathbf{\Phi }$ related loss in $L_2$. Table 1 reports the average posterior loss of the MLE, posterior mean, and the Bayes estimator for the data sample generated in the example. The Bayes estimator improves $L_2$-related risk with a tradeoff of larger $L_1$-related risk. Table 1 shows that the Bayes estimator induces lower posterior risk than the posterior mean by making the $L_2$-related risk substantially lower and the $L_1$-related risk only slightly higher. Both posterior mean and the Bayes estimator dominate the MLE.

Table 1. Posterior average loss of the estimates in the Example.

	E $L_1$	E $L_2$	E L =E ( $L_1$+$L_2$)
MLE (${\mathbf{\Phi }}_{\mathrm{M}\mathrm{L}\mathrm{E}},{\mathbf{\Sigma }}_{\mathrm{M}\mathrm{L}\mathrm{E}})$	$11.2940$	$22.6213$	$33.9154$
Posterior mean (${\mathbf{\Phi }}_{\mathrm{M}\mathrm{e}\mathrm{a}\mathrm{n}},{\mathbf{\Sigma }}_{\mathrm{M}\mathrm{e}\mathrm{a}\mathrm{n}})$	$7.5414$	$17.7350$	$25.2764$
Intrinsic Bayes estimator (${\mathbf{\Phi }}_E,{\mathbf{\Sigma }}_E)$	$7.9295$	$13.6053$	$21.5348$

Notes: The loss functions are $L_1(\hat{\mathbf{\Sigma }},\mathbf{\Sigma })=\left(T/2\right)\left\{\mathrm{t}\mathrm{r}\right(\mathbf{\Sigma }{\tilde{\mathbf{\Sigma }}}^{-1})-\log |\mathbf{\Sigma }{\tilde{\mathbf{\Sigma }}}^{-1}|-p\},$ and $L_2(\hat{\mathbf{\Sigma }},\hat{\mathbf{\Phi }},\mathbf{\Sigma },\mathbf{\Phi })=\frac{1}{2}\mathrm{t}\mathrm{r}\left[{\hat{\mathbf{\Sigma }}}^{-1}\right\{(\mathbf{\Phi }-\hat{\mathbf{\Phi }}{)}^{\prime }{\mathrm{E}}_{(\mathit{X}\mid \mathbf{\Phi },\mathbf{\Sigma })}({\mathit{X}}^{\prime }\mathit{X}\left)\right(\mathbf{\Phi }-\hat{\mathbf{\Phi }}\left)\right\}].$ The posterior average loss for $L_1$ and $L_2$ are computed by (15(15) $\begin{array}{rl}& {\hat{\mathrm{E}}}_{\left(\right(\mathbf{\Phi },\mathbf{\Sigma }\left)|\mathbf{Y}\mathbf{\right)}}{L}_1(\hat{\mathbf{\Sigma }};\mathbf{\Sigma })\\ & ={\hat{\mathrm{E}}}_{\left(\right(\mathbf{\Phi },\mathbf{\Sigma }\left)|\mathbf{Y}\mathbf{\right)}}\{\frac{T}{2}\left\{\mathrm{t}\mathrm{r}\right({\hat{\mathbf{\Sigma }}}^{-1}\mathbf{\Sigma })-\mathrm{l}\mathrm{o}\mathrm{g}|{\hat{\mathbf{\Sigma }}}^{-1}\mathbf{\Sigma }|-p\}\\ & =\frac{T}{2}\{\mathrm{t}\mathrm{r}({\hat{\mathbf{\Sigma }}}^{-1}\frac{1}{M}\sum _{k=1}^M{\mathbf{\Sigma }}_k)\\ & +\log |\hat{\mathbf{\Sigma }}|-p-\frac{1}{M}\sum _{k=1}^M\sum _{i=1}^p\log |{\mathit{d}}_{ki}|\};\end{array}$(15) ) and (16(16) $\begin{array}{rl}& {\hat{\mathrm{E}}}_{\left(\right(\mathbf{\Phi },\mathbf{\Sigma }\left)|\mathbf{Y}\mathbf{\right)}}{L}_2(\hat{\mathbf{\Phi }},\hat{\mathbf{\Sigma }};\mathbf{\Phi },\mathbf{\Sigma })\\ & ={\hat{\mathrm{E}}}_{\left(\right(\mathbf{\Phi },\mathbf{\Sigma }\left)|\mathbf{Y}\mathbf{\right)}}\frac{1}{2}\mathrm{t}\mathrm{r}\left[{\hat{\mathbf{\Sigma }}}^{-1}\right\{(\mathbf{\Phi }-\hat{\mathbf{\Phi }}{)}^{\prime }\mathit{W}(\mathbf{\Phi }-\hat{\mathbf{\Phi }}\left)\right\}]\\ & =\frac{1}{2}\{\mathrm{t}\mathrm{r}({\hat{\mathbf{\Sigma }}}^{-1}\frac{1}{M}\sum _{k=1}^M{{\mathbf{\Phi }}_k}^{\prime }{{\mathit{X}}_k^{\ast }}^{\prime }{\mathit{X}}_k^{\ast }{\mathbf{\Phi }}_k)\\ & +\mathrm{t}\mathrm{r}(\hat{\mathbf{\Phi }}{\hat{\mathbf{\Sigma }}}^{-1}{\hat{\mathbf{\Phi }}}^{\prime }\frac{1}{M}\sum _{k=1}^M{{\mathit{X}}_k^{\ast }}^{\prime }{\mathit{X}}_k^{\ast })\}\\ & -\mathrm{t}\mathrm{r}({\hat{\mathbf{\Sigma }}}^{-1}{\hat{\mathbf{\Phi }}}^{\prime }\frac{1}{M}\sum _{k=1}^M{{\mathit{X}}_k^{\ast }}^{\prime }{\mathit{X}}_k^{\ast }{\mathbf{\Phi }}_k),\end{array}$(16) ).

5.2. An empirical study: seasonal effects in a macroeconomic model

We now turn to an empirical application of the Bayesian estimation under the entropy loss. We estimate an LTS consisting of seasonal dummies and four macroeconomic variables: the return of Standard and Poor 500 stock price index (which represents weighted stock prices of large companies), the 3-month Treasury Bill rate, the growth rate of payroll (including government jobs as well as private sector jobs), and the growth rate of industrial production (in that order). All series are measured in percentage terms. The series are monthly data from 1970:1 to 2002:12 and are not seasonally-adjusted. The payroll data are obtained from the Bureau of Labor Statistics, the rest series are obtained from the Federal Reserve Board. There are 12 column dummy variables, representing January to December.

The role of seasonal fluctuations in business cycles has been noted by a number of economists. Barky and Miron (1989) argued that for the U.S. economy the characteristics of seasonal fluctuations are similar to the conventional characterisations of business cycles. Cecchetti et al. (1997) estimated production function of various industries based on how their responses to seasonal shocks vary with the state of business cycle. Ghysels (1988) showed that univariate seasonal adjustment of endogenous variables is not harmless because information on the interactions among endogenous variables will be lost. Miron and Beaulieu (1996) provided a survey on econometric and economic issues on understanding business cycles through seasonal fluctuations. In the finance literature, numerous studies argue that stock returns appear to have a seasonal components. Rozeff and Kinney (1976) documented a celebrated “turn of the year” effect, which refers to the seemingly abnormally high returns in January and July, especially for stocks with small-market capitalisations. A number of theories have been developed to explain the phenomenon. Reinganum (1983) attributed the high stock return in January to the end of year tax-loss selling in December. Chang and Pinegar (1989) found that industrial production trails the seasonal movement of stock returns by one month. The reported point estimates of the seasonal effects in the literature are model dependent and based on OLS or MLE. We will estimate the seasonal effects. Our primary interest lies in comparison of posterior mean with Bayes estimate under the entropy loss.

In this section, we employ an LTS model (22) $\begin{array}{r}{\mathit{y}}_t^{\prime }=({\mathit{x}}_t^{\prime },{\mathit{y}}_{t-1}^{\prime },\dots ,{\mathit{y}}_{t-L}^{\prime })\mathbf{\Phi }+{\mathit{ϵ}}_t,\end{array}$(22) for $t=1,\dots ,T.$ The dimension of the VAR variable ${\mathit{y}}_t$ is four. The exogenous variable ${\mathit{x}}_{0t}=(x_{0t,1},\dots ,x_{0t,12})$ is a 12-dimensional vector representing seasonal cycles, with $x_{0i,1}$ equals 1 if period i is January, and 0 otherwise; $x_{0i,12}=1$ if period i is December and 0 otherwise.

Based on the Schwarz criterion, for each sample period the lag length L of the LTS is 2. The Yang–Berger reference prior is applied to the covariance matrix $\mathbf{\Sigma }$. The prior for the LTS coefficient $\mathit{\phi }$ is a rather diffuse $N(\mathbf{0},{\mathit{M}}_0)$. Here ${\mathit{M}}_0$ is a diagonal matrix with 10.0 being the diagonal element for parameters corresponding to the dummy variables and 2.0 being the diagonal element for parameters corresponding to the lag coefficients. We draw the posterior from M MCMC cycles after $0.1\times M$ burn-in runs. The MCMC length M is set at 50, 000, 100, 000, and 1, 000, 000.

Under algorithm MCMC, reducing the length of MCMC cycles to 100, 000 or 50, 000 from 1, 000, 000 makes little difference. The MLE, posterior mean, and Bayes estimate of the covariance matrix $\mathbf{\Sigma }$ are as follows. As dictated by the theoretical result, the Bayes estimate under the entropy loss, ${\mathbf{\Sigma }}_E$ is larger than the posterior mean ${\mathbf{\Sigma }}_{\mathrm{M}\mathrm{e}\mathrm{a}\mathrm{n}}$. $\begin{array}{r}{\hat{\mathbf{\Sigma }}}_{\mathrm{M}\mathrm{L}\mathrm{E}}=\left(\begin{array}{cccc}19.775& -0.335& -0.094& 0.059\\ -0.335& 0.254& 0.055& 0.011\\ -0.094& 0.055& 0.937& 0.070\\ 0.059& 0.011& 0.070& 0.044\end{array}\right).\end{array}$ With M = 1, 000, 000, the posterior mean and entropy-based Bayes estimates are $\begin{array}{rl}{\hat{\mathbf{\Sigma }}}_{\mathrm{M}\mathrm{e}\mathrm{a}\mathrm{n}}& =\left(\begin{array}{cccc}20.908& -0.354& -0.100& 0.063\\ -0.354& 0.271& 0.058& 0.012\\ -0.100& 0.058& 0.996& 0.075\\ 0.063& 0.012& 0.075& 0.048\end{array}\right).\\ {\hat{\mathbf{\Sigma }}}_E& =\left(\begin{array}{cccc}22.729& -0.384& -0.101& 0.078\\ -0.384& 0.295& 0.057& 0.012\\ -0.101& 0.057& 1.074& 0.081\\ 0.078& 0.012& 0.081& 0.052\end{array}\right).\end{array}$ The posterior standard deviations of the elements of the covariance matrix are $\begin{array}{r}\left(\begin{array}{cccc}1.538& 0.124& 0.235& 0.052\\ 0.124& 0.020& 0.027& 0.006\\ 0.235& 0.027& 0.074& 0.012\\ 0.052& 0.006& 0.012& 0.003\end{array}\right).\end{array}$ The difference between estimates ${\hat{\mathbf{\Sigma }}}_{\mathrm{M}\mathrm{e}\mathrm{a}\mathrm{n}}$ and ${\hat{\mathbf{\Sigma }}}_E$ is large relative to the posterior standard deviations.

The above point estimates and standard deviations are similar to those with M = 50, 000. The MCMC algorithm yields posteriors with few outliers. This is because in Step 2 the MCMC algorithm, $\mathit{\phi }$ is generated from an average of sample data and generated data, instead of sample data alone. A few outliers in the posterior may affect the posterior mean slightly but can change the posterior risk and the Bayesian estimate substantially because the few explosive parameters carry disproportionately large weights in the posterior average loss.

The intrinsic Bayes estimate dominates the posterior mean by a large margin in terms of posterior expected loss. The large difference in the posterior expected loss is mainly due to the difference in the $\mathbf{\Phi }$-related risk, i.e. the quadratic term in (5(5) $\begin{array}{rl}& \frac{T}{2}\{\mathrm{t}\mathrm{r}\left({\tilde{\mathbf{\Sigma }}}^{-1}\mathbf{\Sigma }\right)-\log |{\tilde{\mathbf{\Sigma }}}^{-1}\mathbf{\Sigma }|-p\}\\ & +\frac{1}{2}\mathrm{t}\mathrm{r}\{{\tilde{\mathbf{\Sigma }}}^{-1}(\mathbf{\Phi }-\tilde{\mathbf{\Phi }}{)}^{\prime }\mathit{W}(\mathbf{\Phi }-\tilde{\mathbf{\Phi }})\},\end{array}$(5) ). This difference in risk is approximately $\begin{array}{rl}& \frac{1}{2}\mathrm{t}\mathrm{r}[{{\hat{\mathbf{\Sigma }}}_E}^{-1}\mathrm{E}\left\{\right(\mathbf{\Phi }-{\hat{\mathbf{\Phi }}}_E{)}^{\prime }\mathit{W}\\ & \times (\mathbf{\Phi }-{\hat{\mathbf{\Phi }}}_E)\mid \mathit{Y}\left\}{\hat{\mathbf{\Sigma }}}_E^{-1}\right({\hat{\mathbf{\Sigma }}}_E-{\hat{\mathbf{\Sigma }}}_{\mathrm{M}\mathrm{e}\mathrm{a}\mathrm{n}})],\end{array}$ which is proportional to the frequentist expectation $\mathit{W}$, and the latter is comparable to ${\mathit{X}}^{\prime }\mathit{X}$. ${\mathit{X}}^{\prime }\mathit{X}$ is quite large in this application, largely due to the strong serial correlation of the 3-month Tbill rates. As a result, with a larger $\hat{\mathbf{\Sigma }}$ the Bayesian estimate substantially reduces the posterior risk, compared with the posterior mean. Simulation with M = 1, 000, 000 shows that the posterior average loss of the posterior mean estimate $\left(3230.6\right)$ is larger than that of the intrinsic Bayes estimate $\left(69.7\right)$. The lower overall posterior average loss of the intrinsic Bayes estimate is achieved by substantially lowering the risk of the quadratic part, from 3225.3 for the posterior mean to 61.8. The first term of the loss related to $\mathbf{\Sigma }$ of the intrinsic Bayes estimate is slightly larger $\left(7.9\right)$ compared to that under the posterior mean estimate $\left(5.3\right)$. As is noted earlier, the intrinsic Bayes estimate improves $\mathbf{\Phi }$-related loss with a tradeoff of larger $\mathbf{\Sigma }$-related loss. The empirical result shows that the Bayesian estimate induces lower posterior average loss than the posterior mean by making the $\mathbf{\Phi }$-related loss substantially lower and the $\mathbf{\Sigma }$-related loss only slightly higher.

We now turn to compare the estimates of the regression coefficient $\mathbf{\Phi }$. Table 2 reports the MLE, posterior mean, the intrinsic Bayes estimate, and posterior standard deviations. As a consequence of applying a rather diffuse prior, the posterior mean ${\mathbf{\Phi }}_{\mathrm{M}\mathrm{e}\mathrm{a}\mathrm{n}}$ is quite similar to the MLE. For stock returns of Standard and Poor 500 index, MLE and the posterior mean estimate indicate a moderate positive seasonal factor in January, which is smaller than the seasonal factor of March, April, October, November, and December. Most surprisingly, October registers the largest seasonal gain, despite the fact that the sample included the 1987 October sell off. The data of recent years suggest that the estimates of seasonality in large capitalisation stock returns are quite sensitive to the sample period and regression model. In comparison to the MLE and the posterior mean, the intrinsic Bayesian estimate shows a smaller January effect and much smaller end-of-the-year positive seasonal returns. The sum of the seasonal coefficients of the MLE is above $11\mathrm{\%}$ while that of the intrinsic Bayesian estimates is about half as much. The large discrepancy between the posterior mean and the intrinsic Bayesian estimate s casts doubt on the robustness of the seasonality of returns of Standard and Poor 500 stock.

Table 2. Estimates of three equations.

	Stock return equation			IP growth equation			Payroll growth equation
	MLE	PMean (SD)	BE	MLE	PMean(SD)	BE	MLE	PMean (SD)	BE
Jan	1.507	1.469 (1.093)	0.986	0.701	0.701 (0.243)	0.680	−2.019	−2.019 (0.053)	−2.054
Feb	0.737	0.743 (2.344)	0.784	3.980	3.961 (0.532)	3.977	0.842	0.840 (0.117)	0.892
Mar	1.780	1.749 (2.355)	1.662	0.179	0.186 (0.537)	0.123	0.912	0.912 (0.118)	0.955
Apr	1.534	1.461 (1.170)	0.905	−1.286	−1.281 (0.262)	−1.416	0.518	0.518 (0.057)	0.466
May	0.370	0.304 (1.381)	−0.356	−0.138	−0.134 (0.312)	−0.226	0.485	0.485 (0.068)	0.414
June	0.242	0.195 (1.412)	−0.444	2.092	2.095 (0.318)	2.003	0.455	0.455 (0.070)	0.381
July	0.637	0.585 (1.338)	−0.019	−4.589	−4.584 (0.298)	−4.759	−1.284	−1.284 (0.066)	−1.352
Aug	0.822	0.753 (1.830)	0.485	5.270	5.257 (0.414)	5.251	0.262	0.261 (0.091)	0.255
Sept	−1.881	−1.846 (1.798)	−1.964	1.021	1.025 (0.405)	0.998	1.040	1.041 (0.089)	1.049
Oct	2.162	2.072 (1.377)	1.583	−0.409	−0.403 (0.308)	−0.614	0.155	0.155 (0.067)	0.106
Nov	1.891	1.815 (1.240)	1.271	−1.332	−1.33 (0.278)	−1.464	−0.011	−0.011 (0.061)	−0.075
Dec	1.473	1.419 (1.149)	0.884	−0.792	−0.792 (0.257)	−0.843	0.017	0.017 (0.056)	−0.031
stock1	−0.012	−0.011 (0.052)	0.025	0.013	0.013 (0.011)	0.012	0.000	0.000 (0.002)	0.001
int1	−0.581	−0.580 (0.384)	−0.518	0.152	0.153 (0.084)	0.132	0.030	0.030 (0.018)	0.025
ip1	−0.198	−0.201 (0.243)	−0.208	0.052	0.052 (0.054)	0.076	0.021	0.021 (0.012)	0.021
emp1	0.292	0.310 (1.017)	0.543	1.072	1.063 (0.231)	1.100	0.218	0.217 (0.051)	0.255
stock2	−0.040	−0.039 (0.052)	−0.002	0.009	0.009 (0.011)	0.009	0.000	0.000 (0.002)	0.000
int2	0.589	0.594 (0.386)	0.587	−0.207	−0.207 (0.085)	−0.176	−0.036	−0.036 (0.019)	−0.028
ip2	−0.382	−0.373 (0.241)	−0.365	−0.087	−0.087 (0.053)	−0.064	0.046	0.046 (0.012)	0.045
emp2	0.008	0.019 (0.973)	0.198	−0.038	−0.036 (0.222)	−0.006	0.183	0.183 (0.049)	0.219

Notes: stock lag 1, int lag 1, ip lag 1, and emp lag 1 represent the first lag of stock return, the first lag 3-month Tbill rate, the first lag of industrial production growth, and the first lag of private sector payroll growth. The lag 2 variables are the second lag of the corresponding variables. Estimates are obtained with MCMC cycle size M = 1, 000, 000.

Compared to the stock return, the seasonality of industrial production growth rate is much more robust. The most distinct pattern is a steep decline in July followed by a surge in August; then weakness in the end of the year precedes a strong rebound in February. The strong showing of industrial production in February and August is consistent with the pattern reported in Chang and Pinegar (1989) while the predicted industrial production by Standard and Poor stock returns is quite small. The magnitude of the seasonal effects by the entropy-loss-based Bayesian estimates is on average slightly larger than that of the posterior mean.

Lastly, we examine the estimates of the employment growth rate equation. The most prominent seasonal patterns are the decline in January followed by a rebound in February and March and the weakness in July followed by a recovery in September. The estimated seasonality in payroll growth is somewhat different from that of the industrial production growth. Note that in year 2002 about $83\mathrm{\%}$ of the payroll consists of service sector jobs while $85\mathrm{\%}$ of industrial production concerns the manufacture sector. The subject of interest is the point estimates. Similar to the industrial production equation, the entropy-loss-based Bayes estimates for the payroll growth rate equation are similar to the posterior mean.

In summary, the Bayesian estimates based on the entropy loss show a qualitatively similar seasonal pattern to that of the posterior mean estimates for industrial production and employment growth but a distinctly different one for stock returns. The posterior average loss of the Bayes estimates with respect to the entropy loss is substantially smaller than that of the posterior mean.

6. Concluding remarks

In this paper we investigate properties of Bayes estimators of LTS model $(\mathbf{\Phi },\mathbf{\Sigma })$ derived from the entropy loss function. These estimators are distinctly different from the multivariate iid model because of the serial correlation of the time series variables. Bayesian computation under the entropy loss requires simulating a frequentist moment of the regressors. We propose a data-augmenting algorithm for simulation of posteriors and computation of Bayes estimators under the entropy loss and a normal-reference prior. The algorithm that draws from the full conditional posterior is shown to be quite efficient. A novel approach taken in this paper concerns generating data in an MCMC as latent parameters. This idea may be useful for other contexts for simulating complicated posterior moments. Our empirical application to a macroeconomic problem shows that the Bayes estimates under the entropy loss can differ substantially from the posterior mean.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Notes on contributors

Shawn Ni

Dr. Shawn Ni holds a PhD in Economics from University of Minnesota. He is currently Middlebush Professor of Economics and Adjunct Professor of Statistics at University of Missouri-Columbia. He conducts research on a wide range of empirical economics topics and Bayesian statistics.

Dongchu Sun

Dr. Dongchu Sun holds a PhD in Statistics from Purdue University. He is a research professor of statistics at the University of Nebraska-Lincoln and East China Normal University. His research interests includes Bayesian analysis, small area estimation, decision theory, business and econometrics, space-time and longitudinal models, and smoothing splines.

Appendix. Proof of Fact 2.1 and posterior properties

A.1. Proof of Fact  2.1

Proof.

Let $(\tilde{\mathit{\mu }},\tilde{\mathbf{\Sigma }})$ denote an arbitrary estimator of $(\mathit{\mu },\mathbf{\Sigma })$. For the entropy loss function L and posterior $\pi (\mathit{\mu },\mathbf{\Sigma }\mid \mathit{Y})$, the expected posterior loss is $\begin{array}{rl}R(\tilde{\mathit{\mu }},\tilde{\mathbf{\Sigma }}\mid \mathit{Y})& =\frac{T}{2}\mathrm{E}\{(\tilde{\mathit{\mu }}-\mathit{\mu }{)}^{\prime }{\tilde{\mathbf{\Sigma }}}^{-1}(\tilde{\mathit{\mu }}-\mathit{\mu })\\ & +\mathrm{t}\mathrm{r}\left(\mathbf{\Sigma }{\tilde{\mathbf{\Sigma }}}^{-1}\right)-\log |\mathbf{\Sigma }{\tilde{\mathbf{\Sigma }}}^{-1}|-p\}.\end{array}$ The Bayes estimator, which minimises the expected posterior loss, can be derived through conditions on first-order derivatives. Note that because $\tilde{\mathbf{\Sigma }}$ is symmetric, $\begin{array}{r}\frac{\mathrm{\partial }R(\tilde{\mathit{\mu }},\tilde{\mathbf{\Sigma }})}{\mathrm{\partial }\tilde{\mathit{\mu }}}=2\mathrm{E}\left(\tilde{\mathbf{\Sigma }}(\tilde{\mathit{\mu }}-\mathit{\mu })\right).\end{array}$ Let the derivative be 0 yields ${\hat{\mathit{\mu }}}_{\mathrm{i}\mathrm{i}\mathrm{d}}=\mathrm{E}(\mathit{\mu }\mid \mathit{Y})$.

The following identities are known (e.g. Harville, 1998, p. 327) for symmetric matrices $\mathit{A}$ and $\mathit{B}$, $\begin{array}{rl}\frac{\mathrm{\partial }\log \left(|\mathit{A}|\right)}{\mathrm{\partial }\mathit{A}}& =2{\mathit{A}}^{-1}-\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left({\mathit{A}}^{-1}\right),\\ \frac{\mathrm{\partial }\mathrm{t}\mathrm{r}\left(\mathit{A}{\mathit{B}}^{-1}\right)}{\mathrm{\partial }\mathit{B}}& =-2\left({\mathit{B}}^{-1}\mathit{A}{\mathit{B}}^{-1}\right)+\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left({\mathit{B}}^{-1}\mathit{A}{\mathit{B}}^{-1}\right).\end{array}$ Here $\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\mathit{A}\right)$ is a diagonal matrix, whose diagonal elements are these from $\mathit{A}$. Using the conclusion that the estimator for $\mathit{\mu }$ is the posterior mean, we have $\mathrm{E}\left\{\right({\hat{\mathit{\mu }}}_{\mathrm{i}\mathrm{i}\mathrm{d}}-\mathit{\mu }\left)\right({\hat{\mathit{\mu }}}_{\mathrm{i}\mathrm{i}\mathrm{d}}-\mathit{\mu }{)}^{\prime }\mid \mathit{Y}\}=\mathrm{v}\mathrm{a}\mathrm{r}(\mathit{\mu }\mid \mathit{Y})$. Taking this result to the derivative with respect to $\tilde{\mathbf{\Sigma }},$ we have $\begin{array}{rl}& \frac{\mathrm{\partial }R({\hat{\mathit{\mu }}}_{\mathrm{i}\mathrm{i}\mathrm{d}},\tilde{\mathbf{\Sigma }})}{\mathrm{\partial }\tilde{\mathbf{\Sigma }}}\\ & =\mathrm{E}\frac{T}{2}\{-2{\tilde{\mathbf{\Sigma }}}^{-1}\mathrm{v}\mathrm{a}\mathrm{r}\left(\mathit{\mu }\mid \mathit{Y}\right){\tilde{\mathbf{\Sigma }}}^{-1}+\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left({\tilde{\mathbf{\Sigma }}}^{-1}\mathrm{v}\mathrm{a}\mathrm{r}\right(\mathit{\mu }\mid \mathit{Y}\left){\tilde{\mathbf{\Sigma }}}^{-1}\right)\\ & -2{\tilde{\mathbf{\Sigma }}}^{-1}\mathbf{\Sigma }{\tilde{\mathbf{\Sigma }}}^{-1}+\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left({\tilde{\mathbf{\Sigma }}}^{-1}\mathbf{\Sigma }{\tilde{\mathbf{\Sigma }}}^{-1}\right)+2{\tilde{\mathbf{\Sigma }}}^{-1}-\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left({\tilde{\mathbf{\Sigma }}}^{-1}\right)\}.\\ & =\frac{T}{2}\{2{\tilde{\mathbf{\Sigma }}}^{-1}(\mathit{I}-[\mathrm{E}(\mathbf{\Sigma }\mid \mathit{Y})+\mathrm{v}\mathrm{a}\mathrm{r}(\mathit{\mu }\mid \mathit{Y})\left]{\tilde{\mathbf{\Sigma }}}^{-1}\right)+\\ & -\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left({\tilde{\mathbf{\Sigma }}}^{-1}\right(\mathit{I}-\left[\mathrm{E}\right(\mathbf{\Sigma }\mid \mathit{Y}\left)+\mathrm{v}\mathrm{a}\mathrm{r}\right(\mathit{\mu }\mid \mathit{Y}\left)\right]{\tilde{\mathbf{\Sigma }}}^{-1})\}.\end{array}$ The derivative is 0 when $\begin{array}{r}\mathit{I}=\left[\mathrm{E}\right(\mathbf{\Sigma }\mid \mathit{Y})+\mathrm{v}\mathrm{a}\mathrm{r}(\mathit{\mu }\mid \mathit{Y}\left)\right]{\tilde{\mathbf{\Sigma }}}^{-1},\end{array}$ which yields ${\hat{\mathbf{\Sigma }}}_{\mathrm{i}\mathrm{i}\mathrm{d}}=\mathrm{E}(\mathbf{\Sigma }\mid \mathit{Y})+\mathrm{v}\mathrm{a}\mathrm{r}(\mathit{\mu }\mid \mathit{Y}).$

A.2. Posterior properties

A.2.1. $(\mathit{\theta }\mid {\mathit{X}}^{\ast },\mathit{X})$ under normal prior for $\mathbf{\Phi }$ and Jeffreys type prior for $\mathbf{\Sigma }$

A commonly used noninformative prior for $\mathbf{\Sigma }$ is the Jeffreys prior ${\pi }_J\left(\mathbf{\Sigma }\right)\propto |\mathbf{\Sigma }{|}^{-(p+1)/2}.$ The prior for $\mathbf{\Sigma }$ in the RATS statistical package is a modified version of the Jeffreys prior, ${\pi }_A\left(\mathbf{\Sigma }\right)\propto |\mathbf{\Sigma }{|}^{-(L+1)p/2-1}.$ Zellner's maximum data informative (MDI) prior, ${\pi }_M\left(\mathbf{\Sigma }\right)\propto |\mathbf{\Sigma }{|}^{-1/2}.$ For analysis on prior choice in VAR models see Ni and Sun (2003), and Sun and Ni (2004). We consider a class of joint priors, (A1) $\begin{array}{r}{\pi }_b(\mathit{\phi },\mathbf{\Sigma })={\pi }_N\left(\mathit{\phi }\right){\pi }_b^{\ast }\left(\mathbf{\Sigma }\right),\end{array}$(A1) where ${\pi }_N\left(\mathit{\phi }\right)$ is the normal prior for $\mathit{\phi }$ given by (10(10) $\begin{array}{r}{\pi }_N\left(\mathit{\phi }\right)\propto |{\mathit{M}}_0{|}^{-1/2}\exp \left\{-\frac{1}{2}(\mathit{\phi }-{\mathit{\phi }}_0{)}^{\prime }{\mathit{M}}_0^{-1}(\mathit{\phi }-{\mathit{\phi }}_0)\right\}.\end{array}$(10) ), and ${\pi }_b^{\ast }\left(\mathbf{\Sigma }\right)$ ($b\in \text{I R}$) is given by (A2) $\begin{array}{r}{\pi }_b^{\ast }\left(\mathbf{\Sigma }\right)\propto \frac{1}{|\mathbf{\Sigma }{|}^{b/2}}.\end{array}$(A2) Note that ${\pi }_{NJ},{\pi }_{NA}$, and ${\pi }_{NM}$ are special cases of (A2(A2) $\begin{array}{r}{\pi }_b^{\ast }\left(\mathbf{\Sigma }\right)\propto \frac{1}{|\mathbf{\Sigma }{|}^{b/2}}.\end{array}$(A2) ) when b equals to p + 1, $(L+1)p+2$ and 1, respectively.

We propose using the posterior quantities conditional on the simulated data ${\mathit{X}}^{\ast }$, $h(\mathit{\theta }\mid {\mathit{X}}^{\ast },\mathit{X})$, instead of the marginal posterior, $h(\mathit{\theta }\mid \mathit{X})$, as the estimator of $h\left(\mathit{\theta }\right)$. Note that the posterior $\begin{array}{rl}f(\mathbf{\Phi },\mathbf{\Sigma }\mid {\mathit{X}}^{\ast },\mathit{X})& \propto \frac{1}{|\mathbf{\Sigma }{|}^{T+b/2}}\\ & \times \mathrm{e}\mathrm{t}\mathrm{r}\left\{-\frac{1}{2}\left(\mathit{Y}-\mathit{X}\mathbf{\Phi }\right){\mathbf{\Sigma }}^{-1}(\mathit{Y}-\mathit{X}\mathbf{\Phi }{)}^{\prime }\right\}\\ & \times \mathrm{e}\mathrm{t}\mathrm{r}\left\{-\frac{1}{2}\left({\mathit{Y}}^{\ast }-{\mathit{X}}^{\ast }\mathbf{\Phi }\right){\mathbf{\Sigma }}^{-1}({\mathit{Y}}^{\ast }-{\mathit{X}}^{\ast }\mathbf{\Phi }{)}^{\prime }\right\}\\ & \times \exp \left\{-\frac{1}{2}(\mathit{\phi }-{\mathit{\phi }}_0{)}^{\prime }{\mathit{M}}_0^{-1}(\mathit{\phi }-{\mathit{\phi }}_0)\right\}.\end{array}$ We would like to express the posterior of $\mathbf{\Sigma }$ in terms of ${\mathit{X}}^{\ast }$. Integrating out $\mathbf{\Phi }$ results in $\begin{array}{rl}& f(\mathbf{\Sigma }\mid {\mathit{X}}^{\ast },\mathit{X})\\ & =\frac{|{\mathit{M}}_0+\mathbf{\Sigma }\otimes ({\mathit{X}}^{\prime }\mathit{X}+{{\mathit{X}}^{\ast }}^{\prime }{\mathit{X}}^{\ast }{)}^{-1}{|}^{1/2}}{|\mathbf{\Sigma }{|}^{T+b/2}}\\ & \times \mathrm{e}\mathrm{t}\mathrm{r}\left\{-\frac{1}{2}{\mathbf{\Sigma }}^{-1}\left[\right(\mathit{Y}-\mathit{X}{\hat{\mathbf{\Phi }}}_M{)}^{\prime }(\mathit{Y}-\mathit{X}{\hat{\mathbf{\Phi }}}_M)\right.\\ & +\left.({\mathit{Y}}^{\ast }-{\mathit{X}}^{\ast }{{\hat{\mathbf{\Phi }}}_M}^{\ast }{)}^{\prime }({\mathit{Y}}^{\ast }-{\mathit{X}}^{\ast }{{\hat{\mathbf{\Phi }}}_M}^{\ast }\left)\right]\right\}\\ & \times \exp \left\{-\frac{1}{2}(\tilde{\mathit{\phi }}-{\mathit{\phi }}_0{)}^{\prime }[{\mathit{M}}_0+\mathbf{\Sigma }\otimes ({\mathit{X}}^{\prime }\mathit{X}\right.\\ & +\left.{{\mathit{X}}^{\ast }}^{\prime }{\mathit{X}}^{\ast }{)}^{-1}{]}^{-1}(\tilde{\mathit{\phi }}-{\mathit{\phi }}_0)\right\}\\ & \times \exp \left\{-\frac{1}{2}({\hat{\mathit{\phi }}}_M-{{\hat{\mathit{\phi }}}_M}^{\ast }{)}^{\prime }{\mathbf{\Sigma }}^{-1}\otimes [({\mathit{X}}^{\prime }\mathit{X}{)}^{-1}\right.\\ & +\left.\left({{\mathit{X}}^{\ast }}^{\prime }{\mathit{X}}^{\ast }{)}^{-1}{]}^{-1}\right({\hat{\mathit{\phi }}}_M-{{\hat{\mathit{\phi }}}_M}^{\ast })\right\},\end{array}$ where the lower case defines the vec operator: $\mathit{\phi }=\mathrm{v}\mathrm{e}\mathrm{c}\left(\mathbf{\Phi }\right)$, $\tilde{\mathbf{\Phi }}=({\mathit{X}}^{\prime }\mathit{X}+{{\mathit{X}}^{\ast }}^{\prime }{\mathit{X}}^{\ast }{)}^{-1}({\mathit{X}}^{\prime }\mathit{Y}+{{\mathit{X}}^{\ast }}^{\prime }{\mathit{Y}}^{\ast }),$ ${\hat{\mathbf{\Phi }}}_M=({\mathit{X}}^{\prime }\mathit{X}{)}^{-1}{\mathit{X}}^{\prime }\mathit{Y}$, and ${{\hat{\mathbf{\Phi }}}_M}^{\ast }=({{\mathit{X}}^{\ast }}^{\prime }{\mathit{X}}^{\ast }{)}^{-1}{{\mathit{X}}^{\ast }}^{\prime }{\mathit{Y}}^{\ast }.$

The posterior $f(\mathbf{\Phi },\mathbf{\Sigma }\mid {\mathit{X}}^{\ast },\mathit{X})$ has a closed form when the prior for $\mathit{\phi }$ can be written as $N({\mathit{\phi }}_0,{\mathit{M}}_0)$ with (A3) $\begin{array}{r}{\mathit{M}}_0=\mathbf{\Sigma }\otimes {\mathbf{\Omega }}_0,\end{array}$(A3) where ${\mathbf{\Omega }}_0$ is a Lp + 1 by Lp + 1 known covariance matrix. In the extreme case of ${{\mathbf{\Omega }}_0}^{-1}\to 0$, ${\mathit{M}}_0^{-1}\to 0$ and the prior approaches a constant prior. Under the assumption (A3(A3) $\begin{array}{r}{\mathit{M}}_0=\mathbf{\Sigma }\otimes {\mathbf{\Omega }}_0,\end{array}$(A3) ), we have $\begin{array}{rl}& f(\mathbf{\Sigma }\mid {\mathit{X}}^{\ast },\mathit{X})\propto |{\mathbf{\Omega }}_0+({\mathit{X}}^{\prime }\mathit{X}+{{\mathit{X}}^{\ast }}^{\prime }{\mathit{X}}^{\ast }{)}^{-1}{|}^{p/2}\\ & |\mathbf{\Sigma }{|}^{-T-\frac{b-Lp-1}{2}}\mathrm{e}\mathrm{t}\mathrm{r}\left(-\frac{1}{2}\mathit{V}{\mathbf{\Sigma }}^{-1}\right),\end{array}$ where $\begin{array}{rl}& \mathit{V}(\mathit{X},{\mathit{X}}^{\ast })\\ & =({\hat{\mathbf{\Phi }}}_M-{{\hat{\mathbf{\Phi }}}_M}^{\ast }{)}^{\prime }[({\mathit{X}}^{\prime }\mathit{X}{)}^{-1}+({{\mathit{X}}^{\ast }}^{\prime }{\mathit{X}}^{\ast }{)}^{-1}\left]\right({\hat{\mathbf{\Phi }}}_M-{{\hat{\mathbf{\Phi }}}_M}^{\ast })\\ & +(\hat{\mathbf{\Phi }}-{\mathbf{\Phi }}_0{)}^{\prime }[{\mathbf{\Omega }}_0+({\mathit{X}}^{\prime }\mathit{X}+{{\mathit{X}}^{\ast }}^{\prime }{\mathit{X}}^{\ast }{)}^{-1}{]}^{-1}(\tilde{\mathbf{\Phi }}-{\mathbf{\Phi }}_0)\\ & +(\mathit{Y}-\mathit{X}{\hat{\mathbf{\Phi }}}_M{)}^{\prime }(\mathit{Y}-\mathit{X}{\hat{\mathbf{\Phi }}}_M)\\ & +({\mathit{Y}}^{\ast }-{\mathit{X}}^{\ast }{{\hat{\mathbf{\Phi }}}_M}^{\ast }{)}^{\prime }({\mathit{Y}}^{\ast }-{\mathit{X}}^{\ast }{{\hat{\mathbf{\Phi }}}_M}^{\ast }).\end{array}$$f(\mathbf{\Sigma }\mid {\mathit{X}}^{\ast },\mathit{X})\sim IW\{\mathit{V},2T+b-(L+1)p-2\}$. The mean of $\mathbf{\Sigma }\mid {\mathit{X}}^{\ast },\mathit{X}$ is $\mathit{V}(\mathit{X},{\mathit{X}}^{\ast })/(2T+b-(L+1)p-2-p-1)=\mathit{V}(\mathit{X},{\mathit{X}}^{\ast })/(2T+b-(L+2)p-3)$. The posterior mean $\mathbf{\Sigma }\mid \mathit{X}$ is estimated by $\sum _{j=1}^M(\mathbf{\Sigma }\mid {\mathit{X}}_j^{\ast },\mathit{X})/M$.

The marginal posterior of $f(\mathbf{\Phi }\mid {\mathit{X}}^{\ast },\mathit{X})$ can be obtained by integrating out $\mathbf{\Sigma }$ in $f(\mathbf{\Phi },\mathbf{\Sigma }\mid {\mathit{X}}^{\ast },\mathit{X})$. It is easy to verify that (A4) $\begin{array}{rl}& f(\mathbf{\Phi },\mathbf{\Sigma }\mid {\mathit{X}}^{\ast },\mathit{X})\propto |\mathbf{\Sigma }{|}^{-T-b/2}\mathrm{e}\mathrm{t}\mathrm{r}\left(-\frac{1}{2}\mathit{U}{\mathbf{\Sigma }}^{-1}\right)\\ & =|\mathbf{\Sigma }{|}^{-(\nu +p+1)/2}\mathrm{e}\mathrm{t}\mathrm{r}\left(-\frac{1}{2}\mathit{U}{\mathbf{\Sigma }}^{-1}\right),\end{array}$(A4) where the degree of freedom $\nu =2T+b-p-1$ and (A5) $\begin{array}{rl}\mathit{U}& =\left(\mathbf{\Phi }-{\mathbf{\Phi }}^u{)}^{\prime }\right({{\mathbf{\Omega }}_0}^{-1}+{\mathit{X}}^{\prime }\mathit{X}+{{\mathit{X}}^{\ast }}^{\prime }{\mathit{X}}^{\ast }\left)\right(\mathbf{\Phi }-{\mathbf{\Phi }}^u\left)+\mathit{V}\right(\mathit{X},{\mathit{X}}^{\ast }),\end{array}$(A5) (A6) $\begin{array}{rl}{\mathbf{\Phi }}^u& =\left({{\mathbf{\Omega }}_0}^{-1}+{\mathit{X}}^{\prime }\mathit{X}+{{\mathit{X}}^{\ast }}^{\prime }{\mathit{X}}^{\ast }{)}^{-1}\right({{\mathbf{\Omega }}_0}^{-1}{\mathbf{\Phi }}_0+{\mathit{X}}^{\prime }\mathit{Y}+{{\mathit{X}}^{\ast }}^{\prime }{\mathit{Y}}^{\ast }).\end{array}$(A6)

It follows from (A4(A4) $\begin{array}{rl}& f(\mathbf{\Phi },\mathbf{\Sigma }\mid {\mathit{X}}^{\ast },\mathit{X})\propto |\mathbf{\Sigma }{|}^{-T-b/2}\mathrm{e}\mathrm{t}\mathrm{r}\left(-\frac{1}{2}\mathit{U}{\mathbf{\Sigma }}^{-1}\right)\\ & =|\mathbf{\Sigma }{|}^{-(\nu +p+1)/2}\mathrm{e}\mathrm{t}\mathrm{r}\left(-\frac{1}{2}\mathit{U}{\mathbf{\Sigma }}^{-1}\right),\end{array}$(A4) ) that (A7) $\begin{array}{rl}& f(\mathbf{\Phi }\mid {\mathit{X}}^{\ast },\mathit{X})\propto |\mathit{U}{|}^{-\nu /2}\\ & =\left|(\mathbf{\Phi }-{\mathbf{\Phi }}^u{)}^{\prime }({{\mathbf{\Omega }}_0}^{-1}\right.\\ & {\left.+{\mathit{X}}^{\prime }\mathit{X}+{{\mathit{X}}^{\ast }}^{\prime }{\mathit{X}}^{\ast }\left)\right(\mathbf{\Phi }-{\mathbf{\Phi }}^u)+\mathit{V}(\mathit{X},{\mathit{X}}^{\ast })\right|}^{-\nu /2},\end{array}$(A7) which follows a matrix version of the Student-t distribution. The mean of $\mathbf{\Phi }\mid {\mathit{X}}^{\ast },\mathit{X}$ is ${\mathbf{\Phi }}^u$. To calculate the intrinsic Bayes estimator, note that the frequentist expectation $\hat{\mathit{W}}$ can be estimated by $\sum _{j=1}^M{\mathit{X}}_j^{\ast \mathrm{\prime }}{\mathit{X}}_j^{\ast }/M$. The posterior mean $\mathrm{E}(\mathit{W}\mathbf{\Phi }\mid \mathit{X})$ is estimated by $\sum _{j=1}^M{\mathit{X}}_j^{\ast \mathrm{\prime }}{\mathit{X}}_j^{\ast }\mathbf{\Phi }({\mathit{X}}_j^{\ast },\mathit{X})/M.$

A.2.2. Conditional posteriors under normal prior for $\mathbf{\Phi }$ and the Yang-Berger reference prior for $\mathbf{\Sigma }$

Fact A.1

Consider the normal prior for $\mathit{\phi }$ given in (10(10) $\begin{array}{r}{\pi }_N\left(\mathit{\phi }\right)\propto |{\mathit{M}}_0{|}^{-1/2}\exp \left\{-\frac{1}{2}(\mathit{\phi }-{\mathit{\phi }}_0{)}^{\prime }{\mathit{M}}_0^{-1}(\mathit{\phi }-{\mathit{\phi }}_0)\right\}.\end{array}$(10) ). The conditional density of $\mathit{\phi }$ given $(\mathbf{\Sigma };\mathit{Y})$ is $N_J({\mathit{\mu }}_M,{\mathit{V}}_M),$ where (A8) $\begin{array}{rl}{\mathit{\mu }}_M& ={\hat{\mathit{\phi }}}_{\mathrm{M}\mathrm{L}\mathrm{E}}+\{{\mathit{M}}_0^{-1}+{\mathbf{\Sigma }}^{-1}\otimes ({\mathit{X}}^{\prime }\mathit{X}\left){\}}^{-1}{\mathit{M}}_0^{-1}\right({\mathit{\phi }}_0-{\hat{\mathit{\phi }}}_{\mathrm{M}\mathrm{L}\mathrm{E}});\end{array}$(A8) (A9) $\begin{array}{rl}{\mathit{V}}_M& =\{{\mathit{M}}_0^{-1}+{\mathbf{\Sigma }}^{-1}\otimes ({\mathit{X}}^{\prime }\mathit{X}){\}}^{-1},\end{array}$(A9) where (A10) $\begin{array}{r}{\hat{\mathit{\phi }}}_{\mathrm{M}\mathrm{L}\mathrm{E}}=\mathrm{v}\mathrm{e}\mathrm{c}\left\{\right({\mathit{X}}^{\prime }\mathit{X}{)}^{-1}{\mathit{X}}^{\prime }\mathit{Y}\}.\end{array}$(A10)

Fact A.2

Consider the normal prior for $\mathit{\phi }$ given in (10(10) $\begin{array}{r}{\pi }_N\left(\mathit{\phi }\right)\propto |{\mathit{M}}_0{|}^{-1/2}\exp \left\{-\frac{1}{2}(\mathit{\phi }-{\mathit{\phi }}_0{)}^{\prime }{\mathit{M}}_0^{-1}(\mathit{\phi }-{\mathit{\phi }}_0)\right\}.\end{array}$(10) ). The conditional density of $\mathit{\phi }$ given $(\mathbf{\Sigma };\mathit{Y},{\mathit{Y}}^{\ast })$ is $\sim N({\mathit{\phi }}_M^{\ast },{\mathit{V}}_M^{\ast }),$ where $\begin{array}{rl}{\mathit{V}}_M^{\ast }& =\{{\mathit{M}}_0^{-1}+{\mathbf{\Sigma }}^{-1}\otimes ({\mathit{X}}^{\prime }\mathit{X}+{{\mathit{X}}^{\ast }}^{\prime }{\mathit{X}}^{\ast }){\}}^{-1},\\ {\mathit{\phi }}_M^{\ast }& ={\mathit{V}}_M^{\ast }\left[\right\{{\mathbf{\Sigma }}^{-1}\otimes ({\mathit{X}}^{\prime }\mathit{X}+{{\mathit{X}}^{\ast }}^{\prime }{\mathit{X}}^{\ast })\}{\hat{\mathit{\phi }}}_{\mathrm{M}\mathrm{L}\mathrm{E}}^{\ast }+{\mathit{M}}_0^{-1}{\mathit{\phi }}_0];\\ {\hat{\mathit{\phi }}}_{\mathrm{M}\mathrm{L}\mathrm{E}}^{\ast }& =\mathrm{v}\mathrm{e}\mathrm{c}\{({\mathit{X}}^{\prime }\mathit{X}+{{\mathit{X}}^{\ast }}^{\prime }{\mathit{X}}^{\ast }{)}^{-1}({\mathit{X}}^{\prime }\mathit{Y}+{{\mathit{X}}^{\ast }}^{\prime }{\mathit{Y}}^{\ast })\}.\end{array}$

Fact A.3

The conditional density of $\mathbf{\Sigma }$ given $(\mathit{\phi },\mathit{Y})$ is (A11) $\begin{array}{r}\pi (\mathbf{\Sigma }|\mathbf{\Phi },\mathit{Y})\propto \frac{\mathrm{e}\mathrm{t}\mathrm{r}\{-\frac{1}{2}{\mathbf{\Sigma }}^{-1}\mathit{S}(\mathbf{\Phi }\left)\right\}}{|\mathbf{\Sigma }{|}^{T/2+1}{\displaystyle \prod _{1\le i<j\le p}(d_i-d_j)}},\end{array}$(A11) where $\mathit{S}\left(\mathbf{\Phi }\right)=(\mathit{Y}-\mathit{X}\mathbf{\Phi }{)}^{\prime }(\mathit{Y}-\mathit{X}\mathbf{\Phi })$.

Fact A.4

The conditional density of $\mathbf{\Sigma }$ given $(\mathit{\phi },\mathit{Y},{\mathit{Y}}^{\ast })$ is (A12) $\begin{array}{rl}& \pi (\mathbf{\Sigma }|\mathbf{\Phi },\mathit{Y},{\mathit{Y}}^{\ast })\propto \frac{\mathrm{e}\mathrm{t}\mathrm{r}\{-\frac{1}{2}{\mathbf{\Sigma }}^{-1}(\mathit{S}\left(\mathbf{\Phi }\right)+{\mathit{S}}^{\ast }\left(\mathbf{\Phi }\right)\left)\right\}}{|\mathbf{\Sigma }{|}^{T+1}{\displaystyle \prod _{1\le i<j\le p}(d_i-d_j)}},\end{array}$(A12) where ${\mathit{S}}^{\ast }\left(\mathbf{\Phi }\right)=({\mathit{Y}}^{\ast }-{\mathit{X}}^{\ast }\mathbf{\Phi }{)}^{\prime }({\mathit{Y}}^{\ast }-{\mathit{X}}^{\ast }\mathbf{\Phi })$.

We have shown how simulated data facilitate computation of a frequentist moment in the Bayes estimator. In the appendix, we will show that the simulated data can also be used to reduced the variance of MCMC, making the simulation more efficient.

For simulation of $\mathbf{\Sigma }$, we adopt a hit-and-run algorithm used in Yang and Berger (1994). In implementing the algorithm, we consider a one-to-one transformation ${\mathbf{\Sigma }}^{\mathrm{\#}}=\log \left(\mathbf{\Sigma }\right)$, or $\mathbf{\Sigma }=\exp \left({\mathbf{\Sigma }}^{\mathrm{\#}}\right)$ in the sense that $\begin{array}{r}\mathbf{\Sigma }=\sum _{j=0}^{\mathrm{\infty }}\frac{({\mathbf{\Sigma }}^{\mathrm{\#}}{)}^j}{j!}.\end{array}$ The reason for simulating $\mathbf{\Sigma }$ as $\exp \left({\mathbf{\Sigma }}^{\mathrm{\#}}\right)$ is to ensure the generated $\mathbf{\Sigma }$ matrices are positive definite. It can be shown that the conditional posterior density of ${\mathbf{\Sigma }}^{\mathrm{\#}}$ given $(\mathbf{\Phi },\mathit{Y})$ is (A13) $\begin{array}{rl}& \pi ({\mathbf{\Sigma }}^{\mathrm{\#}}|\mathbf{\Phi },\mathit{Y})\\ & =\pi \left({\mathbf{\Sigma }}^{\mathrm{\#}}|\mathit{S}\right(\mathbf{\Phi }\left)\right)\propto \frac{\mathrm{e}\mathrm{t}\mathrm{r}\{-\frac{T}{2}|{\mathbf{\Sigma }}^{\mathrm{\#}}|-\frac{1}{2}(\mathrm{e}\mathrm{x}\mathrm{p}{\mathbf{\Sigma }}^{\ast }{)}^{-1}\mathit{S}\left(\mathbf{\Phi }\right)\}}{\prod _{i<j}(d_i^{\mathrm{\#}}-d_j^{\mathrm{\#}})},\end{array}$(A13) and that the conditional posterior density of ${\mathbf{\Sigma }}^{\ast }$ given $(\mathit{\phi },\mathit{Y},{\mathit{Y}}^{\ast })$ is then (A14) $\begin{array}{rl}& \pi ({\mathbf{\Sigma }}^{\mathrm{\#}}|\mathbf{\Phi },\mathit{Y},{\mathit{Y}}^{\ast })=\pi \left({\mathbf{\Sigma }}^{\mathrm{\#}}|\mathit{S}\right(\mathbf{\Phi }),{\mathit{S}}^{\ast }(\mathbf{\Phi }\left)\right)\\ & \propto \frac{\mathrm{e}\mathrm{t}\mathrm{r}[-T|{\mathbf{\Sigma }}^{\mathrm{\#}}|-\frac{1}{2}(\mathrm{e}\mathrm{x}\mathrm{p}{\mathbf{\Sigma }}^{\ast }{)}^{-1}\left\{\mathit{S}\right(\mathbf{\Phi })+{\mathit{S}}^{\ast }(\mathbf{\Phi }\left)\right\}]}{\prod _{i<j}(d_i^{\mathrm{\#}}-d_j^{\mathrm{\#}})},\end{array}$(A14) where ${\mathbf{\Sigma }}^{\mathrm{\#}}={\mathit{O}}^{\prime }{\mathit{D}}^{\mathrm{\#}}\mathit{O}$, $\mathit{O}$ is an orthogonal matrix, ${\mathit{D}}^{\mathrm{\#}}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(d_1^{\mathrm{\#}},\dots ,d_p^{\mathrm{\#}})$, with $d_1^{\mathrm{\#}}\ge \dots \ge d_p^{\mathrm{\#}}$. Note that $\mathrm{e}\mathrm{x}\mathrm{p}\left({\mathbf{\Sigma }}^{\mathrm{\#}}\right)={\mathit{O}}^{\prime }\exp \left({\mathit{D}}^{\mathrm{\#}}\right)\mathit{O}$.