# Derive Marginal Posterior to set up Gibbs-Sampler

I am currently trying to replicate a Hierarchical Model for multivariate returns proposed in the paper Portfolio selection using hierarchical Bayesian analysis and MCMC methods. However, in order to set up a Gibbs-Scheme I need the marginal posterior of the parameters. I did the calculations, however, my Gibbs sampler gives me some values I did not expect, so maybe someone can point me to possible mistakes in my calculations, especially for the scalar parameter $\mu_0$?

In the paper the authors propose a hierarchical model for the distribution of the return time series of the following form: \begin{alignat*}{2} \text{Likelihood: } & r|\mu,\Sigma &&\sim N\left(\mu,\Sigma\right)\\ \text{Prior: } &\mu|\mu_0,\Sigma&&\sim N\left(\mu_0\iota,\Sigma/\kappa_0\right)\\ &\Sigma &&\sim IW_{v}\left(\tau_0 ^{-2} P_0 ^{-1}\right)\\ \text{Hyperparameter: } &\mu_0 &&\sim \text{U}\left(-\infty,\infty\right)\\ &\tau_0 ^2 &&\sim IG\left(\varepsilon,\varepsilon\right) \end{alignat*} where $\varepsilon:=0.0001,\kappa_0:=0.1h$ and $P_0:=\left(\begin{array}{cccc} 1 & 0.5 &\ldots & 0.5 \\ 0.5 & 1 &\ldots & 0.5 \\ \vdots &\vdots &\ddots & \vdots \\ 0.5 & 0.5 &\ldots & 1 \end{array} \right).$ $h$ denotes the size of the estimation window.

