# EM Estimation of Conditional Matrix Variate $t$ Distributions

Battulga Gankhuu111Department of Applied Mathematics, National University of
Mongolia; E-mail<EMAIL_ADDRESS>Phone Number: 976–99246036

###### Abstract

Conditional matrix variate student $t$ distribution was introduced by
Battulga24g. In this paper, we propose a new version of the conditional matrix
variate student $t$ distribution. The paper provides EM algorithms, which
estimate parameters of the conditional matrix variate student $t$
distributions, including general cases and special cases with Minnesota prior.

##  1 Introduction

An expectation–maximization (EM) algorithm was proposed and given its name by
Dempster77. The EM algorithm is an iterative method to obtain (local) maximum
likelihood estimate of parameters of distribution functions, which depend on
unobserved (latent) variables. The EM algorithm alternates an expectation (E)
step and a maximization (M) step. In E–step, one considers that conditional on
an available data and current estimate of the parameters, expectation of
augmented log–likelihood of the data and unobserved (latent) variables. The
E–Step defines a objective function. In M–step, to obtain parameter estimate
of next iteration, one maximizes the objective function with respect to the
parameters. Alternating between these steps, the EM algorithm produces
improved parameter estimates at each step (in the sense that the value of the
original log–likelihood is continually increased) and it converges to the
maximum likelihood (ML) estimates of the parameters.

The EM algorithm is widely used in econometrics. In particular, Hamilton90
introduced a parameter estimation method for general regime–switching model.
The regime–switching model assumes that a discrete unobservable Markov process
randomly switches among a finite set of regimes and that a particular
parameter set defines each regime. A Kalman filtering, which was introduced by
Kalman60 is an algorithm that provides estimates of some observed and
unobserved (state) processes. The Kalman filtering is an example of EM
algorithm. To estimate the parameters of the state–space model and to make
inferences about the state–space model, the Kalman filtering is used, see
Hamilton94 and Lutkepohl05. In finance, to value a private company, whose
market price is unobservable, Battulga24c applied the EM algorithm. McNeil05
provide EM algorithm to estimate parameters of generalized hyperbolic
distribution. Also, EM algorithm has been used in classifications. For
example, to estimate parameters of matrix variate $t$ distribution, Thompson20
used the EM algorithm. Sun10 provided EM algorithm to estimate parameters of
Pearson VII distribution.

Classic Vector Autoregressive (VAR) process was proposed by Sims (1980) who
criticize large–scale macro–econometric models, which are designed to model
interdependencies of economic variables. Besides Sims80, there are some other
important works on multiple time series modeling, see, e.g., Tiao81, where a
class of vector autoregressive moving average models was studied. For the VAR
process, a variable in the process is modeled by its past values and the past
values of other variables in the process. After the work of Sims80, VARs have
been used for macroeconomic forecasting and policy analysis. However, if the
number of variables in the system increases or the time lag is chosen high,
then too many parameters need to be estimated. This will reduce the degrees of
freedom of the model and entail a risk of over–parametrization.

Therefore, to reduce the number of parameters in a high–dimensional VAR
process, Litterman79 introduced probability distributions for coefficients
that are centered at the desired restrictions but that have a small and
nonzero variance. Those probability distributions are known as Minnesota prior
in Bayesian VAR (BVAR) literature, which is widely used in practice. Due to
over–parametrization, the generally accepted result is that the forecast of
the BVAR model is better than the VAR model estimated by the frequentist
technique. The BVAR relies on Monte–Carlo simulation. Recently, for Bayesian
Markov–Switching VAR process, Battulga24g introduced a new Monte–Carlo
simulation method that removes duplication in a regime vector. Also, the
author introduced importance sampling method to estimate probability of rare
event, which corresponds to endogenous variables. Research works have shown
that BVAR is an appropriate tool for modeling large data sets; for example,
see Banbura10.

The rest of the paper is organized as follows: In Section 2, for type I
conditional matrix variate distributions, including general case and special
case with Minnesota prior we develop EM algorithms. Section 3 is dedicated to
studying EM algorithms for type I conditional matrix variate distributions,
including general case and special case with Minnesota prior. Finally, Section
4 concludes the study.

##  2 EM Estimation of Type I Conditional Matrix Variate $t$ Distribution

We consider a Bayesian Vector Autoregressive process of $p$ order (BVAR($p$)),
which is given by the following equation

$y_{t}=A_{0,t}\psi_{t}+A_{1,t}y_{t-1}+\dots+A_{p,t}y_{t-p}+\xi_{t},~{}t=1,\dots,T,$
(1)

where $y_{t}=(y_{1,t},\dots,y_{n,t})^{\prime}$ is an $(n\times 1)$ vector of
endogenous variables, $\psi_{t}=(1,\psi_{2,t},\dots,\psi_{l,t})^{\prime}$ is
an $(l\times 1)$ vector of exogenous variables,
$\xi_{t}=(\xi_{1,t},\dots,\xi_{n,t})^{\prime}$ is an $(n\times 1)$ residual
process, $A_{0,t}$ is an $(n\times l)$ random coefficient matrix,
corresponding to the vector of exogenous variables, for $i=1,\dots,p$,
$A_{i,t}$ are $(n\times n)$ random coefficient matrices, corresponding to
$y_{t-1},\dots,y_{t-p}$. Equation (1) can be written by

$y_{t}=\Pi_{t}\mathsf{Y}_{t}+\xi_{t},~{}t=1,\dots,T,$ (2)

where $\Pi_{t}:=[A_{0,t}:A_{1,t}:\dots:A_{p,t}]$ is an $(n\times d)$ random
coefficient matrix with $d:=l+np$, which consist of all the random coefficient
matrices and
$\mathsf{Y}_{t}:=(\psi_{t}^{\prime},y_{t-1}^{\prime},\dots,y_{t-p}^{\prime})^{\prime}$
is a $(d\times 1)$ vector, which consist of exogenous variable $\psi_{t}$ and
last $p$ lagged values of the process $y_{t}$. The process $\mathsf{Y}_{t}$ is
measurable with respect to a $\sigma$–field $\mathcal{F}_{t-1}$, which is
defined below. Let us collect endogenous variables into ($[nT]\times 1$)
vector $y$, i.e., $y:=(y_{1}^{\prime},\dots,y_{T}^{\prime})^{\prime}$.

For the residual process $\xi_{t}$, we assume that it has
$\xi_{t}:=\Sigma_{t}^{1/2}\varepsilon_{t}$, $t=1,\dots,T$ representation,
where $\Sigma_{t}^{1/2}$ is a Cholesky factor of a positive definite $(n\times
n)$ random matrix $\Sigma_{t}$ and $\varepsilon_{1},\dots,\varepsilon_{T}$ is
a random sequence of independent identically multivariate normally distributed
random vectors with means of 0 and covariance matrices of $n$ dimensional
identity matrix $I_{n}$. We also assume that the strong white noise process
$\\{\varepsilon_{t}\\}_{t=1}^{T}$ is independent of the random coefficient
matrices $(\Pi_{1},\dots,\Pi_{T})$ and $(\Sigma_{1},\dots,\Sigma_{T})$
conditional on initial information
$\mathcal{F}_{0}:=\\{y_{1-p},\dots,y_{0},\psi_{1},\dots,\psi_{T}\\}$, where
$\psi_{1},\dots,\psi_{T}$ are values of exogenous variables and they are known
at time zero. We also denote available information at time $t$ by
$\mathcal{F}_{t}:=\\{\mathcal{F}_{0},y_{1},\dots,y_{t}\\}$.

###  2.1 General Type I Conditional Matrix Variate $t$ Distribution

Let us assume that for $t=1,\dots,T$, the random coefficient matrices
$\Pi_{t}$ and random covariance matrices $\Sigma_{t}$ are equals, that is,
$\Pi:=\Pi_{1}=\dots=\Pi_{T}$ and $\Sigma:=\Sigma_{1}=\dots=\Sigma_{T}$. Then,
by using the Kronecker product, the BVAR($p$) process can be written by the
following equation

$y_{t}=\Pi\mathsf{Y}_{t}+\xi_{t}=(\mathsf{Y}_{t}^{\prime}\otimes
I_{n})\pi+\xi_{t},~{}~{}~{}t=1,\dots,T,$ (3)

where $\otimes$ is the Kronecker product of two matrices and
$\pi:=\text{vec}(\Pi)$ is an $(nd\times 1)$ vectorization of the random
coefficient matrix $\Pi$. Now we define distributions of the random
coefficient vector $\pi$ and covariance matrix $\Sigma$. We assume that
conditional on the initial information $\mathcal{F}_{0}$, a distribution of
the random covariance matrix $\Sigma$ is given by

$\Sigma~{}|~{}\mathcal{F}_{0}\sim\mathcal{IW}(\nu_{0},V_{0})$ (4)

where the notation $\mathcal{IW}$ denotes the Inverse–Wishart distribution,
$\nu_{0}>n-1$ is a degrees of freedom and $V_{0}$ is a positive definite scale
matrix. Consequently, a distribution of the residual vector $\xi_{t}$ equals

$\xi_{t}~{}|~{}\Sigma,\mathcal{F}_{0}\sim\mathcal{N}\big{(}0,\Sigma\big{)},$
(5)

where $\mathcal{N}$ denotes the normal distribution. Also, we assume that
conditional on the covariance matrix $\Sigma$ and initial information
$\mathcal{F}_{0}$, a distribution of the random coefficient vector $\pi$ is
given by

$\pi~{}|~{}\Sigma,\mathcal{F}_{0}\sim\mathcal{N}\Big{(}\pi_{0},\Lambda_{0}\otimes\Sigma\Big{)},$
(6)

where $\pi_{0}$ is an $(nd\times 1)$ vector and
$\Lambda_{0}=\text{diag}\\{\lambda_{1},\dots,\lambda_{d}\\}$ is a diagonal
$(d\times d)$ matrix. Then, according to Battulga24g, the following
Proposition holds.

######  Proposition 1.

Let
$\pi~{}|~{}\Sigma,\mathcal{F}_{0}\sim\mathcal{N}\big{(}\pi_{0},\Lambda_{0}\otimes\Sigma\big{)}$,
and $\Sigma~{}|~{}\mathcal{F}_{0}\sim\mathcal{IW}(\nu_{0},V_{0})$. Then,
first, conditional on the initial information $\mathcal{F}_{0}$, a joint
density function of the random vector $\bar{y}_{t}$ is given by

$\displaystyle
f(y|\mathcal{F}_{0})=\frac{1}{\pi^{nT/2}}\frac{|\Lambda_{0}^{-1}|^{n/2}\Gamma_{n}\big{(}(\nu_{0}+T)/2\big{)}|V_{0}|^{\nu_{0}/2}}{|\Lambda_{0|T}|^{n/2}\Gamma_{n}(\nu_{0}/2)\big{|}B_{T}+V_{0}\big{|}^{(\nu_{0}+T)/2}},$
(7)

where $\Gamma_{n}(\cdot)$ is the multivariate gamma function,
$\Lambda_{0|T}^{-1}:=\mathsf{Y}^{\circ}(\mathsf{Y}^{\circ})^{\prime}+\Lambda_{0}^{-1}$
is a $(d\times d)$ matrix, $y^{\circ}:=[y_{1}:\dots:y_{T}]$ is an $(n\times
T)$ matrix, $\mathsf{Y}^{\circ}:=[\mathsf{Y}_{0}:\dots:\mathsf{Y}_{T-1}]$ is a
$(d\times T)$ matrix,
$\pi_{0}^{\circ}:=\mathbb{E}\big{(}\Pi\big{|}\Sigma,\mathcal{F}_{0}\big{)}$ is
an $(n\times d)$ matrix, and $B_{T}$ is an $(n\times n)$ positive
semi–definite matrix and equals

$\displaystyle
B_{T}:=\big{(}y^{\circ}-\pi_{0}^{\circ}\mathsf{Y}^{\circ}\big{)}\big{(}I_{T}+(\mathsf{Y}^{\circ})^{\prime}\Lambda_{0}\mathsf{Y}^{\circ}\big{)}^{-1}\big{(}y^{\circ}-\pi_{0}^{\circ}\mathsf{Y}^{\circ}\big{)}^{\prime}.$
(8)

Second, conditional on the random covariance matrix $\Sigma$ and information
$\mathcal{F}_{t}$, a joint density function of the random coefficient vector
$\pi$ is given by

$\displaystyle
f(\pi|\Sigma,\mathcal{F}_{T})=\frac{1}{(2\pi)^{nd/2}|\Lambda_{0|T}|^{n/2}|\Sigma|^{d/2}}\exp\bigg{\\{}-\frac{1}{2}\Big{(}\pi-\pi_{0|T}\Big{)}^{\prime}\big{(}\Lambda_{0|T}^{-1}\otimes\Sigma^{-1}\big{)}\Big{(}\pi-\pi_{0|T}\Big{)}\bigg{\\}},$
(9)

where $\pi_{0|T}:=\big{(}(\Lambda_{0|T}\mathsf{Y}^{\circ})\otimes
I_{n}\big{)}y+\big{(}(\Lambda_{0|T}\Lambda_{0}^{-1})\otimes
I_{n}\big{)}\pi_{0}$ is an $([nd]\times 1)$ vector. Finally, conditional on
the information $\mathcal{F}_{T}$, a joint density function of the random
coefficient matrix $\Sigma$ is given by

$\displaystyle f(\Sigma|\mathcal{F}_{T})$ $\displaystyle=$
$\displaystyle\frac{\big{|}B_{T}+V_{0}\big{|}^{(\nu_{0}+T)/2}}{\Gamma_{n}\big{(}(\nu_{0}+T)/2\big{)}2^{n(\nu_{0}+T)/2}}|\Sigma|^{-(\nu_{0}+T+n+1)/2}\exp\bigg{\\{}-\frac{1}{2}\mathrm{tr}\Big{(}\big{(}B_{T}+V_{0}\big{)}\Sigma^{-1}\Big{)}\bigg{\\}}.$
(10)

The joint density function (7) is called conditional matrix variate student
$t$ density, see Battulga24g. To differentiate from a conditional matrix $t$
distribution, which arise in the following section, we refer to the
distribution as type I conditional matrix variate student $t$ distribution.
From equations (9) and (10), one recognizes that posterior distributions of
the random coefficient vector $\pi$ and covariance matrix $\Sigma$ are
multivariate normal
$\mathcal{N}\big{(}\pi_{0|t},\Lambda_{0|T}\otimes\Sigma\big{)}$ and
inverse–Wishart $\mathcal{IW}\big{(}\nu_{0}+T,B_{T}+V_{0}\big{)}$,
respectively. Let us denote a vector of all the parameters of the joint
density function by $\theta:=\text{vec}(\pi_{0},\Lambda_{0},\nu_{0},V_{0})$.
In this section, we develop Expectation–Maximization (EM) algorithm to
estimate parameters of the density function. In E–Step, we consider that
conditional on the full information $\mathcal{F}_{T}$ and parameter at
iteration $k$, $\theta^{[k]}$, expectation of augmented log–likelihood of the
data $y$ and unobserved (latent) variables $\pi$ and $\Sigma$. The E–Step
defines a objective function $\mathcal{L}$, namely,

$\displaystyle\mathcal{L}$ $\displaystyle=$
$\displaystyle\mathbb{E}\bigg{[}-\frac{n(T+d)}{2}\ln(2\pi)+\frac{T+d+\nu_{0}+n+1}{2}\ln|\Sigma^{-1}|$
$\displaystyle-$
$\displaystyle\frac{1}{2}\sum_{t=1}^{T}\big{(}y_{t}-(\mathsf{Y}_{t}^{\prime}\otimes
I_{n})\pi\big{)}^{\prime}\Sigma^{-1}\big{(}y_{t}-(\mathsf{Y}_{t}^{\prime}\otimes
I_{n})\pi\big{)}$ $\displaystyle+$
$\displaystyle\frac{n}{2}\ln|\Lambda_{0}^{-1}|-\frac{1}{2}(\pi-\pi_{0})^{\prime}(\Lambda_{0}^{-1}\otimes\Sigma^{-1})(\pi-\pi_{0})$
$\displaystyle+$
$\displaystyle\frac{\nu_{0}}{2}\ln|V_{0}|-\ln\bigg{(}\Gamma_{n}\bigg{(}\frac{\nu_{0}}{2}\bigg{)}\bigg{)}-\frac{n\nu_{0}}{2}\ln(2)-\frac{1}{2}\mathrm{tr}\big{(}V_{0}\Sigma^{-1}\big{)}\bigg{|}\mathcal{F}_{T};\theta^{[k]}\bigg{]}$

In M–Step, to obtain parameter estimate of next iteration $\theta^{[k+1]}$,
one maximizes the objective function with respect to the parameter $\theta$.
First, let us consider partial derivative from the objective function with
respect to the parameter $\lambda_{i}$ for $i=1,\dots,d$. Since an inverse of
the matrix $\Lambda_{0}$ is
$\Lambda_{0}^{-1}=\text{diag}\\{1/\lambda_{1},\dots,1/\lambda_{d}\\}$, it is
clear that

$\displaystyle\frac{\partial(\pi-\pi_{0})^{\prime}(\Lambda_{0}^{-1}\otimes\Sigma^{-1})(\pi-\pi_{0})}{\partial\lambda_{i}}$
$\displaystyle=$
$\displaystyle-\frac{1}{\lambda_{i}^{2}}\mathrm{tr}\big{\\{}(\pi-\pi_{0})(\pi-\pi_{0})^{\prime}(\mathsf{E}_{ii}^{d}\otimes\Sigma^{-1})\big{\\}},$
(12)

where $\mathsf{E}_{ij}^{d}$ is a $(d\times d)$ matrix and its $(i,j)$–th
element equals 1 and others 0. In general case, to obtain parameter estimation
of the matrix $\Lambda_{0}$, one may use the following partial derivative

$\frac{\partial\mathrm{vec}\big{(}\Lambda_{0}^{-1}\otimes\Sigma^{-1}\big{)}}{\partial\mathrm{vec}(\Lambda_{0})^{\prime}}=-\Big{[}\mathrm{vec}\big{(}\Lambda_{0}^{-1}\mathsf{E}_{11}^{d}\Lambda_{0}^{-1}\otimes\Sigma^{-1}\big{)}:\dots:\mathrm{vec}\big{(}\Lambda_{0}^{-1}\mathsf{E}_{dd}^{d}\Lambda_{0}^{-1}\otimes\Sigma^{-1}\big{)}\Big{]}.$
(13)

It follows from equation (9) that

$\mathbb{E}\Big{[}(\pi-\pi_{0})(\pi-\pi_{0})^{\prime}\Big{|}\Sigma,\mathcal{F}_{T};\theta^{[k]}\Big{]}=\Big{(}\Lambda_{0|T}^{[k]}\otimes\Sigma\Big{)}+\Theta_{|T}^{[k]},$
(14)

where

$\Theta_{|T}^{[k]}:=\Big{(}\pi_{0|T}^{[k]}-\pi_{0}^{[k]}\Big{)}\Big{(}\pi_{0|T}^{[k]}-\pi_{0}^{[k]}\Big{)}^{\prime}$
(15)

is a positive semi–definite $(d\times d)$ matrix,

$\Lambda_{0|T}^{[k]}:=\Big{(}\mathsf{Y}^{\circ}(\mathsf{Y}^{\circ})^{\prime}+\big{(}\Lambda_{0}^{[k]}\big{)}^{-1}\Big{)}^{-1}$
(16)

is a positive semi–definite $(d\times d)$ matrix and

$\pi_{0|T}^{[k]}:=\Big{(}\Lambda_{0|T}^{[k]}\mathsf{Y}^{\circ}\otimes
I_{n}\Big{)}y+\Big{(}\Lambda_{0|T}^{[k]}\big{(}\Lambda_{0}^{[k]}\big{)}^{-1}\otimes
I_{n}\Big{)}\pi_{0}^{[k]}$ (17)

is an $([nd]\times 1)$ vector. $\Lambda_{0|T}^{[k]}$ and $\pi_{0|T}^{[k]}$ are
Bayesian estimators of the parameter matrix $\Lambda_{0}$ and the parameter
vector $\pi$ at iteration $k$, respectively. Note that the vector
$\pi_{0|T}^{[k]}$ does not depend on the random covariance matrix $\Sigma$.
According to the iterated expectation formula and expectation formula of
Wishart distributed random matrix $\Sigma^{-1}$, we have that

$\displaystyle\Psi_{i|T}^{[k]}$ $\displaystyle:=$
$\displaystyle\mathbb{E}\big{[}(\pi-\pi_{0})(\pi-\pi_{0})^{\prime}(\mathsf{E}_{ii}^{d}\otimes\Sigma^{-1})\big{|}\mathcal{F}_{T};\theta^{[k]}\big{]}$
(18) $\displaystyle=$
$\displaystyle\Big{(}\Lambda_{0|T}^{[k]}\mathsf{E}_{ii}^{d}\otimes
I_{n}\Big{)}+\Big{(}\nu_{0}^{[k]}+T\Big{)}\Theta_{|T}^{[k]}\Big{(}\mathsf{E}_{ii}^{d}\otimes\big{(}V_{0}^{[k]}+B_{T}^{[k]}\big{)}^{-1}\Big{)}.$

It can be shown that the matrix $\Psi_{i|T}^{[k]}$ is a positive semi–definite
matrix. Consequently, as $\ln|\Lambda_{0}|=\sum_{i=1}^{d}\ln(\lambda_{i})$,
for $i=1,\dots,d$, an estimator at iteration $(k+1)$ of the parameter
$\lambda_{i}$ is given by

$\displaystyle\lambda_{i}^{[k+1]}=\frac{1}{n}\text{tr}\Big{\\{}\Psi_{i|T}^{[k]}\Big{\\}}.$
(19)

Because $\Psi_{i|T}^{[k]}$ is a positive semi–definite matrix,
$\lambda_{i}^{[k+1]}$ takes a non–negative value. For $i=1,\dots,d$, we
collect $\lambda_{i}^{[k+1]}$ into a matrix
$\Lambda_{0}^{[k+1]}:=\text{diag}\Big{\\{}\lambda_{1}^{[k+1]},\dots,\lambda_{d}^{[k+1]}\Big{\\}}$.

Second, we consider partial derivative from the objective function
$\mathcal{L}$ with respect to the parameter vector $\pi_{0}$. The partial
derivative is given by

$\frac{\partial\mathcal{L}}{\partial\pi_{0}}=\mathbb{E}\big{[}(\Lambda_{0}^{-1}\otimes\Sigma^{-1})(\pi-\pi_{0})\big{|}\mathcal{F}_{T};\theta^{[k]}\big{]}.$
(20)

According to the iterated expectation formula and equations (9) and (10), one
obtains estimator at iteration $(k+1)$ of the parameter vector $\pi_{0}$

$\pi_{0}^{[k+1]}=\pi_{0|T}^{[k]}.$ (21)

Third, we consider partial derivative from the objective function
$\mathcal{L}$ with respect to the parameter matrix $V_{0}$. The partial
derivative is given by

$\frac{\partial\mathcal{L}}{\partial
V_{0}}=\mathbb{E}\bigg{[}\frac{\nu_{0}}{2}V_{0}^{-1}-\frac{1}{2}\Sigma^{-1}\bigg{|}\mathcal{F}_{T};\theta^{[k]}\bigg{]}.$
(22)

Consequently, we have that

$V_{0}^{[k+1]}=\frac{\nu_{0}}{\nu_{0}+T}\Big{(}V_{0}^{[k]}+B_{T}^{[k]}\Big{)}.$
(23)

Fourth, we consider partial derivative from the objective function
$\mathcal{L}$ with respect to the parameter $\nu_{0}$. The partial derivative
is given by