The Likelihood function for a sample $D$ of observations $y_i$ $(i=1,\ldots,h)$ is defined as: \begin{align*} \mathcal{L}\left(D|\mu,\Sigma\right)\propto&\left|\Sigma\right|^{-\frac{h}{2}}\exp\left(-\frac{1}{2}\sum\limits_{i=1}^{h}\left(y_i-\mu\right)'\Sigma^{-1}\left(y_i-\mu\right)\right)\\ =&|\Sigma|^{-\frac{h}{2}}\exp\left(-\frac{1}{2}\sum\limits_{i=1}^{h}\left(y_i-\bar{y}+\bar{y}-\mu\right)'\Sigma^{-1}\left(y_i-\bar{y}+\bar{y}-\mu\right)\right)\\ =&|\Sigma|^{-\frac{h}{2}}\exp\left(-\frac{1}{2}\text{tr}\left(S\Sigma^{-1}\right)-\frac{h}{2}\left(\bar{y}-\mu\right)'\Sigma^{-1}\left(\bar{y}-\mu\right)\right) \end{align*} where $S=\sum\limits_{i=1}^{h}\left(y_i-\bar{y}\right)\left(y_i-\bar{y}\right)'$. The joint posterior has the following form: $$\pi(\mu,\Sigma,\mu_0,\tau_0 ^2|D)= \mathcal{L}(D|\mu,\Sigma)\pi\left(\mu|\mu_0,\frac{\Sigma}{\kappa_0}\right)\pi\left(\Sigma|\tau_0 ^2\right)\pi\left(\mu_0\right)\pi\left(\tau_0 ^2\right).$$ Calculation of $\pi(\mu|\Sigma,\mu_0,\tau_0 ^2,D)$: \begin{align*} \pi(\mu|\Sigma,\mu_0,\tau_0 ^2,D)= & \mathcal{L}(D|\mu,\Sigma)\pi(\mu|\mu_0,\frac{\Sigma}{\kappa_0})\\ \propto&\exp\left(-\frac{N}{2}(\bar{y}-\mu)'\Sigma^{-1}(\bar{y}-\mu)-\frac{\kappa_0}{2}(\mu-\iota\mu_0)'\Sigma^{-1}(\mu-\iota\mu_0)\right):=\exp(Z)\\ Z=&-\frac{1}{2}\left(N\text{ tr}((\mu-\bar{y})(\mu-\bar{y})'\Sigma^{-1})+\kappa_0 \text{ tr}((\mu-\iota\mu_0)(\mu-\iota\mu_0)'\Sigma^{-1})\right)\\ =&-\frac{1}{2}\text{ tr}\left((N(\mu\mu'-2\mu\bar{y}'+\bar{y}\bar{y}')+\kappa_0 (\mu\mu'-2\mu_0\mu \iota'+\mu_0^2 \iota \iota'))\Sigma^{-1}\right)\\ =&-\frac{1}{2}\text{ tr}\left((N+\kappa_0)\left[\mu\mu' -2\mu(\frac{n}{n+\kappa_0}\bar{y}'+\mu_0\frac{\kappa_0}{n+\kappa_0}\iota ')+\ldots\right]\Sigma^{-1}\right)\\ \Rightarrow \pi(\mu|\Sigma,\mu_0,\tau_0 ^2,D)\propto &\exp(-\frac{1}{2}(\mu-\mu_p)'(\frac{\Sigma}{n+\kappa_0})^{-1}(\mu-\mu_p))\\ \sim&N(\mu_p,\frac{\Sigma}{n+\kappa_0}) \end{align*} where $\mu_p:=\frac{n}{n+\kappa_0}\bar{y}+\mu_0\frac{\kappa_0}{n+\kappa_0}\iota.$

Calculation of $\pi(\Sigma|\mu,\mu_0,\tau_0 ^2,D)$: Plug-in leads to the density kernel: \begin{align*}\pi(\Sigma|\mu,\mu_0,\tau_0 ^2,D)&\propto \mathcal{L}(D|\mu,\Sigma)\pi(\mu|\mu_0,\frac{\Sigma}{\kappa_0})\pi(\Sigma|\tau_0 ^2)\\ &\propto|\Sigma|^{-\frac{v_p+k+1}{2}}\exp(-\frac{1}{2}\text{ tr}((\tau_0 ^{-2}P_0 ^{-1}+S+\kappa_0 (\mu-\mu_0\iota)(\mu-\mu_0\iota)'+ n(\bar{y}-\mu)(\bar{y}-\mu)')\Sigma^{-1}))\\ &\sim IW(v_p,\tau_0 ^{-2}P_0 ^{-1}+S+\kappa_0 (\mu-\mu_0\iota)(\mu-\mu_0\iota)'+n(\bar{y}-\mu)(\bar{y}-\mu)')\end{align*} where $v_p=v+n+1$. I used the linearity of the trace operator. Furthermore, $|\frac{\Sigma}{\kappa_0} |^{\frac{1}{2}}\propto |\Sigma| ^{\frac{1}{2}}$ as $\kappa_0$ is known.

Calculation of $\pi(\mu_0|\mu,\Sigma,\tau_0 ^2,D)$: \begin{align*} \pi(\mu_0|\mu,\Sigma,\tau_0 ^2,D)\propto& \pi(\mu|\mu_0,\frac{\Sigma}{\kappa_0})\pi(\mu_0)\\ \propto& \exp(-\frac{\kappa_0}{2}(\mu-\mu_0\iota)'\Sigma^{-1}(\mu-\mu_0\iota))\\ \propto&\exp\left(-\frac{\kappa_0}{2}(\mu_0^2(\iota'\Sigma^{-1}\iota)-2\mu_0(\mu'\Sigma^{-1}\iota))\right)\\ \propto&\exp\left(-\frac{\kappa_0}{2}\left(\mu_0\sqrt{\iota'\Sigma^{-1}\iota}-\frac{\mu'\Sigma^{-1}\iota}{\sqrt{\iota'\Sigma^{-1}\iota}}\right)^2\right)\\ \propto&\exp\left(-\frac{\kappa_0\iota'\Sigma^{-1}\iota}{2}\left(\mu_0-\frac{\mu'\Sigma^{-1}\iota}{(\iota'\Sigma^{-1}\iota)^\frac{3}{2}}\right)^2\right)\\ \sim&N\left(\frac{\mu'\Sigma^{-1}\iota}{(\iota'\Sigma^{-1}\iota)^\frac{3}{2}},\frac{1}{\kappa_0\iota'\Sigma^{-1}\iota}\right) \end{align*}

Calculation of $\pi(\tau_0 ^2|\mu,\Sigma,\mu_0,D)$ \begin{align*} \pi(\tau_0 ^2|\mu,\Sigma,\mu_0,D)\propto& \pi(\tau_0 ^2)\pi(\Sigma|\tau_0^2)\\ \propto&\tau_0^{-2\epsilon-2} \exp(-\frac{\epsilon}{\tau_0 ^2})|\tau_0 ^{-2}P_0 ^{-1}|^{\frac{v}{2}}\exp(-\frac{1}{2}\text{ tr}(\tau_0 ^{-2}P_0 ^{-1}\Sigma^{-1}))\\ \propto&\tau_0 ^{-2 (\epsilon+1+\frac{kv}{2})}\exp(-\frac{1}{\tau_0 ^2}(\frac{\text{ tr}(P_0 ^{-1}\Sigma^{-1})}{2}+\epsilon))\\ \sim& IG(\epsilon +\frac{kv}{2},\epsilon+\frac{\text{ tr}(P_0 ^{-1}\Sigma^{-1})}{2}) \end{align*} where I used the linearity of the trace-operator as well as $\det(rA)=r^k \det(A)$ for $r\in\mathbb{R}$, $A\in\mathbb{R}^{k x k}$.