$\displaystyle\frac{\partial\mathcal{L}}{\partial\nu_{0}}=\mathbb{E}\bigg{[}\frac{1}{2}\ln|\Sigma^{-1}|+\frac{1}{2}\ln|V_{0}|-\frac{1}{2}\psi_{n}\bigg{(}\frac{\nu_{0}}{2}\bigg{)}-\frac{n}{2}\ln(2)\bigg{|}\mathcal{F}_{T};\theta^{[k]}\bigg{]},$
(24)

where $\psi_{n}(\cdot)$ is the multivariate digamma function, which is defined
by derivative of log of the multivariate gamma function. As a result, since
$\mathbb{E}\big{[}\ln|\Sigma^{-1}|\big{|}\mathcal{F}_{T}\big{]}=\psi_{n}(\frac{\nu_{0}+T}{2})+n\ln(2)-\ln|V_{0}+B_{T}|$,
see Nguyen23, estimator at iteration $(k+1)$ of the parameter $\nu_{0}$ is a
solution of the following equation

$\psi_{n}\bigg{(}\frac{\nu_{0}+T}{2}\bigg{)}-\psi_{n}\bigg{(}\frac{\nu_{0}}{2}\bigg{)}-\ln\Big{|}V_{0}^{[k]}+B_{T}^{[k]}\Big{|}+\ln\Big{|}V_{0}^{[k]}\Big{|}=0.$
(25)

In the following algorithm, we give EM algorithm for parameters of conditional
matrix variate $t$ density function

######  Algorithm 1.

(EM Estimation of General Type I Conditional Matrix Variate $t$ Distribution).

  * (1)

Set $k=0$ and initial value of the parameter vector
$\theta^{[0]}=\big{(}\pi_{0}^{[0]},\Lambda_{0}^{[0]},\nu_{0}^{[0]},V_{0}^{[0]}\big{)}$.

  * (2)

Calculate $\pi_{0}^{[k+1]}$ using equation (21).

  * (3)

For $i=1,\dots,d$, calculate $\lambda_{i}^{[k+1]}$ using equation (19) and
collect $\lambda_{i}^{[k+1]}$ into a diagonal matrix
$\Lambda_{0}^{[k+1]}:=\text{diag}\Big{\\{}\lambda_{1}^{[k+1]},\dots,\lambda_{d}^{[k+1]}\Big{\\}}.$

  * (4)

To obtain $\nu_{0}^{[k+1]}$, solve nonlinear equation (25) for $\nu_{0}$.

  * (5)

Calculate $V_{0}^{[k+1]}$ using equation (23).

  * (6)

Increase iteration count $k=k+1$ and go to step (2).

###  2.2 Type I Conditional Matrix Variate $t$ Distribution with Minnesota
Prior

In this subsection, we follow Battulga24g and we develop EM algorithm for
parameters of type I conditional matrix variate $t$ distribution with
Minnesota prior. In practice, one usually adopts the Minnesota prior to
estimating the parameters of the VAR$(p)$ process. The first version of
Minnesota prior was introduced by Litterman79. Also, Banbura10 used Minnesota
prior for large Bayesian VAR and showed that the forecast of large Bayesian
VAR is better than small Bayesian VAR. However, there are many different
variants of the Minnesota prior, we consider a prior, which is included in
Miranda18. The idea of Minnesota prior is that it shrinks diagonal elements of
the matrix $A_{1}$ toward $\phi_{i}$ and off–diagonal elements of $A_{1}$ and
all elements of other matrices $A_{0},A_{2},\dots,A_{p}$ toward 0, where
$\phi_{i}$ is 0 for a stationary variable $y_{i,t}$ and 1 for a variable with
unit root $y_{i,t}$. However, we adopt a different prior condition for the
random coefficient matrix $A_{0}$. For the prior, it is assumed that
conditional on $\Sigma$ and $\mathcal{F}_{0}$, $A_{0},A_{1},\dots,A_{p}$ are
jointly normally distributed, and for $(i,j)$–th element of the matrix
$A_{\ell}$ $(\ell=0,\dots,p)$, it holds that for $i=1,\dots,n$ and
$j=1,\dots,l$,

$\mathbb{E}\big{(}(A_{0})_{i,j}\big{|}\Sigma,\mathcal{F}_{0}\big{)}=(1-\phi_{i})C_{i,j},$
(26)

and

$\text{Var}\big{(}(A_{0})_{i,j}\big{|}\Sigma,\mathcal{F}_{0}\big{)}=(\sigma_{i}/\varepsilon_{j})^{2}$
(27)

and for $i,j=1,\dots,n$,

$\mathbb{E}\big{(}(A_{\ell})_{i,j}\big{|}\Sigma,\mathcal{F}_{0}\big{)}=\begin{cases}\phi_{i}&\text{if}~{}~{}~{}i=j,~{}\ell=1,\\\
0&\text{if}~{}~{}~{}\text{otherwise}\end{cases},~{}~{}~{}\text{for}~{}\ell=1,\dots,p$
(28)

and

$\text{Var}\big{(}(A_{\ell})_{i,j}\big{|}\Sigma,\mathcal{F}_{0}\big{)}=\begin{cases}\displaystyle\bigg{(}\frac{\sigma_{i}}{\alpha\ell^{\beta}\gamma_{i}}\bigg{)}^{2}&\text{if}~{}~{}~{}i=j,\\\
\displaystyle\bigg{(}\frac{\sigma_{i}}{\alpha\ell^{\beta}\gamma_{j}}\bigg{)}^{2}&\text{if}~{}~{}~{}\text{otherwise}\end{cases}~{}~{}~{}\text{for}~{}\ell=1,\dots,p,$
(29)

where $C_{i,j}$ is $(i,j)$–th element of an $(n\times l)$ parameter matrix $C$
and $\sigma_{i}^{2}$ is an $(i,i)$–th element of the random covariance matrix
$\Sigma$. It follows from equation (26) that the expectation of
$(A_{0})_{i,j}$ equals $C_{i,j}$ for stationary variable $y_{i,t}$ and zero
for non stationary variable with unit root $y_{i,t}$. A small
$\varepsilon_{i}^{2}$ corresponds to an uninformative diffuse prior for
$(A_{0})_{i,j}$, the parameter $\alpha$ controls the overall tightness of the
prior distribution, the parameter $\beta$ controls amount of information prior
information at higher lags, and $\gamma_{i}$ is a scaling parameter, see
Miranda18. Thus, the factor $1/\ell^{2\beta}$ represents a rate at which prior
variance decreases with increasing lag length.

According to Banbura10, it can be shown that the following equation satisfies
the prior conditions (26)–(29)

$\hat{y}^{\circ}=\Pi\hat{\mathsf{Y}}^{\circ}+\hat{\xi}^{\circ},$ (30)

where $\hat{y}^{\circ}$ and $\hat{\mathsf{Y}}^{\circ}$ are $(n\times d)$ and
$(d\times d)$ matrices of dummy variables and are defined by

$\hat{y}^{\circ}:=\big{[}\text{diag}\\{1-\phi_{1},\dots,1-\phi_{n}\\}C\text{diag}\\{\varepsilon_{1},\dots,\varepsilon_{l}\\}:\lambda_{1}\text{diag}\\{\phi_{1}\gamma_{1},\dots,\phi_{n}\gamma_{n}\\}:0_{[n\times
n(p-1)]}\big{]}$ (31)

and

$\hat{\mathsf{Y}}^{\circ}:=\begin{bmatrix}\text{diag}\\{\varepsilon_{1},\dots,\varepsilon_{l}\\}&0_{[l\times
np]}\\\ 0_{[np\times
l]}&\lambda_{1}\big{(}J_{\beta}\otimes\text{diag}\\{\gamma_{1},\dots,\gamma_{n}\\}\big{)}\end{bmatrix}$
(32)

with $J_{\beta}:=\text{diag}\\{1^{\beta},\dots,p^{\beta}\\}$, respectively,
and $\hat{\xi}^{\circ}:=[\xi_{1}:\dots:\xi_{d}]$ is an $(n\times d)$ matrix of
residual process. Note that one can add constraints for elements of the
coefficient matrix $\Pi$ to the matrices of dummy variables, see Miranda18. It
is worth mentioning that the matrices of dummy variables $\hat{y}_{t}$ and
$\hat{\mathsf{Y}}_{t}$ should not depend on the covariance matrix $\Sigma$. If
the dummy variables depend on the covariance matrix, an OLS estimator, and
matrix $\Lambda_{0}$ depend on the covariance matrix $\Sigma$, see below.
Consequently, in this case, one can not use the results of Proposition 1. For
this reason, we choose the prior condition (26)–(29). Equation (30) can be
written by

$\hat{y}=\big{(}(\hat{\mathsf{Y}}^{\circ})^{\prime}\otimes
I_{n}\big{)}\pi+\hat{\xi},$ (33)

where $\hat{y}$ and $\hat{\xi}$ are $([nd]\times 1)$ vectors and are
vectorizations of the matrix of dummy variables $\hat{y}^{\circ}$ and matrix
of the residual process $\hat{\xi}^{\circ}$, respectively, i.e.,
$\hat{y}:=\text{vec}(\hat{y}^{\circ})$ and
$\hat{\xi}:=\text{vec}(\hat{\xi}^{\circ})$. It follows from equation (33) that

$\pi\overset{d}{=}\big{(}((\hat{\mathsf{Y}}^{\circ}(\hat{\mathsf{Y}}^{\circ})^{\prime})^{-1}\hat{\mathsf{Y}}^{\circ})\otimes
I_{n}\big{)}\hat{y}+\big{(}((\hat{\mathsf{Y}}^{\circ}(\hat{\mathsf{Y}}^{\circ})^{\prime})^{-1}\hat{\mathsf{Y}}^{\circ})\otimes
I_{n}\big{)}\hat{\xi},$ (34)

where $d$ denotes equal distribution. It should be noted that the first term
of the right–hand side of the above equation is a vecorization of the ordinary
least square (OLS) estimator of the coefficient matrix $\Pi$, namely,

$\displaystyle\pi_{0}^{\circ}$ $\displaystyle:=$
$\displaystyle\hat{y}_{t}^{\circ}(\hat{\mathsf{Y}}^{\circ})^{\prime}(\hat{\mathsf{Y}}^{\circ}(\hat{\mathsf{Y}}^{\circ})^{\prime})^{-1}$
(35) $\displaystyle=$
$\displaystyle\big{[}\text{diag}\\{1-\phi_{1},\dots,1-\phi_{n}\\}C:\text{diag}\\{\phi_{1},\dots,\phi_{n}\\}:0_{[n\times
n(p-1)]}\big{]}.$

Note that the OLS estimator of the coefficient matrix $\Pi$ are same as the
prior conditions (26) and (28). Consequently, conditional on $\Sigma$ and
$\mathcal{F}_{0}$, a distribution of the coefficient vector $\pi$ is given by

$\pi~{}|~{}\Sigma,\mathcal{F}_{0}\sim\mathcal{N}\Big{(}\pi_{0},\big{(}\Lambda_{0}\otimes\Sigma\big{)}\Big{)}.$
(36)

where
$\Lambda_{0}:=(\hat{\mathsf{Y}}^{\circ}(\hat{\mathsf{Y}}^{\circ})^{\prime})^{-1}$
is a $(d\times d)$ diagonal matrix and its inverse equals

$\Lambda_{0}^{-1}:=\hat{\mathsf{Y}}^{\circ}(\hat{\mathsf{Y}}^{\circ})^{\prime}=\begin{bmatrix}\text{diag}\\{\varepsilon_{1}^{2},\dots,\varepsilon_{l}^{2}\\}&0_{[l\times
np]}\\\ 0_{[np\times
l]}&\text{diag}\\{1^{2\beta},\dots,p^{2\beta}\\}\otimes\text{diag}\\{\alpha^{2}\gamma_{1}^{2},\dots,\alpha^{2}\gamma_{n}^{2}\\}\end{bmatrix}.$
(37)

Since the matrix $\Lambda_{0}^{-1}$ is a diagonal matrix, its determinant,
which appears in the objective function (2.1) is

$\ln|\Lambda_{0}^{-1}|=2\sum_{j=1}^{l}\ln(\varepsilon_{j})+2np\ln(\alpha)+\beta\sum_{\ell=1}^{p}\ln(\ell)+2p\sum_{i=1}^{n}\ln(\gamma_{i}).$
(38)

Let $c:=\text{vec}(C)$ be a $(nl\times 1)$ vectorization of the parameter
matrix $C$. Then, partial derivative from the objective function $\mathcal{L}$
with respect to the parameter $c$ is

$\displaystyle\frac{\partial\mathcal{L}}{\partial c^{\prime}}$
$\displaystyle=$
$\displaystyle\mathbb{E}\big{[}(\pi-\pi_{0})^{\prime}(\Lambda_{0}^{-1}\otimes\Sigma^{-1})J_{\phi}^{\prime}\big{|}\mathcal{F}_{T};\theta^{[k]}\big{]}$
(39) $\displaystyle=$
$\displaystyle\varepsilon^{2}\mathbb{E}\big{[}(\pi_{c}-c)^{\prime}\big{(}I_{l}\otimes\Sigma^{-1}\big{)}\big{|}\mathcal{F}_{T};\theta^{[k]}\big{]},$

where
$J_{\phi}:=\big{[}I_{l}\otimes\text{diag}\\{1-\phi_{1},\dots,1-\phi_{n}\\}:0_{[ln\times
n^{2}p]}\big{]}$ is an $(ln\times dn)$ matrix and $\pi_{c}:=J_{\phi}\pi$ is an
$(ln\times 1)$ vector. Then, according to the iterated expectation formula, we
get an estimator at iteration $(k+1)$ of the parameter vector $c$

$c^{[k+1]}:=J_{\phi}\pi_{0|T}^{[k]}.$ (40)

where $\pi_{0|T}^{[k]}$ is given by equation (17). To obtain estimators at
iteration $(k+1)$ of the parameters $\alpha$, $\beta$, $\gamma_{i}$, and
$\varepsilon_{j}$ for $i=1,\dots,n$ and $j=1,\dots,l$, we define the following
$(d\times d)$ matrices:

$\Delta_{\varepsilon_{j}}:=\begin{bmatrix}\mathsf{E}_{jj}^{l}&0_{[l\times
np]}\\\ 0_{[np\times l]}&0_{[np\times np]}\end{bmatrix},$ (41)
$\Delta_{\alpha}^{[k]}:=\begin{bmatrix}0_{l\times l}&0_{[l\times np]}\\\
0_{[np\times
l]}&\text{diag}\big{\\{}1^{2\beta^{[k]}},\dots,p^{2\beta^{[k]}}\big{\\}}\otimes\text{diag}\big{\\{}\big{(}\gamma_{1}^{[k]}\big{)}^{2},\dots,\big{(}\gamma_{n}^{[k]}\big{)}^{2}\big{\\}}\end{bmatrix},$
(42)
$\Delta_{\beta}^{[k]}:=\big{(}\alpha^{[k]}\big{)}^{2}\begin{bmatrix}0_{l\times
l}&0_{[l\times np]}\\\ 0_{[np\times
l]}&\text{diag}\big{\\{}1^{2\beta}\ln(1),\dots,p^{2\beta}\ln(p)\big{\\}}\otimes\text{diag}\big{\\{}\big{(}\gamma_{1}^{[k]}\big{)}^{2},\dots,\big{(}\gamma_{n}^{[k]}\big{)}^{2}\big{\\}}\end{bmatrix},$
(43)

and

$\Delta_{\gamma_{i}}^{[k]}:=\big{(}\alpha^{[k]}\big{)}^{2}\begin{bmatrix}0_{l\times
l}&0_{[l\times np]}\\\ 0_{[np\times
l]}&\text{diag}\big{\\{}1^{2\beta^{[k]}}\ln(1),\dots,p^{2\beta^{[k]}}\ln(p)\big{\\}}\otimes\mathsf{E}_{ii}^{n}\end{bmatrix}.$
(44)

It follows from determinant equation (38) and the objective function (2.1)
that similarly to equation (19), one obtains that for $i=1,\dots,n$ and
$j=1,\dots,l$, estimators at iteration $(k+1)$ of the parameters
$\varepsilon_{j}$, $\alpha$, and $\gamma_{i}$ are given by

$\displaystyle\Big{(}\varepsilon_{j}^{[k+1]}\Big{)}^{2}:=\frac{n}{\text{tr}\Big{\\{}\Big{(}\Lambda_{0|T}^{[k]}\Delta_{\varepsilon_{j}}\otimes
I_{n}\Big{)}+\Big{(}\nu_{0}^{[k]}+T\Big{)}\Theta_{|T}^{[k]}\Big{(}\Delta_{\varepsilon_{j}}\otimes\big{(}V_{0}^{[k]}+B_{T}^{[k]}\big{)}^{-1}\Big{)}\Big{\\}}},$
(45)
$\displaystyle\Big{(}\alpha^{[k+1]}\Big{)}^{2}:=\frac{n^{2}p}{\text{tr}\Big{\\{}\Big{(}\Lambda_{0|T}^{[k]}\Delta_{\alpha}^{[k]}\otimes
I_{n}\Big{)}+\Big{(}\nu_{0}^{[k]}+T\Big{)}\Theta_{|T}^{[k]}\Big{(}\Delta_{\alpha}^{[k]}\otimes\big{(}V_{0}^{[k]}+B_{T}^{[k]}\big{)}^{-1}\Big{)}\Big{\\}}},$
(46)

and

$\displaystyle\Big{(}\gamma_{i}^{[k+1]}\Big{)}^{2}:=\frac{np}{\text{tr}\Big{\\{}\Big{(}\Lambda_{0|T}^{[k]}\Delta_{\gamma_{i}}^{[k]}\otimes
I_{n}\Big{)}+\Big{(}\nu_{0}^{[k]}+T\Big{)}\Theta_{|T}^{[k]}\Big{(}\Delta_{\gamma_{i}}^{[k]}\otimes\big{(}V_{0}^{[k]}+B_{T}^{[k]}\big{)}^{-1}\Big{)}\Big{\\}}}.$
(47)

For $p\geq 2$, an estimator at iteration $(k+1)$ of the parameter $\beta$,
$\beta^{[k+1]}$ is obtained from the following nonlinear equation

$\displaystyle\sum_{\ell=1}^{p}\ln(\ell)=\frac{1}{n}\text{tr}\Big{\\{}\Big{(}\Lambda_{0|T}^{[k]}\Delta_{\beta}^{[k]}\otimes
I_{n}\Big{)}+\Big{(}\nu_{0}^{[k]}+T\Big{)}\Theta_{|T}^{[k]}\Big{(}\Delta_{\beta}^{[k]}\otimes\big{(}V_{0}^{[k]}+B_{T}^{[k]}\big{)}^{-1}\Big{)}\Big{\\}}.$
(48)

It should be noted that it is not difficult to show that the right–hand sides
of equations (45)–(48) take non negative values. Consequently, an EM algorithm
for parameters of conditional matrix variate $t$ distribution with Minnesota
prior is given by the following algorithm.

######  Algorithm 2.

(EM Estimation of Type I Conditional Matrix Variate $t$ Distribution with
Minnesota Prior).

  * (1)

Set $k=0$ and initial value of the parameter vector

$\theta^{[0]}=\Big{(}C^{[0]},\varepsilon_{1}^{[0]},\dots,\varepsilon_{l}^{[0]},\alpha^{[0]},\beta^{[0]},\gamma_{1}^{[0]},\dots,\gamma_{n}^{[0]},\nu_{0}^{[0]},V_{0}^{[0]}\Big{)}.$
(49)

  * (2)

Calculate $C^{[k+1]}$ using equation (40).

  * (3)

For $i=1,\dots,n$ and $j=1,\dots,l$, calculate
$\Big{(}\varepsilon_{j}^{[k+1]}\Big{)}^{2}$,
$\Big{(}\alpha^{[k+1]}\Big{)}^{2}$, and $\Big{(}\gamma_{i}^{[k+1]}\Big{)}^{2}$
using equations (45)–(47).

  * (4)

For $p\geq 2$, to obtain $\beta^{[k+1]}$, solve nonlinear equation (48) for
$\beta$.

  * (5)

To obtain $\nu_{0}^{[k+1]}$, solve nonlinear equation (25) for $\nu_{0}$.

  * (6)

Calculate $V_{0}^{[k+1]}$ using equation (23).

  * (7)

Increase iteration count $k=k+1$ and go to step (2).

##  3 EM Estimation of Type II Conditional Matrix Variate $t$ Distribution

Let us reconsider a BVAR($p$) process, given in equation (1), namely,

$y_{t}=\Pi_{t}\mathsf{Y}_{t}+\xi_{t}=(\mathsf{Y}_{t}^{\prime}\otimes
I_{n})\pi_{t}+\xi_{t},~{}~{}~{}t=1,\dots,T,$ (50)

where $\pi_{t}:=\text{vec}(\Pi_{t})$ is an $(nd\times 1)$ vectorization of the
random coefficient matrix $\Pi_{t}$. For the random coefficient vectors and
covariance matrices, we assume the following assumption holds.

######  Assumption 1.

Conditional on the initial information $\mathcal{F}_{0}$, the random
coefficient vectors and covariance matrices
$(\pi_{1},\Sigma_{1}),\dots,(\pi_{T},\Sigma_{T})$ are independent and
identically distributed.

Note that from the assumption, we can conclude that conditional on initial
information $\mathcal{F}_{0}$, for each $t=2,\dots,T$, $(\pi_{t},\Sigma_{t})$
is independent of a random vector
$(y_{1}^{\prime},\dots,y_{t-1}^{\prime})^{\prime}$. Now we define
distributions of the random coefficient vector $\pi_{t}$ and covariance matrix
$\Sigma_{t}$. We suppose that conditional on the initial information
$\mathcal{F}_{0}$, a distribution of the random covariance matrix $\Sigma_{t}$
is given by

$\Sigma_{t}~{}|~{}\mathcal{F}_{0}\sim\mathcal{IW}(\nu_{0},V_{0}).$ (51)

where $\nu_{0}>n-1$ is a degrees of freedom and $V_{0}$ is a positive definite
scale matrix. Hence, a distribution of the residual vector $\xi_{t}$ equals

$\xi_{t}~{}|~{}\Sigma_{t},\mathcal{F}_{0}\sim\mathcal{N}\big{(}0,\Sigma_{t}\big{)}.$
(52)

Also, we suppose that conditional on the covariance matrix $\Sigma_{t}$ and
initial information $\mathcal{F}_{0}$, a distribution of the random
coefficient vector $\pi_{t}$ is given by

$\pi_{t}~{}|~{}\Sigma_{t},\mathcal{F}_{0}\sim\mathcal{N}\Big{(}\pi_{0},\Lambda_{0}\otimes\Sigma_{t}\Big{)},$
(53)

where $\pi_{0}$ is an $(nd\times 1)$ vector and $\Lambda_{0}$ is a symmetric
positive definite $(d\times d)$ matrix. Then, the following Proposition holds.

######  Proposition 2.

Let for $t=1,\dots,T$,
$\pi_{t}~{}|~{}\Sigma_{t},\mathcal{F}_{0}\sim\mathcal{N}\big{(}\pi_{0},\Lambda_{0}\otimes\Sigma_{t}\big{)}$,
$\Sigma_{t}~{}|~{}\mathcal{F}_{0}\sim\mathcal{IW}(\nu_{0},V_{0})$, and
Assumption 1 holds. Then, first, conditional on the information
$\mathcal{F}_{t-1}$, a joint density function of the random vector $y_{t}$ is
given by

$\displaystyle
f(y_{t}|\mathcal{F}_{t-1})=\frac{1}{\pi^{n/2}}\frac{\Gamma_{n}\big{(}(\nu_{0}+1)/2\big{)}|V_{0}|^{\nu_{0}/2}}{\big{(}1+\mathsf{Y}_{t}^{\prime}\Lambda_{0}\mathsf{Y}_{t}\big{)}^{n/2}\Gamma_{n}(\nu_{0}/2)\big{|}\tilde{B}_{t}+V_{0}\big{|}^{(\nu_{0}+1)/2}},$
(54)

where

$\displaystyle\tilde{B}_{t}:=\frac{1}{1+\mathsf{Y}_{t}^{\prime}\Lambda_{0}\mathsf{Y}_{t}}\big{(}y_{t}-\pi_{0}^{\circ}\mathsf{Y}_{t}\big{)}\big{(}y_{t}-\pi_{0}^{\circ}\mathsf{Y}_{t}\big{)}^{\prime}$
(55)

is a positive semi–define $(n\times n)$ matrix. Second, conditional on the
random covariance matrix $\Sigma_{t}$ and information $\mathcal{F}_{t}$, a
joint density function of the random coefficient vector $\pi_{t}$ is given by

$\displaystyle
f(\pi_{t}|\Sigma_{t},\mathcal{F}_{t})=\frac{1}{(2\pi)^{nd/2}|\tilde{\Lambda}_{0|t}|^{n/2}|\Sigma|^{d/2}}\exp\bigg{\\{}-\frac{1}{2}\Big{(}\pi_{t}-\tilde{\pi}_{0|t}\Big{)}^{\prime}\big{(}\tilde{\Lambda}_{0|t}^{-1}\otimes\Sigma_{t}^{-1}\big{)}\Big{(}\pi_{t}-\tilde{\pi}_{0|t}\Big{)}\bigg{\\}},$
(56)

where
$\tilde{\Lambda}_{0|t}^{-1}:=\mathsf{Y}_{t}\mathsf{Y}_{t}+\Lambda_{0}^{-1}$ is
a $(d\times d)$ matrix and
$\tilde{\pi}_{0|t}:=\big{(}(\tilde{\Lambda}_{0|T}\mathsf{Y})\otimes
I_{n}\big{)}y_{t}+\big{(}(\tilde{\Lambda}_{0|t}\Lambda_{0}^{-1})\otimes
I_{n}\big{)}\pi_{0}$ is an $([nd]\times 1)$ vector. Finally, conditional on
the information $\mathcal{F}_{t}$, a joint density function of the random
coefficient matrix $\Sigma$ is given by

$\displaystyle f(\Sigma_{t}|\mathcal{F}_{t})$ $\displaystyle=$
$\displaystyle\frac{\big{|}\tilde{B}_{t}+V_{0}\big{|}^{(\nu_{0}+1)/2}}{\Gamma_{n}\big{(}(\nu_{0}+1)/2\big{)}2^{n(\nu_{0}+1)/2}}|\Sigma_{t}|^{-(\nu_{0}+n+2)/2}\exp\bigg{\\{}-\frac{1}{2}\mathrm{tr}\Big{(}\big{(}\tilde{B}_{t}+V_{0}\big{)}\Sigma_{t}^{-1}\Big{)}\bigg{\\}}.$
(57)

We refer to distribution function, which is given in equation (54) as type II
conditional matrix variate $t$ distribution function. In this section, we
introduce an EM algorithm to estimate the parameters of the distribution
function. Similarly to equation (2.1), an conditional expectation (objective
function) in E–Step is given by

$\displaystyle\mathcal{L}$ $\displaystyle=$
$\displaystyle\mathbb{E}\bigg{[}-\frac{nT(1+d)}{2}\ln(2\pi)+\sum_{t=1}^{T}\frac{d+\nu_{0}+n+2}{2}\ln|\Sigma_{t}^{-1}|$
$\displaystyle-$
$\displaystyle\frac{1}{2}\sum_{t=1}^{T}\big{(}y_{t}-(\mathsf{Y}_{t}^{\prime}\otimes
I_{n})\pi_{t}\big{)}^{\prime}\Sigma_{t}^{-1}\big{(}y_{t}-(\mathsf{Y}_{t}^{\prime}\otimes
I_{n})\pi_{t}\big{)}$ $\displaystyle-$
$\displaystyle\frac{nT}{2}\ln|\Lambda_{0}|-\frac{1}{2}\sum_{t=1}^{T}(\pi_{t}-\pi_{0})^{\prime}(\Lambda_{0}^{-1}\otimes\Sigma_{t}^{-1})(\pi_{t}-\pi_{0})$
$\displaystyle+$
$\displaystyle\frac{\nu_{0}T}{2}\ln|V_{0}|-\sum_{t=1}^{T}\ln\bigg{(}\Gamma_{n}\bigg{(}\frac{\nu_{0}}{2}\bigg{)}\bigg{)}-\sum_{t=1}^{T}\frac{n\nu_{0}}{2}\ln(2)-\frac{1}{2}\sum_{t=1}^{T}\mathrm{tr}\big{(}V_{0}\Sigma_{t}^{-1}\big{)}\bigg{|}\mathcal{F}_{T};\theta^{[k]}\bigg{]}.$

In M–Step, to obtain ML estimators at iteration $(k+1)$ of the parameters of
the type II conditional matrix variate $t$ distribution, we need the following
Lemma.

######  Lemma 1.

Under Assumption 1, conditional on the full information $\mathcal{F}_{T}$, a
joint density function of coefficient vector $\pi_{t}$ and covariance matrix
$\Sigma_{t}$ is given by

$f(\pi_{t},\Sigma_{t}|\mathcal{F}_{T})=f(\pi_{t},\Sigma_{t}|\mathcal{F}_{t}).$
(59)

By using Proposition 2, Lemma 1, and ideas in subsection 2.1, one can arrive
the following EM Algorithm, which estimates parameters of the general type II
conditional matrix variate $t$ distribution.

######  Algorithm 3.

(EM Estimation of General Type II Conditional Matrix Variate $t$
Distribution).

  * (1)

Set $k=0$ and initial value of the parameter vector
$\theta^{[0]}=\big{(}\pi_{0}^{[0]},\Lambda_{0}^{[0]},\nu_{0}^{[0]},V_{0}^{[0]}\big{)}$.

  * (2)

Calculate $\pi_{0}^{[k+1]}$ using equation

$\displaystyle\pi_{0}^{[k+1]}$ $\displaystyle=$
$\displaystyle\bigg{(}\sum_{t=1}^{T}\Big{(}I_{d}\otimes\big{(}V_{0}^{[k]}+\tilde{B}_{t}^{[k]}\big{)}^{-1}\Big{)}\bigg{)}^{-1}\sum_{t=1}^{T}\Big{(}I_{d}\otimes\big{(}V_{0}^{[k]}+\tilde{B}_{t}^{[k]}\big{)}^{-1}\Big{)}\tilde{\pi}_{0|t}^{[k]},$
(60)

where the matrix and vector are given by

$\displaystyle\tilde{B}_{t}^{[k]}:=\frac{1}{1+\mathsf{Y}_{t}^{\prime}\Lambda_{0}^{[k]}\mathsf{Y}_{t}}\Big{(}y_{t}-\big{(}\pi_{0}^{\circ}\big{)}^{[k]}\mathsf{Y}_{t}\Big{)}\Big{(}y_{t}-\big{(}\pi_{0}^{\circ}\big{)}^{[k]}\mathsf{Y}_{t}\Big{)}^{\prime}$
(61)

and

$\tilde{\pi}_{0|t}^{[k]}:=\Big{(}\tilde{\Lambda}_{0|t}^{[k]}\mathsf{Y}_{t}\otimes
I_{n}\Big{)}y_{t}+\Big{(}\tilde{\Lambda}_{0|t}^{[k]}\big{(}\Lambda_{0}^{[k]}\big{)}^{-1}\otimes
I_{n}\Big{)}\pi_{0}^{[k]}.$ (62)

  * (3)

For $i=1,\dots,d$, calculate $\lambda_{i}^{[k+1]}$ using equation

$\displaystyle\lambda_{i}^{[k+1]}=\frac{1}{nT}\sum_{t=1}^{T}\text{tr}\Big{\\{}\Big{(}\tilde{\Lambda}_{0|t}^{[k]}\mathsf{E}_{ii}^{d}\otimes
I_{n}\Big{)}+\Big{(}\nu_{0}^{[k]}+1\Big{)}\tilde{\Theta}_{|t}^{[k]}\Big{(}\mathsf{E}_{ii}^{d}\otimes\big{(}V_{0}^{[k]}+\tilde{B}_{t}^{[k]}\big{)}^{-1}\Big{)}\Big{\\}},$
(63)

where the matrices are given by

$\tilde{\Theta}_{|t}^{[k]}:=\Big{(}\tilde{\pi}_{0|t}^{[k]}-\pi_{0}^{[k]}\Big{)}\Big{(}\tilde{\pi}_{0|t}^{[k]}-\pi_{0}^{[k]}\Big{)}^{\prime}~{}~{}~{}\text{and}~{}~{}~{}\tilde{\Lambda}_{0|t}^{[k]}:=\Big{(}\mathsf{Y}_{t}\mathsf{Y}_{t}^{\prime}+\big{(}\Lambda_{0}^{[k]}\big{)}^{-1}\Big{)}^{-1}.$
(64)

Collect $\lambda_{i}^{[k+1]}$ into a diagonal matrix
$\Lambda_{0}^{[k+1]}:=\text{diag}\Big{\\{}\lambda_{1}^{[k+1]},\dots,\lambda_{d}^{[k+1]}\Big{\\}}.$

  * (4)

To obtain $\nu_{0}^{[k+1]}$, solve the following nonlinear equation for
$\nu_{0}$

$\sum_{t=1}^{T}\Bigg{(}\psi_{n}\bigg{(}\frac{\nu_{0}+1}{2}\bigg{)}-\psi_{n}\bigg{(}\frac{\nu_{0}}{2}\bigg{)}-\ln\Big{|}V_{0}^{[k]}+\tilde{B}_{t}^{[k]}\Big{|}+\ln\Big{|}V_{0}^{[k]}\Big{|}\Bigg{)}=0.$
(65)

  * (5)

Calculate $V_{0}^{[k+1]}$ using equation

$V_{0}^{[k+1]}=\frac{T\nu_{0}^{[k]}}{\nu_{0}^{[k]}+1}\Bigg{(}\sum_{t=1}^{T}\Big{(}V_{0}^{[k]}+\tilde{B}_{t}^{[k]}\Big{)}^{-1}\Bigg{)}^{-1}.$
(66)

  * (6)

Increase iteration count $k=k+1$ and go to step (2).

Now we consider type II conditional matrix variate $t$ distribution with
Minnesota prior. For random coefficient matrices
$A_{0,t},A_{1,t},\dots,A_{p,t}$, we assume that same prior conditions hold as
subsection 2.2, namely, for $i=1,\dots,n$ and $j=1,\dots,l$,

$\mathbb{E}\big{(}(A_{0,t})_{i,j}\big{|}\Sigma_{t},\mathcal{F}_{0}\big{)}=(1-\phi_{i})C_{i,j},$
(67)

and

$\text{Var}\big{(}(A_{0,t})_{i,j}\big{|}\Sigma_{t},\mathcal{F}_{0}\big{)}=(\sigma_{i,t}/\varepsilon_{j})^{2}$
(68)

and for $i,j=1,\dots,n$,

$\mathbb{E}\big{(}(A_{\ell,t})_{i,j}\big{|}\Sigma_{t},\mathcal{F}_{0}\big{)}=\begin{cases}\phi_{i}&\text{if}~{}~{}~{}i=j,~{}\ell=1,\\\
0&\text{if}~{}~{}~{}\text{otherwise}\end{cases},~{}~{}~{}\text{for}~{}\ell=1,\dots,p$
(69)

and

$\text{Var}\big{(}(A_{\ell,t})_{i,j}\big{|}\Sigma_{t},\mathcal{F}_{0}\big{)}=\begin{cases}\displaystyle\bigg{(}\frac{\sigma_{i,t}}{\alpha\ell^{\beta}\gamma_{i}}\bigg{)}^{2}&\text{if}~{}~{}~{}i=j,\\\
\displaystyle\bigg{(}\frac{\sigma_{i,t}}{\alpha\ell^{\beta}\gamma_{j}}\bigg{)}^{2}&\text{if}~{}~{}~{}\text{otherwise}\end{cases}~{}~{}~{}\text{for}~{}\ell=1,\dots,p,$
(70)

where $\sigma_{i,t}^{2}$ is an $(i,i)$–th element of the random covariance
matrix $\Sigma_{t}$. Then, similarly to the Algorithm 2, one obtains the
following EM algorithm, which estimate parameters of type II conditional
matrix variate $t$ distribution with Minnesota prior.

######  Algorithm 4.

(EM Estimation of Type II Conditional Matrix Variate $t$ Distribution with
Minnesota Prior).

  * (1)

Set $k=0$ and initial value of the parameter vector

$\theta^{[0]}=\Big{(}C^{[0]},\varepsilon_{1}^{[0]},\dots,\varepsilon_{l}^{[0]},\alpha^{[0]},\beta^{[0]},\gamma_{1}^{[0]},\dots,\gamma_{n}^{[0]},\nu_{0}^{[0]},V_{0}^{[0]}\Big{)}.$
(71)

  * (2)

Calculate $C^{[k+1]}$ using equation

$c^{[k+1]}:=\bigg{(}\sum_{t=1}^{T}\Big{(}I_{l}\otimes\big{(}\tilde{B}_{t}^{[k]}+V_{0}^{[k]}\big{)}^{-1}\Big{)}\bigg{)}^{-1}\sum_{t=1}^{T}\Big{(}I_{l}\otimes\big{(}\tilde{B}_{t}^{[k]}+V_{0}^{[k]}\big{)}^{-1}\Big{)}J_{\phi}\tilde{\pi}_{0|t}^{[k]}.$
(72)

  * (3)

For $i=1,\dots,n$ and $j=1,\dots,l$, calculate
$\Big{(}\varepsilon_{j}^{[k+1]}\Big{)}^{2}$,
$\Big{(}\alpha^{[k+1]}\Big{)}^{2}$, and $\Big{(}\gamma_{i}^{[k+1]}\Big{)}^{2}$
using equations

$\displaystyle\Big{(}\varepsilon_{j}^{[k+1]}\Big{)}^{2}:=\frac{nT}{\displaystyle\sum_{t=1}^{T}\text{tr}\Big{\\{}\Big{(}\tilde{\Lambda}_{0|t}^{[k]}\Delta_{\varepsilon_{j}}\otimes
I_{n}\Big{)}+\Big{(}\nu_{0}^{[k]}+1\Big{)}\tilde{\Theta}_{|t}^{[k]}\Big{(}\Delta_{\varepsilon_{j}}\otimes\big{(}V_{0}^{[k]}+\tilde{B}_{t}^{[k]}\big{)}^{-1}\Big{)}\Big{\\}}},$
(73)
$\displaystyle\Big{(}\alpha^{[k+1]}\Big{)}^{2}:=\frac{n^{2}pT}{\displaystyle\sum_{t=1}^{T}\text{tr}\Big{\\{}\Big{(}\tilde{\Lambda}_{0|t}^{[k]}\Delta_{\alpha}^{[k]}\otimes
I_{n}\Big{)}+\Big{(}\nu_{0}^{[k]}+1\Big{)}\tilde{\Theta}_{|t}^{[k]}\Big{(}\Delta_{\alpha}^{[k]}\otimes\big{(}V_{0}^{[k]}+\tilde{B}_{t}^{[k]}\big{)}^{-1}\Big{)}\Big{\\}}},$
(74)

and

$\displaystyle\Big{(}\gamma_{i}^{[k+1]}\Big{)}^{2}:=\frac{npT}{\displaystyle\sum_{t=1}^{T}\text{tr}\Big{\\{}\Big{(}\tilde{\Lambda}_{0|t}^{[k]}\Delta_{\gamma_{i}}^{[k]}\otimes
I_{n}\Big{)}+\Big{(}\nu_{0}^{[k]}+T\Big{)}\tilde{\Theta}_{|t}^{[k]}\Big{(}\Delta_{\gamma_{i}}^{[k]}\otimes\big{(}V_{0}^{[k]}+\tilde{B}_{t}^{[k]}\big{)}^{-1}\Big{)}\Big{\\}}}.$
(75)

  * (4)

For $p\geq 2$, to obtain $\beta^{[k+1]}$, solve the following nonlinear
equation for $\beta$

$\displaystyle\sum_{\ell=1}^{p}\ln(\ell)=\frac{1}{nT}\sum_{t=1}^{T}\text{tr}\Big{\\{}\Big{(}\tilde{\Lambda}_{0|t}^{[k]}\Delta_{\beta}^{[k]}\otimes
I_{n}\Big{)}+\Big{(}\nu_{0}^{[k]}+1\Big{)}\tilde{\Theta}_{|t}^{[k]}\Big{(}\Delta_{\beta}^{[k]}\otimes\big{(}V_{0}^{[k]}+\tilde{B}_{t}^{[k]}\big{)}^{-1}\Big{)}\Big{\\}}.$
(76)

  * (5)

To obtain $\nu_{0}^{[k+1]}$, solve nonlinear equation (65) for $\nu_{0}$.

  * (6)

Calculate $V_{0}^{[k+1]}$ using equation (66).

  * (7)

Increase iteration count $k=k+1$ and go to step (2).

##  4 Conclusion

Conditional matrix variate student $t$ distribution was introduced by
Battulga24g. In this paper, we provide EM algorithms, which estimate
parameters of the conditional matrix variate student $t$ distributions,
including general case and special case with Minnesota prior. Also, we
introduce a new conditional matrix variate student $t$ distribution, which is
closely related to the Battulga24g’s conditional matrix variate student $t$
distribution.

##  5 Proofs of Results

Here we provide proofs of Proposition 2 and Lemma 1.

###### Proof of Proposition 2.

The proof follows Battulga24g. For given coefficient vector $\pi_{t}$,
covariance matrix $\Sigma_{t}$, and information $\mathcal{F}_{t-1}$,
conditional density functions of the random vector of endogenous variables
$y_{t}$ is

$\displaystyle
f(y_{t}|\pi_{t},\Sigma_{t},\mathcal{F}_{t-1})=\frac{|\Sigma_{t}|^{-1/2}}{(2\pi)^{n/2}}\exp\bigg{\\{}-\frac{1}{2}\big{(}y_{t}-(\mathsf{Y}_{t}^{\prime}\otimes
I_{n})\pi_{t}\big{)}^{\prime}\Sigma_{t}^{-1}\big{(}y_{t}-(\mathsf{Y}_{t}^{\prime}\otimes
I_{n}\big{)}\pi_{t})\bigg{\\}}.$ (77)

Since conditional on $\mathcal{F}_{0}$, $(\pi_{t},\Sigma_{t})$ is independent
of $\bar{y}_{t-1}$, conditional density functions of the random coefficient
vector $\pi_{t}$ and covariance matrix $\Sigma_{t}$ are given by

$\displaystyle
f(\pi_{t}|\Sigma_{t},\mathcal{F}_{t-1})=\frac{|\Sigma_{t}|^{-d/2}}{(2\pi)^{nd/2}|\Lambda_{0}|^{n/2}}\exp\bigg{\\{}-\frac{1}{2}\big{(}\pi_{t}-\pi_{0}\big{)}^{\prime}\big{(}\Lambda_{0}^{-1}\otimes\Sigma_{t}^{-1}\big{)}\big{(}\pi_{t}-\pi_{0}\big{)}\bigg{\\}}$
(78)

and

$f(\Sigma_{t}|\mathcal{F}_{t-1})=\frac{|V_{0}|^{\nu_{0}/2}}{\Gamma_{n}(\nu_{0}/2)2^{n\nu_{0}/2}}|\Sigma_{t}|^{-(\nu_{0}+n+1)/2}\exp\bigg{\\{}-\frac{1}{2}\text{tr}\Big{(}V_{0}\Sigma_{t}^{-1}\Big{)}\bigg{\\}}.$
(79)

By the completing square method, a joint conditional density function of the
random vectors $y_{t}$ and $\pi_{t}$ is

$\displaystyle f(y_{t},\pi_{t}|\Sigma_{t},\mathcal{F}_{t-1})$ $\displaystyle=$
$\displaystyle
c_{1}|\Sigma_{t}|^{-(1+d)/2}\exp\bigg{\\{}-\frac{1}{2}\Big{(}\pi_{t}-\tilde{\pi}_{0|t}\Big{)}^{\prime}\big{(}\tilde{\Lambda}_{0|t}^{-1}\otimes\Sigma_{t}^{-1}\big{)}\Big{(}\pi_{t}-\tilde{\pi}_{0|t}\Big{)}\bigg{\\}}$
$\displaystyle\times$
$\displaystyle\exp\bigg{\\{}-\frac{1}{2}\Big{(}y_{t}^{\prime}\Sigma_{t}^{-1}y_{t}+\pi_{0}^{\prime}\big{(}\Lambda_{0}^{-1}\otimes\Sigma_{t}^{-1}\big{)}\pi_{0}-\tilde{\pi}_{0|t}^{\prime}\big{(}\tilde{\Lambda}_{0|t}^{-1}\otimes\Sigma_{t}^{-1}\big{)}\tilde{\pi}_{0|t}\Big{)}\bigg{\\}},$

where
$\tilde{\Lambda}_{0|t}^{-1}:=\mathsf{Y}_{t}\mathsf{Y}_{t}^{\prime}+\Lambda_{0}^{-1}$
is a $(d\times d)$ matrix,
$\tilde{\pi}_{0|t}:=\big{(}(\tilde{\Lambda}_{0|t}\mathsf{Y}_{t})\otimes\Sigma_{t}^{-1}\big{)}y_{t}+\big{(}(\tilde{\Lambda}_{0|t}\Lambda_{0}^{-1})\otimes\Sigma_{t}^{-1}\big{)}\pi_{0}$
is an $([nd]\times 1)$ vector, and normalizing coefficient equals

$c_{1}:=\frac{1}{(2\pi)^{n(1+d)/2}|\tilde{\Lambda}_{0|t}|^{n/2}}.$ (81)

If we integrate the above joint density function with respect to the vector
$\pi_{t}$, then an integral, corresponding to the first exponential is
proportional to
$|\tilde{\Lambda}_{0|t}\otimes\Sigma_{t}|^{1/2}=|\mathsf{Y}_{t}\mathsf{Y}_{t}^{\prime}+\Lambda_{0|t}^{-1})|^{n/2}|\Sigma_{t}|^{d/2}$.
Therefore, we have that

$\displaystyle
f(y_{t}|\Sigma_{t},\mathcal{F}_{t-1})=c_{2}|\Sigma_{t}|^{-1/2}\exp\bigg{\\{}-\frac{1}{2}\Big{(}y_{t}^{\prime}\Sigma_{t}^{-1}y_{t}+\pi_{0}^{\prime}\big{(}\Lambda_{0}^{-1}\otimes\Sigma_{t}^{-1}\big{)}\pi_{0}-\tilde{\pi}_{0|t}^{\prime}\big{(}\tilde{\Lambda}_{0|t}^{-1}\otimes\Sigma_{t}^{-1}\big{)}\tilde{\pi}_{0|t}\Big{)}\bigg{\\}},$
(82)

where the normalizing coefficient equals

$c_{2}:=\frac{1}{(2\pi)^{n/2}|\Lambda_{0}|^{n/2}|\tilde{\Lambda}_{0|t}^{-1}|^{n/2}}.$
(83)

Hence, according to the well–known formula that for suitable matrices
$A,B,C,D$,

$\text{vec}(A)^{\prime}(B\otimes
C)\text{vec}(D)=\text{tr}(DB^{\prime}A^{\prime}C),$ (84)

we find that

$f(y_{t}|\Sigma_{t},\mathcal{F}_{t-1})=c_{2}|\Sigma_{t}|^{-1/2}\exp\bigg{\\{}-\frac{1}{2}\text{tr}\big{(}\tilde{B}_{t}\Sigma_{t}^{-1}\big{)}\bigg{\\}}.$
(85)

Thus, it follows from equations (77) and (79) that a joint conditional density
of the random vector $y_{t}$ and random matrix $\Sigma_{t}$ is

$f(y_{t},\Sigma_{t}|\mathcal{F}_{t-1})=c_{3}|\Sigma_{t}|^{-(\nu_{0}+n+2)/2}\exp\bigg{\\{}-\frac{1}{2}\text{tr}\Big{(}\big{(}\tilde{B}_{t}+V_{0}\big{)}\Sigma_{t}^{-1}\Big{)}\bigg{\\}},$
(86)

where the normalizing coefficient equals

$c_{3}:=\frac{1}{(2\pi)^{n/2}}\frac{|V_{0}|^{\nu_{0}/2}}{|\Lambda_{0}|^{n/2}|\tilde{\Lambda}_{0|t}^{-1}|^{n/2}\Gamma_{n}(\nu_{0}/2)2^{n\nu_{0}/2}}.$
(87)

Consequently, a density function of the random vector $y_{t}$ is given by

$\displaystyle
f(\bar{y}_{t}|\bar{s}_{t},\mathcal{F}_{0})=\int_{\Sigma_{t}>0}f(y_{t},\Sigma_{t}|\mathcal{F}_{t-1})d\Sigma_{t}=c_{3}\prod_{k=1}^{r_{t}}\frac{\Gamma_{n}\big{(}\nu_{0t|t}/2\big{)}2^{n(\nu_{0}+1)/2}}{\big{|}\tilde{B}_{t}+V_{0}\big{|}^{(\nu_{0}+1)/2}}.$
(88)

By the completing square method, the matrix $\tilde{B}_{t}$ can be written by

$\displaystyle\tilde{B}_{t}$ $\displaystyle=$
$\displaystyle\big{(}y_{t}-\pi_{0}^{\circ}\Lambda_{0}^{-1}\tilde{\Lambda}_{0|t}\mathsf{Y}_{t}\phi_{t}^{-1}\big{)}\phi_{t}^{-1}\big{(}y_{t}-\pi_{0}^{\circ}\Lambda_{0}^{-1}\tilde{\Lambda}_{0|t}\mathsf{Y}_{t}\phi_{t}^{-1}\big{)}^{\prime}$
(89) $\displaystyle-$
$\displaystyle\pi_{0}^{\circ}\Lambda_{0}^{-1}\tilde{\Lambda}_{0|t}\mathsf{Y}_{t}\phi_{t}^{-1}(\mathsf{Y}_{t})^{\prime}\tilde{\Lambda}_{0|t}\Lambda_{0}^{-1}(\pi_{0}^{\circ})^{\prime}+\pi_{0}^{\circ}\Lambda_{0}^{-1}(\pi_{0}^{\circ})^{\prime}-\pi_{0}^{\circ}\Lambda_{0}^{-1}\tilde{\Lambda}_{0|t}\Lambda_{0}^{-1}(\pi_{0}^{\circ})^{\prime},$

where
$\phi_{t}:=1-\mathsf{Y}_{t}^{\prime}\tilde{\Lambda}_{0|t}\mathsf{Y}_{t}$. We
consider the following product

$1_{t}:=\big{(}1+\mathsf{Y}_{t}^{\prime}\Lambda_{0}\mathsf{Y}_{t}\big{)}\big{(}1-\mathsf{Y}_{t}^{\prime}\tilde{\Lambda}_{0|t}\mathsf{Y}_{t}\big{)}.$
(90)

It equals

$\displaystyle 1_{t}$ $\displaystyle=$ $\displaystyle
1+\mathsf{Y}_{t}^{\prime}\Lambda_{0}\mathsf{Y}_{t}-\mathsf{Y}_{t}^{\prime}\tilde{\Lambda}_{0|t}\mathsf{Y}_{t}-\mathsf{Y}_{t}^{\prime}\Lambda_{0}\mathsf{Y}_{t}\mathsf{Y}_{t}^{\prime}\tilde{\Lambda}_{0|t}\mathsf{Y}_{t}.$
(91)

If we add and subtract the matrix $\Lambda_{0}^{-1}$ into the term
$\mathsf{Y}_{t}\mathsf{Y}_{t}^{\prime}$ in the last line of the above
equation, then $1_{t}$ equals 1. Consequently,
$1+\mathsf{Y}_{t}^{\prime}\Lambda_{0}\mathsf{Y}_{t}$ is a reciprocal of
$\phi_{t}$, that is,

$\phi_{t}^{-1}=1+\mathsf{Y}_{t}^{\prime}\Lambda_{0}\mathsf{Y}_{t}.$ (92)

Since it takes a positive value, the matrix $B_{t}$ is a positive
semi–definite matrix. Now, we consider the term
$\Lambda_{0}^{-1}\tilde{\Lambda}_{0|t}\mathsf{Y}_{t}\phi_{t}^{-1}$ in the
first line in equation (89). Similarly as before, by adding and subtracting
$\Lambda_{0}^{-1}$ into the term $\mathsf{Y}_{t}\mathsf{Y}_{t}^{\prime}$, one
obtains that

$\Lambda_{0}^{-1}\tilde{\Lambda}_{0|t}\mathsf{Y}_{t}\phi_{t}^{-1}=\mathsf{Y}_{t}.$
(93)

Consequently, the second line of equation (89) equals zero. Let
$\Lambda_{0t}^{1/2}$ be the Cholesky factor of the matrix $\Lambda_{0}$, i.e.,
$\Lambda_{0}=\big{(}\Lambda_{0}^{1/2}\big{)}^{\prime}\Lambda_{0}^{1/2}$. Then,
according to the Sylvester’s determinant theorem, see Lutkepohl05,
$\phi_{t}^{-1}$ equals

$|\phi_{t}^{-1}|=\big{|}1+\big{(}\Lambda_{0}^{1/2}\mathsf{Y}_{t}\big{)}^{\prime}\Lambda_{0}^{1/2}\mathsf{Y}_{t}\big{|}=\big{|}I_{d}+\Lambda_{0}^{1/2}\mathsf{Y}_{t}\mathsf{Y}_{t}^{\prime}\big{(}\Lambda_{0}^{1/2}\big{)}^{\prime}\big{|}.$
(94)

That completes the proof of the Proposition. ∎

###### Proof of Lemma 1.

By the conditional probability formula and law of total probability, the
density function $f(\pi_{t},\Sigma_{t}|\mathcal{F}_{T})$ is represented by

$\displaystyle
f(\pi_{t},\Sigma_{t}|\mathcal{F}_{T})=\frac{\int_{\pi_{-t},\Sigma_{-t}}f(y,\pi,\Sigma|\mathcal{F}_{0})d\pi_{-t}d\Sigma_{-t}}{f(y|\mathcal{F}_{0})}$
(95)

for $t=1,\dots,T$, where $\pi_{-t}$ is a $(d\times[T-1])$ matrix, which
excludes the vector $\pi_{t}$ from a matrix $[\pi_{1}:\dots:\pi_{T}]$ and
$\Sigma_{-t}$ is an $(n\times[(T-1)n])$ matrix, which excludes the matrix
$\Sigma_{t}$ from a matrix $[\Sigma_{1}:\dots:\Sigma_{T}$. Due to the
conditional probability formula and the assumption that for given initial
information $\mathcal{F}_{0}$,
$(\pi_{1},\Sigma_{1}),\dots,(\pi_{T},\Sigma_{T})$ are independent, the
numerator of the above equation equals

$\displaystyle\int_{\pi_{-t},\Sigma_{-t}}\prod_{i=1}^{T}f(y_{i}|\pi_{i},\Sigma_{i},\mathcal{F}_{i-1})f(\pi_{i},\Sigma_{i}|\mathcal{F}_{0})d\pi_{-t}d\Sigma_{-t}$
$\displaystyle=\prod_{i=1,i\neq
t}^{T}f(y_{i}|\mathcal{F}_{i-1})f(y_{t}|\pi_{t},\Sigma_{t},\mathcal{F}_{t-1})f(\pi_{t},\Sigma_{t}|\mathcal{F}_{0}).$
(96)

On the other hand, by the conditional probability formula, the denominator of
equation (95) equals $\prod_{i=1}^{T}f(y_{i}|\mathcal{F}_{i-1})$.
Consequently, since conditional on the initial information $\mathcal{F}_{0}$,
$(\pi_{t},\Sigma_{t})$ is independent of the random vector $\bar{y}_{t-1}$,
one obtains equation (59). ∎

## References

  * Bańbura . (2010) Banbura10Bańbura, M., Giannone, D. Reichlin, L.  2010\.  Large Bayesian Vector Autoregressions Large Bayesian Vector Autoregressions.  Journal of Applied Econometrics25171–92.
  * Battulga (20241) Battulga24gBattulga, G.  20241\.  Bayesian Markov-Switching Vector Autoregressive Process Bayesian Markov-Switching Vector Autoregressive Process.  Available at: https://arxiv.org/abs/2404.11235.
  * Battulga (20242) Battulga24cBattulga, G.  20242\.  The Log Private Company Valuation Model The Log Private Company Valuation Model.  to appear in Numerical Algebra, Control & Optimization.  Available at: https://arxiv.org/abs/2206.09666.
  * Dempster . (1977) Dempster77Dempster, AP., Laird, NM. Rubin, DB.  1977\.  Maximum likelihood from incomplete data via the EM algorithm Maximum likelihood from incomplete data via the EM algorithm.  Journal of the royal statistical society: series B (methodological)3911–22.
  * Hamilton (1990) Hamilton90Hamilton, JD.  1990\.  Analysis of Time Series Subject to Changes in Regime Analysis of Time Series Subject to Changes in Regime.  Journal of Econometrics451-239–70.
  * Hamilton (1994) Hamilton94Hamilton, JD.  1994\.  Time Series Econometrics Time Series Econometrics.  Princeton University Press, Princeton.
  * Kalman (1960) Kalman60Kalman, RE.  1960\.  A New Approach to Linear Filtering and Prediction Problems A New Approach to Linear Filtering and Prediction Problems.  Journal of Basic Engineering82135–44.
  * Litterman (1979) Litterman79Litterman, R.  1979\.  Techniques of forecasting using vector autoregressions Techniques of forecasting using vector autoregressions.  Federal Reserve of Minneapolis Working Paper 115.
  * Lütkepohl (2005) Lutkepohl05Lütkepohl, H.  2005\.  New Introduction to Multiple Time Series Analysis New Introduction to Multiple Time Series Analysis (2 ).  Springer Berlin Heidelberg.
  * McNeil . (2005) McNeil05McNeil, AJ., Frey, R. Embrechts, P.  2005\.  Quantitative Risk Management: Concepts, Techniques and Tools Quantitative Risk Management: Concepts, Techniques and Tools.  Princeton University Press.
  * Miranda-Agrippino Ricco (2018) Miranda18Miranda-Agrippino, S. Ricco, G.  2018\.  Bayesian vector autoregressions Bayesian vector autoregressions.
  * Nguyen (2023) Nguyen23Nguyen, D.  2023\.  An in depth introduction to variational bayes note An in depth introduction to variational bayes note.  Available at SSRN.
  * Sims (1980) Sims80Sims, CA.  1980\.  Macroeconomics and Reality Macroeconomics and reality.  Econometrica4811–48.
  * Sun . (2010) Sun10Sun, J., Kabán, A. Garibaldi, JM.  2010\.  Robust mixture clustering using Pearson type VII distribution Robust mixture clustering using Pearson type VII distribution.  Pattern Recognition Letters31162447–2454.
  * Thompson . (2020) Thompson20Thompson, GZ., Maitra, R., Meeker, WQ. Bastawros, AF.  2020\.  Classification With the Matrix-Variate-t Distribution Classification With the Matrix-Variate-t Distribution.  Journal of Computational and Graphical Statistics293668–674.
  * Tiao Box (1981) Tiao81Tiao, GC. Box, GE.  1981\.  Modeling multiple time series with applications Modeling multiple time series with applications.  Journal of the American Statistical Association76376802–816.


MIDLMedical Imaging with Deep Learning 2019 MIDL 2019 – Extended Abstract
Track Christian Payer1,2<EMAIL_ADDRESS>
1 Institute of Computer Graphics and Vision, Graz University of Technology,
Graz, Austria
2 Ludwig Boltzmann Institute for Clinical Forensic Imaging, Graz, Austria and
Darko Štern2<EMAIL_ADDRESS>
Horst Bischof1<EMAIL_ADDRESS>
Martin Urschler2,1<EMAIL_ADDRESS>
# Integrating Spatial Configuration into Heatmap Regression Based CNNs for
Landmark Localization

###### Abstract

In many medical image analysis applications, often only a limited amount of
training data is available, which makes training of convolutional neural
networks (CNNs) challenging. In this work on anatomical landmark localization,
we propose a CNN architecture that learns to split the localization task into
two simpler sub-problems, reducing the need for large training datasets. Our
fully convolutional SpatialConfiguration-Net (SCN) dedicates one component to
locally accurate but ambiguous candidate predictions, while the other
component improves robustness to ambiguities by incorporating the spatial
configuration of landmarks. In our experimental evaluation, we show that the
proposed SCN outperforms related methods in terms of landmark localization
error on size-limited datasets.

###### keywords:

anatomical landmarks, localization, heatmap regression, spatial configuration

##  1 Introduction

Localization of anatomical landmarks is an important step in medical image
analysis, e.g., in segmentation [Beichel et al. (2005)], or registration
[Johnson and Christensen (2002)]. Unfortunately, locally similar structures
often introduce difficulties due to ambiguity into landmark localization. To
deal with these difficulties, machine learning based approaches often combine
local landmark predictions with explicit handcrafted graphical models, aiming
to restrict predictions to feasible spatial configurations. Thus, the landmark
localization problem is simplified by separating the task into two successive
steps. The first step is dedicated to locally accurate but potentially
ambiguous predictions, while in the second step graphical models [Cootes et
al. (1995), Felzenszwalb and Huttenlocher (2005)] eliminate ambiguities.

Figure 1: Landmark localization by regressing heatmaps for each landmark in
our end-to-end trained fully convolutional SpatialConfiguration-Net (SCN).

Recent advances in computer vision and medical imaging have mainly been driven
by convolutional neural networks (CNNs) due to their superior capabilities to
automatically learn important image features [LeCun et al. (2015)].
Unfortunately, CNNs typically need large amounts of training data. Especially
in medical imaging, this requirement is hard to fulfill, due to ethical and
financial concerns as well as time consuming expert annotations.

In this work, we show that the amount of required training data can be reduced
with our proposed two-component SpatialConfiguration-Net (SCN), which follows
the idea of handcrafted graphical models to split landmark localization into
two successive steps. This extended abstract gives a short overview of the key
concepts of our journal paper published in [Payer et al. (2019)], while we
refer the reader to the full paper for more detailed descriptions and more
extensive evaluations on a variety of datasets.

##  2 Method

Our method for landmark localization is based on regressing heatmap images
(Tompson et al., 2014), which encode the pseudo-probability of a landmark
being located at a certain pixel position. With $N$ being the total number of
landmarks, we define the target heatmap image of a landmark $L_{i}$,
$i=\\{1,...,N\\}$ as the $d$-dimensional Gaussian function
${g_{i}(\vec{x}):\mathbb{R}^{d}\rightarrow\mathbb{R}}$ centered at the target
landmark’s groundtruth coordinate
${\vec{\overset{\ast}{x}}_{i}\in\mathbb{R}^{d}}$.

The network is set up to regress $N$ heatmaps simultaneously by minimizing the
differences between predicted heatmaps $h_{i}(\vec{x})$ and the corresponding
target heatmaps $g_{i}(\vec{x})$ in an end-to-end manner Ronneberger et al.
(2015); Shelhamer et al. (2017). In network inference, we obtain the predicted
coordinate $\vec{\hat{x}}_{i}\in\mathbb{R}^{d}$ of each landmark $L_{i}$ by
taking the coordinate, where the heatmap has its highest value.

###  2.1 SpatialConfiguration-Net

The fundamental concept of the SpatialConfiguration-Net (SCN) is the
interaction between its two components (see Fig. 1). The first component takes
the image as input to generate locally accurate but potentially ambiguous
local appearance heatmaps $h^{\text}{LA}_{i}(\vec{x})$. Motivated by
handcrafted graphical models for eliminating these potential ambiguities, the
second component takes the predicted candidate heatmaps
$h^{\text}{LA}_{i}(\vec{x})$ as input to generate inaccurate but unambiguous
spatial configuration heatmaps $h^{\text}{SC}_{i}(\vec{x})$.

For $N$ landmarks, the set of predicted heatmaps
$\mathbb{H}=\\{h_{i}(\vec{x})\;|\;i=1\dots N\\}$ is obtained by element-wise
multiplication $\odot$ of the corresponding heatmap outputs
$h^{\text}{LA}_{i}(\vec{x})$ and $h^{\text}{SC}_{i}(\vec{x})$ of the two
components:

$h_{i}(\vec{x})=h^{\text}{LA}_{i}(\vec{x})\odot h^{\text}{SC}_{i}(\vec{x}).$
(1)

This multiplication is crucial for the SCN, as it forces both of its
components to generate a response on the location of the target landmark
$\vec{\overset{\ast}{x}}_{i}$, i.e., both $h^{\text}{LA}_{i}(\vec{x})$ and
$h^{\text}{SC}_{i}(\vec{x})$ deliver responses for $\vec{x}$ close to
$\vec{\overset{\ast}{x}}_{i}$, while on all other locations one component may
have a response as long as the other one does not have one.

##  3 Experiments and Results

We evaluate our proposed SCN on a dataset of 895 radiographs of left hands
with 37 annotated characteristic landmarks on finger tips and bone joints. We
compare our SCN to state-of-the-art random regression forests Ebner et al.
(2014); Lindner et al. (2015); Štern et al. (2016); Urschler et al. (2018),
our previous CNN-based method of Payer et al. (2016), and our implementation
of a localization U-Net for heatmap regression. Results of the image-specific
point-to-point errors for three-fold cross validation of the 895 radiographs
are shown in Fig. 2. When using all training images, our SCN outperforms all
other compared methods. Additionally, when drastically reducing the number of
training images to 100, 50, and 10, respectively, our SCN greatly outperforms
the localization U-Net. This confirms that splitting the localization task
into predicting accurate but potentially ambiguous local appearance heatmaps
and inaccurate but unambiguous spatial configuration heatmaps is especially
useful when dealing with only limited amounts of training data.

[all training images] [100/50/10 training images]

Figure 2: Cumulative distributions of the point-to-point error for 895
radiographs. (a) shows results compared with other state-of-the-art methods.
(b) shows results of SCN and localization U-Net for reduced numbers of
training images.

##  4 Conclusion

In conclusion, we have shown how to combine information of local appearance
and spatial configuration into a single end-to-end trained network for
landmark localization. Our generic architecture achieves state-of-the-art
results in terms of localization error, even when only limited amounts of
training images are available. We are currently looking into extending our SCN
regarding occluded structures and multi-object localization, and into adapting
our SCN for semantic segmentation problems (see Payer et al. (2018) for
preliminary results), where structural constraints may be used in a similar
manner.

## References

  * Beichel et al. (2005) Reinhard Beichel, Horst Bischof, Franz Leberl, and Milan Sonka.  Robust Active Appearance Models and Their Application to Medical Image Analysis.  _IEEE Trans. Med. Imaging_ , 24(9):1151–1169, sep 2005.  10.1109/TMI.2005.853237.
  * Cootes et al. (1995) Tim F. Cootes, Christopher J. Taylor, David H. Cooper, and Jim Graham.  Active Shape Models-Their Training and Application.  _Comput. Vis. Image Underst._ , 61(1):38–59, jan 1995.  10.1006/cviu.1995.1004.
  * Ebner et al. (2014) Thomas Ebner, Darko Štern, René Donner, Horst Bischof, and Martin Urschler.  Towards Automatic Bone Age Estimation from MRI: Localization of 3D Anatomical Landmarks.  In _Proc. Med. Image Comput. Comput. Interv._ , pages 421–428. Springer, 2014.  10.1007/978-3-319-10470-6_53.
  * Felzenszwalb and Huttenlocher (2005) Pedro F. Felzenszwalb and Daniel P. Huttenlocher.  Pictorial Structures for Object Recognition.  _Int. J. Comput. Vis._ , 61(1):55–79, 2005.  10.1023/B:VISI.0000042934.15159.49.
  * Johnson and Christensen (2002) Hans J. Johnson and Gary E. Christensen.  Consistent Landmark and Intensity-Based Image Registration.  _IEEE Trans. Med. Imaging_ , 21(5):450–461, 2002\.  10.1109/TMI.2002.1009381.
  * LeCun et al. (2015) Yann LeCun, Yoshua Bengio, and Geoffrey Hinton.  Deep Learning.  _Nature_ , 521(7553):436–444, 2015.  10.1038/nature14539.
  * Lindner et al. (2015) Claudia Lindner, Paul A. Bromiley, Mircea C. Ionita, and Tim F. Cootes.  Robust and Accurate Shape Model Matching Using Random Forest Regression-Voting.  _IEEE Trans. Pattern Anal. Mach. Intell._ , 37(9):1862–1874, sep 2015.  10.1109/TPAMI.2014.2382106.
  * Payer et al. (2016) Christian Payer, Darko Štern, Horst Bischof, and Martin Urschler.  Regressing Heatmaps for Multiple Landmark Localization Using CNNs.  In _Proc. Med. Image Comput. Comput. Interv._ , pages 230–238. Springer, 2016.  10.1007/978-3-319-46723-8_27.
  * Payer et al. (2018) Christian Payer, Darko Štern, Horst Bischof, and Martin Urschler.  Multi-label Whole Heart Segmentation Using CNNs and Anatomical Label Configurations.  In _Stat. Atlases Comput. Model. Hear. ACDC MMWHS Challenges. STACOM 2017._ , pages 190–198. Springer, 2018.  10.1007/978-3-319-75541-0_20.
  * Payer et al. (2019) Christian Payer, Darko Štern, Horst Bischof, and Martin Urschler.  Integrating Spatial Configuration into Heatmap Regression Based CNNs for Landmark Localization.  _Med. Image Anal._ , 54:207–219, may 2019.  10.1016/j.media.2019.03.007.
  * Ronneberger et al. (2015) Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-Net: Convolutional Networks for Biomedical Image Segmentation.  In _Proc. Med. Image Comput. Comput. Interv._ , pages 234–241. Springer, 2015.  10.1007/978-3-319-24574-4_28.
  * Shelhamer et al. (2017) Evan Shelhamer, Jonathan Long, and Trevor Darrell.  Fully Convolutional Networks for Semantic Segmentation.  _IEEE Trans. Pattern Anal. Mach. Intell._ , 39(4):640–651, apr 2017.  10.1109/TPAMI.2016.2572683.
  * Štern et al. (2016) Darko Štern, Thomas Ebner, and Martin Urschler.  From Local to Global Random Regression Forests: Exploring Anatomical Landmark Localization.  In _Proc. Med. Image Comput. Comput. Interv._ , pages 221–229. Springer, 2016.  10.1007/978-3-319-46723-8_26.
  * Tompson et al. (2014) Jonathan Tompson, Arjun Jain, Yann LeCun, and Christoph Bregler.  Joint Training of a Convolutional Network and a Graphical Model for Human Pose Estimation.  In _Adv. Neural Inf. Process. Syst._ , pages 1799–1807, 2014.
  * Urschler et al. (2018) Martin Urschler, Thomas Ebner, and Darko Štern.  Integrating Geometric Configuration and Appearance Information into a Unified Framework for Anatomical Landmark Localization.  _Med. Image Anal._ , 43:23–36, jan 2018.  10.1016/j.media.2017.09.003.