# Bayesian Markov–Switching Vector Autoregressive Process

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

###### Abstract

This study introduces marginal density functions of the general Bayesian
Markov–Switching Vector Autoregressive (MS–VAR) process. In the case of the
Bayesian MS–VAR process, we provide closed–form density functions and
Monte–Carlo simulation algorithms, including the importance sampling method.
The Monte–Carlo simulation method departs from the previous simulation methods
because it removes the duplication in a regime vector.

Keywords: Bayesian MS–VAR process, Monte–Carlo simulation method, Stochastic
DDM.

##  1 Introduction

Classic Vector Autoregressive (VAR) process was proposed by Sims80 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. Research works have shown that BVAR is an appropriate tool for
modeling large data sets; for example, see Banbura10.

Sudden and dramatic changes in the financial market and economy are caused by
events such as wars, market panics, or significant changes in government
policies. To model those events, some authors used regime–switching models.
The regime–switching model was introduced by seminal works of
Hamilton89,Hamilton90 (see also books of Hamilton94 and Krolzig97), and the
model is hidden Markov model with dependencies; see Zucchini16. However,
Markov regime–switching models have been introduced before Hamilton (1989),
see, Goldfeld73, Quandt58, and Tong83. 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. The model
fits some financial data well and has become popular in financial modeling,
including equity options, bond prices, and others.

A model that considers all of the above is the Bayesian Markov–Switching VAR
(MS–VAR) process. Its applications in finance can be found in Battulga23b,
Battulga24f, and Battulga24a. In some existing option pricing models, the
underlying asset price is governed by some stochastic process, and economic
variables such as GDP, inflation, unemployment rate, and so on are not taken
into account. For this reason, the author has developed option pricing models,
depending on economic variables. Applying the Bayesian MS–VAR process, with
direct calculation and change of probability measure for some frequently used
options, Battulga24a derived pricing formulas. Also, the author used the
Bayesian MS–VAR process to price equity–linked life insurance products and
rainbow options, see Battulga24f and Battulga23b.

Monte–Carlo simulation methods using the Gibbs sampling algorithm for Bayesian
MS–VAR process are proposed by some authors. In particular, the Monte–Carlo
simulation method of the Bayesian MS–AR(p) process is provided by Albert93,
and its multidimensional extension is given by Krolzig97. In this paper, we
introduce a new Monte–Carlo simulation method that removes duplication in a
regime vector. We also introduce importance sampling method to estimate
probability of rare event, which corresponds to endogenous variables.
Importance sampling is an effective variance reduction technique for studying
the rare events. Glasserman00 used the importance sampling method to model
portfolio loss random variable by using approximation. Also, Glasserman05b
study a loss random variable of credit portfolio applying the method, see also
McNeil05.

Dividend discount models (DDMs), first introduced by Williams38, are common
methods for stock valuation. The basic idea is that a stock price of a firm is
equal to a sum of dividend paid by the firm and the stock price of the firm,
which correspond to the next period and which are discounted at required rate
of return on stock. As the outcome of DDMs depends crucially on dividend
payment forecasts, most research in the last few decades has been around the
proper estimations of dividend development. To model the dividends of a firm,
Battulga22a used the compound Poisson process. Also, parameter estimation of
DDMs is a challenging task. Battulga22a introduced parameter estimation
methods for practically popular DDMs. Battulga23a provided parameter
estimation methods of the required rate of returns for public and private
companies. Under the normal framework, Battulga22b obtained pricing and
hedging formulas for the European options and equity–linked life insurance
products by introducing a DDM with regime–switching process. A review of some
existing DDMs, including deterministic and stochastic models can be found in
dAmico20a.

The rest of the paper is organized as follows: In Section 2, for the general
Bayesian MS–VAR process, we obtain some conditional density functions, which
are helpful for general Monte–Carlo simulation. Section 3 is dedicated to
studying a special case of the process, where we obtain closed-form
conditional density functions of our model’s random components. Some of the
conditional density functions have not been explored before. In Section 3, we
provide Monte–Carlo simulation methods, including the importance sampling
method. Section 4 gives numerical results on three companies, listed in the
S&P 500 index. Finally, Section 5 concludes the study.

##  2 Bayesian MS–VAR($p$) process

Let $(\Omega,\mathcal{H}_{T},\mathbb{P})$ be a complete probability space,
where $\mathbb{P}$ is a given physical or real–world probability measure.
Other elements of the probability space will be defined below. To introduce a
regime–switching, we assume that $\\{s_{t}\\}_{t=1}^{T}$ is a homogeneous
Markov chain with $N$ state and $\mathsf{P}:=\\{p_{ij}\\}_{i=0,j=1}^{N}$ is a
random transition probability matrix, including an initial probability vector,
where $\\{p_{0j}\\}_{j=1}^{N}$ is the initial probability vector. We consider
a Bayesian Markov–Switching Vector Autoregressive process of $p$ order
(MS–VAR($p$)), which is given by the following equation

$y_{t}=A_{0,s_{t}}\psi_{t}+A_{1,s_{t}}y_{t-1}+\dots+A_{p,s_{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,s_{t}}$ is an $(n\times l)$ random coefficient matrix at regime
$s_{t}$, corresponding to the vector of exogenous variables, for
$i=1,\dots,p$, $A_{i,s_{t}}$ are $(n\times n)$ random coefficient matrices at
regime $s_{t}$, corresponding to $y_{t-1},\dots,y_{t-p}$. Equation (1) can be
written by

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

where $\Pi_{s_{t}}:=[A_{0,s_{t}}:A_{1,s_{t}}:\dots:A_{p,s_{t}}]$ is an
$(n\times d)$ random coefficient matrix with $d:=l+np$ at regime $s_{t}$,
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.

For the residual process $\xi_{t}$, we assume that it has
$\xi_{t}:=\Sigma_{s_{t}}^{1/2}\varepsilon_{t}$, $t=1,\dots,T$ representation,
see Lutkepohl05 and McNeil05, where $\Sigma_{s_{t}}^{1/2}$ is a Cholesky
factor of a positive definite $(n\times n)$ random matrix $\Sigma_{s_{t}}$,
which is measurable with respect to $\sigma$–field $\mathcal{H}_{t-1}$,
defined below and depends on $(n_{*}\times d_{*})$ random coefficient matrix
$\Gamma_{s_{t}}:=[B_{0,s_{t}}:B_{1,s_{t}}:\dots:B_{p_{*}+q_{*},s_{t}}]$ with
$d_{*}:=l_{*}+n_{*}(p_{*}+q_{*})$. Here $B_{0,s_{t}}$ is an $(n_{*}\times
l_{*})$ random matrix, for $i=1,\dots,p_{*}+q_{*}$, $B_{i,s_{t}}$ are
$(n_{*}\times n_{*})$ random matrices, 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}$. Then, in
particular, for multivariate GARCH process of $(p_{*},q_{*})$ order,
dependence of $\Sigma_{s_{t}}^{1/2}$ on $\Gamma_{s_{t}}$ is given by

$\text{vech}\big{(}\Sigma_{s_{t}}\big{)}=B_{0,s_{t}}+\sum_{i=1}^{p_{*}}B_{i,s_{t}}\text{vech}\big{(}\xi_{t-i}\xi_{t-i}^{\prime}\big{)}+\sum_{j=1}^{q_{*}}B_{p_{*}+j,s_{t}}\text{vech}(\Sigma_{s_{t}-j}),$
(3)

where $B_{0,s_{t}}$ and $B_{i,s_{t}}$ for $i=1,\dots,p_{*}+q_{*}$ are suitable
$([n(n+1)/2]\times 1)$ random vector and suitable
$([n(n+1)/2]\times[n(n+1)/2])$ matrices, respectively, and the vech is an
operator that stacks elements on and below a main diagonal of a square matrix.

Let us introduce stacked vectors and matrices:
$y:=(y_{1}^{\prime},\dots,y_{T}^{\prime})^{\prime}$,
$s:=(s_{1},\dots,s_{T})^{\prime}$, $\Pi_{s}:=[\Pi_{s_{1}}:\dots:\Pi_{s_{T}}]$,
and $\Gamma_{s}:=[\Gamma_{s_{1}}:\dots:\Gamma_{s_{T}}]$. We also assume that
the strong white noise process $\\{\varepsilon_{t}\\}_{t=1}^{T}$ is
independent of the random coefficient matrices $\Pi_{s}$ and $\Gamma_{s}$,
random transition matrix $\mathsf{P}$, and regime vector $s$ conditional on
initial information
$\mathcal{F}_{0}:=\sigma(y_{1-p}^{\prime},\dots,y_{0}^{\prime},\psi_{1},\dots,\psi_{T},\Sigma_{1-q_{*}},\dots,\Sigma_{0})$.
Here for a generic random vector $X$, $\sigma(X)$ denotes a $\sigma$–field
generated by the random vector $X$, $\Sigma_{1-q_{*}},\dots,\Sigma_{0}$ is
initial values of the random matrix process $\Sigma_{s_{t}}$,
$\psi_{1},\dots,\psi_{T}$ are values of exogenous variables and they are known
at time zero. We further suppose that the transition probability matrix
$\mathsf{P}$ is independent of the random coefficient matrices $\Pi_{s}$ and
$\Gamma_{s}$ given initial information $\mathcal{F}_{0}$ and regime vector
$s$.

To ease of notations, for a generic vector
$o=(o_{1}^{\prime},\dots,o_{T}^{\prime})^{\prime}$, we denote its first $t$
and last $T-t$ sub vectors by $\bar{o}_{t}$ and $\bar{o}_{t}^{c}$,
respectively, that is,
$\bar{o}_{t}:=(o_{1}^{\prime},\dots,o_{t}^{\prime})^{\prime}$ and
$\bar{o}_{t}^{c}:=(o_{t+1}^{\prime},\dots,o_{T}^{\prime})^{\prime}$. We define
$\sigma$–fields: for $t=0,\dots,T$,
$\mathcal{F}_{t}:=\mathcal{F}_{0}\vee\sigma(\bar{y}_{t})$ and
$\mathcal{H}_{t}:=\mathcal{F}_{t}\vee\sigma(\Pi_{s})\vee\sigma(\Gamma_{s})\vee\sigma(s)\vee\sigma(\mathsf{P})$
where for generic sigma fields $\mathcal{O}_{1}$ and $\mathcal{O}_{2}$,
$\mathcal{O}_{1}\vee\mathcal{O}_{2}$ is the minimal $\sigma$–field containing
the $\sigma$–fields $\mathcal{O}_{1}$ and $\mathcal{O}_{2}$. For the
first–order Markov chain, a conditional probability that the regime at time
$t+1$, $s_{t+1}$ equals some particular value conditional on the past regimes
$\bar{s}_{t}$, transition probability matrix $\mathsf{P}$, and initial
information $\mathcal{F}_{0}$ depends only through the most recent regime at
time $t$, $s_{t}$, transition probability matrix $\mathsf{P}$, and initial
information $\mathcal{F}_{0}$, that is,

$p_{s_{t}s_{t+1}}:=\mathbb{P}[s_{t+1}=s_{t+1}|s_{t}=s_{t},\mathsf{P},\mathcal{F}_{0}]=\mathbb{P}\big{[}s_{t+1}=s_{t+1}|\bar{s}_{t}=\bar{s}_{t},\mathsf{P},\mathcal{F}_{0}\big{]}$
(4)

for $t=0,\dots,T-1$, where
$p_{s_{1}}:=p_{s_{0}s_{1}}=\mathbb{P}[s_{1}=s_{1}|\mathsf{P},\mathcal{F}_{0}]$
is an initial probability. A distribution of a residual random vector
$\xi:=(\xi_{1}^{\prime},\dots,\xi_{T}^{\prime})^{\prime}$ is given by

$\xi=(\xi_{1}^{\prime},\dots,\xi_{T}^{\prime})^{\prime}~{}|~{}\mathcal{H}_{0}\sim\mathcal{N}(0,\Sigma_{s}),$
(5)

where $\Sigma_{s}:=\text{diag}\\{\Sigma_{s_{1}},\dots,\Sigma_{s_{T}}\\}$ is a
block diagonal matrix.

To remove duplicates in the random coefficient matrix $(\Pi_{s},\Gamma_{s})$,
for a generic regime vector with length $k$, $o=(o_{1},\dots,o_{k})^{\prime}$,
we define sets

$\mathcal{A}_{\bar{o}_{t}}:=\mathcal{A}_{\bar{o}_{t-1}}\cup\big{\\{}o_{t}\in\\{o_{1},\dots,o_{k}\\}\big{|}o_{t}\not\in\mathcal{A}_{\bar{o}_{t-1}}\big{\\}},~{}~{}~{}t=1,\dots,k,$
(6)

where for $t=1,\dots,k$, $o_{t}\in\\{1,\dots,N\\}$ and an initial set is empty
set, i.e., $\mathcal{A}_{\bar{o}_{0}}=\O$. The final set
$\mathcal{A}_{o}=\mathcal{A}_{\bar{o}_{k}}$ consists of different regimes in
regime vector $o=\bar{o}_{k}$ and $|\mathcal{A}_{o}|$ represents a number of
different regimes in the regime vector $o$.

Let us assume that elements of sets $\mathcal{A}_{s}$,
$\mathcal{A}_{\bar{s}_{t}}$, $\mathcal{A}_{\bar{s}_{t}^{c}}$, intersection set
of the sets $\mathcal{A}_{\bar{s}_{t}}$ and $\mathcal{A}_{\bar{s}_{t}^{c}}$,
and difference sets between the sets $\mathcal{A}_{\bar{s}_{t}^{c}}$ and
$\mathcal{A}_{\bar{s}_{t}}$ are given by
$\mathcal{A}_{s}=\\{\hat{s}_{1},\dots,\hat{s}_{r_{\hat{s}}}\\}$,
$\mathcal{A}_{\bar{s}_{t}}=\\{\alpha_{1},\dots,\alpha_{r_{\alpha}}\\}$,
$\mathcal{A}_{\bar{s}_{t}^{c}}=\\{\beta_{1},\dots,\beta_{r_{\beta}}\\}$,
$\mathcal{A}_{\bar{s}_{t}}\cap\mathcal{A}_{\bar{s}_{t}^{c}}=\\{\gamma_{1},\dots,\gamma_{r_{\gamma}}\\}$,
$\mathcal{A}_{\bar{s}_{t}^{c}}\backslash\mathcal{A}_{\bar{s}_{t}}=\\{\delta_{1},\dots,\delta_{r_{\delta}}\\}$,
and
$\mathcal{A}_{\bar{s}_{t}}\backslash\mathcal{A}_{\bar{s}_{t}^{c}}=\\{\epsilon_{1},\dots,\epsilon_{r_{\epsilon}}\\}$,
respectively, where $r_{\hat{s}}:=|\mathcal{A}_{s}|$,
$r_{\alpha}:=|\mathcal{A}_{\bar{s}_{t}}|$,
$r_{\beta}:=|\mathcal{A}_{\bar{s}_{t}^{c}}|$,
$r_{\gamma}:=|\mathcal{A}_{\bar{s}_{t}}\cap\mathcal{A}_{\bar{s}_{t}^{c}}|$,
$r_{\delta}:=|\mathcal{A}_{\bar{s}_{t}^{c}}\backslash\mathcal{A}_{\bar{s}_{t}}|$,
and
$r_{\epsilon}:=|\mathcal{A}_{\bar{s}_{t}}\backslash\mathcal{A}_{\bar{s}_{t}^{c}}|$
are numbers of elements of the sets, respectively. Note that

$\mathcal{A}_{\bar{s}_{t}}=(\mathcal{A}_{\bar{s}_{t}}\backslash\mathcal{A}_{\bar{s}_{t}^{c}})\cup(\mathcal{A}_{\bar{s}_{t}}\cap\mathcal{A}_{\bar{s}_{t}^{c}}),$
(7)
$\mathcal{A}_{\bar{s}_{t}^{c}}=(\mathcal{A}_{\bar{s}_{t}}\cap\mathcal{A}_{\bar{s}_{t}^{c}})\cup(\mathcal{A}_{\bar{s}_{t}^{c}}\backslash\mathcal{A}_{\bar{s}_{t}}),$
(8)

and

$\mathcal{A}_{s}=\mathcal{A}_{\bar{s}_{t}}\cup\mathcal{A}_{\bar{s}_{t}^{c}}=(\mathcal{A}_{\bar{s}_{t}}\backslash\mathcal{A}_{\bar{s}_{t}^{c}})\cup\mathcal{A}_{\bar{s}_{t}^{c}}=\mathcal{A}_{\bar{s}_{t}}\cup(\mathcal{A}_{\bar{s}_{t}^{c}}\backslash\mathcal{A}_{\bar{s}_{t}})$
(9)

and intersection sets of the sets of right hand sides of equations (7) and
(8), and (9) are empty sets. We introduce the following regime vectors:
$\hat{s}:=(\hat{s}_{1},\dots,\hat{s}_{r_{\hat{s}}})^{\prime}$ is an
$(r_{\hat{s}}\times 1)$ vector,
$\alpha:=(\alpha_{1},\dots,\alpha_{r_{\alpha}})^{\prime}$ is an
$(r_{\alpha}\times 1)$ vector,
$\beta=(\beta_{1},\dots,\beta_{r_{\beta}})^{\prime}$ is an $(r_{\beta}\times
1)$ vector, $\gamma=(\gamma_{1},\dots,\gamma_{r_{\gamma}})^{\prime}$ is an
$(r_{\gamma}\times 1)$ vector,
$\delta=(\delta_{1},\dots,\delta_{r_{\delta}})^{\prime}$ is an
$(r_{\delta}\times 1)$ vector, and
$\epsilon=(\epsilon_{1},\dots,\epsilon_{r_{\epsilon}})^{\prime}$ is an
$(r_{\epsilon}\times 1)$ vector. For the regime vector
$a=(a_{1},\dots,a_{r_{a}})^{\prime}\in\\{\hat{s},\alpha,\beta,\gamma,\delta,\epsilon\\}$,
we also introduce duplication removed random coefficient matrices, whose block
matrices are different: $\Pi_{a}=[\Pi_{a_{1}}:\dots:\Pi_{a_{r_{a}}}]$ is an
$(n\times[dr_{a}])$ matrix,
$\Gamma_{a}=[\Gamma_{a_{1}}:\dots:\Gamma_{a_{r_{a}}}]$ is an
$(n_{*}\times[d_{*}r_{a}])$ matrix, and $(\Pi_{a},\Gamma_{a})$.

We assume that for given duplication removed regime vector $\hat{s}$ and
initial information $\mathcal{F}_{0}$, the coefficient matrices
$(\Pi_{\hat{s}_{1}},\Gamma_{\hat{s}_{1}}),\dots,(\Pi_{\hat{s}_{r_{\hat{s}}}},\Gamma_{\hat{s}_{r_{\hat{s}}}})$
are independent. Under the last assumption, a joint density function of the
random coefficient random matrix $(\Pi_{\hat{s}},\Gamma_{\hat{s}})$ is
represented by

$f\big{(}\Pi_{\hat{s}},\Gamma_{\hat{s}}\big{|}\hat{s},\mathcal{F}_{0}\big{)}=\prod_{t=1}^{r_{\hat{s}}}f\big{(}\Pi_{\hat{s}_{t}},\Gamma_{\hat{s}_{t}}\big{|}\hat{s}_{t},\mathcal{F}_{0}\big{)},$
(10)

where for a generic random vector $X$, we denote its density function by
$f(X)$. Throughout the paper we fix $t=1,\dots,T-1$. For the regime vectors
$\alpha$ and $\beta$, the above joint density function can be written by

$f\big{(}\Pi_{\hat{s}},\Gamma_{\hat{s}}\big{|}\hat{s},\mathcal{F}_{0}\big{)}=f\big{(}\Pi_{\alpha},\Gamma_{\alpha}\big{|}\alpha,\mathcal{F}_{0}\big{)}f_{*}\big{(}\Pi_{\delta},\Gamma_{\delta}\big{|}\delta,\mathcal{F}_{0}\big{)}$
(11)

where the density function
$f_{*}\big{(}\Pi_{\delta},\Gamma_{\delta}\big{|}\delta,\mathcal{F}_{0}\big{)}$
equals

$f_{*}\big{(}\Pi_{\delta},\Gamma_{\delta}\big{|}\delta,\mathcal{F}_{0}\big{)}:=\begin{cases}f\big{(}\Pi_{\delta},\Gamma_{\delta}\big{|}\delta,\mathcal{F}_{0}\big{)},&\text{if}~{}~{}~{}r_{\delta}\neq
0,\\\ 1,&\text{if}~{}~{}~{}r_{\delta}=0.\end{cases}$ (12)

Then, the following Proposition, which is useful for Monte–Carlo simulation
holds, see below.

######  Proposition 1.

Conditional on initial information $\mathcal{F}_{0}$, a joint density function
of the random vectors $\bar{y}_{t}$ and $s$ and random matrices
$\Pi_{\hat{s}}$, $\Gamma_{\hat{s}}$, and $\mathsf{P}$ is given by

$\displaystyle
f(\bar{y}_{t},\Pi_{\hat{s}},\Gamma_{\hat{s}},s,\mathsf{P}|\mathcal{F}_{0})=f\big{(}\bar{y}_{t},\Pi_{\alpha},\Gamma_{\alpha},\bar{s}_{t}\big{|}\mathcal{F}_{0}\big{)}f_{*}\big{(}\Pi_{\delta},\Gamma_{\delta}\big{|}\delta,\mathcal{F}_{0}\big{)}f(s,\mathsf{P}|\mathcal{F}_{0})\big{/}f(\bar{s}_{t}|\mathcal{F}_{0}).$
(13)

In particular, the following relationships holds

$f(\bar{s}_{t}^{c}|\bar{s}_{t},\mathcal{F}_{t})=f(\bar{s}_{t}^{c}|\bar{s}_{t},\mathcal{F}_{0}),$
(14)
$f\big{(}\Pi_{\delta},\Gamma_{\delta}\big{|}\Pi_{\alpha},\Gamma_{\alpha},s,\mathcal{F}_{t}\big{)}=f_{*}\big{(}\Pi_{\delta},\Gamma_{\delta}\big{|}\delta,\mathcal{F}_{0}\big{)},$
(15)
$f\big{(}\Pi_{\beta},\Gamma_{\beta}\big{|}s,\mathcal{F}_{t}\big{)}=f\big{(}\Pi_{\gamma},\Gamma_{\gamma}\big{|}\bar{s}_{t},\mathcal{F}_{t}\big{)}f_{*}\big{(}\Pi_{\delta},\Gamma_{\delta}\big{|}\delta,\mathcal{F}_{0}\big{)},$
(16) $\displaystyle
f(\mathsf{P}|\bar{s}_{t},\mathcal{F}_{t})=f(\mathsf{P}|\bar{s}_{t},\mathcal{F}_{0}),$
(17) $\displaystyle
f(\bar{s}_{t}^{c}|\bar{s}_{t},\mathsf{P},\mathcal{F}_{t})=f(\bar{s}_{t}^{c}|\bar{s}_{t},\mathsf{P},\mathcal{F}_{0}),$
(18)

and

$f(\Pi_{\hat{s}},\Gamma_{\hat{s}}|s,\mathsf{P},\mathcal{F}_{t})=f(\Pi_{\hat{s}},\Gamma_{\hat{s}}|s,\mathcal{F}_{t}).$
(19)

###### Proof.

Let us consider a joint density function of the random vectors $\bar{y}_{t}$
and $s$ and random matrices $\Pi_{\hat{s}}$, $\Gamma_{\hat{s}}$, and
$\mathsf{P}$ for given initial information $\mathcal{F}_{0}$. According to the
conditional probability formula, one gets that

$f(\bar{y}_{t},\Pi_{\hat{s}},\Gamma_{\hat{s}},s,\mathsf{P}|\mathcal{F}_{0})=f(\bar{y}_{t}|\Pi_{\hat{s}},\Gamma_{\hat{s}},s,\mathsf{P},\mathcal{F}_{0})f(\Pi_{\hat{s}},\Gamma_{\hat{s}}|s,\mathsf{P},\mathcal{F}_{0})f(s|\mathsf{P},\mathcal{F}_{0})f(\mathsf{P}|\mathcal{F}_{0}).$
(20)

Since the random vector $\bar{y}_{t}$ depends on $\Pi_{\alpha}$,
$\Gamma_{\alpha}$, $\bar{s}_{t}$, the first joint density function of the
right–hand side of the above equation equals
$f(\bar{y}_{t}|\Pi_{\alpha},\Gamma_{\alpha},\bar{s}_{t},\mathcal{F}_{0})$. As
the random coefficient matrix $(\Pi_{\hat{s}},\Gamma_{\hat{s}})$ is
independent of the transition probability matrix $\mathsf{P}$ conditional on
$s$ and $\mathcal{F}_{0}$, the second joint density function of the right–hand
side can be represented by equation (11). According to the Markov property
(4), the third joint density function equals
$f(\bar{s}_{t}|\mathsf{P},\mathcal{F}_{0})f(\bar{s}_{t}^{c}|\bar{s}_{t},\mathsf{P},\mathcal{F}_{0})$.
Consequently, one obtains equation (13). If we integrate equation (13) by
$\mathsf{P}$, then one finds that

$f(\bar{y}_{t},\Pi_{\hat{s}},\Gamma_{\hat{s}},s|\mathcal{F}_{0})=f\big{(}\bar{y}_{t},\Pi_{\alpha},\Gamma_{\alpha},\bar{s}_{t}\big{|}\mathcal{F}_{0}\big{)}f_{*}\big{(}\Pi_{\delta},\Gamma_{\delta}\big{|}\delta,\mathcal{F}_{0}\big{)}f(\bar{s}_{t}^{c}|\bar{s}_{t},\mathcal{F}_{0}).$
(21)

By integrating the above equation by the random matrix
$(\Pi_{\hat{s}},\Gamma_{\hat{s}})$, one obtains equation (14). Since
$f(\bar{y}_{t}|\Pi_{\alpha},\Gamma_{\alpha},\bar{s}_{t},\bar{s}_{t}^{c},\mathcal{F}_{0})=f(\bar{y}_{t}|\Pi_{\alpha},\Gamma_{\alpha},\bar{s}_{t},\mathcal{F}_{0})$
and
$f(\Pi_{\alpha},\Gamma_{\alpha}|\bar{s}_{t},\bar{s}_{t}^{c},\mathcal{F}_{0})=f(\Pi_{\alpha},\Gamma_{\alpha}|\bar{s}_{t},\mathcal{F}_{0})$,
we have that

$\displaystyle f(\bar{y}_{t},\bar{s}_{t}^{c}|\bar{s}_{t},\mathcal{F}_{0})$
$\displaystyle=$
$\displaystyle\int_{\Pi_{\alpha},\Gamma_{\alpha}}f(\bar{y}_{t}|\Pi_{\alpha},\Gamma_{\alpha},\bar{s}_{t},\bar{s}_{t}^{c},\mathcal{F}_{0})f(\Pi_{\alpha},\Gamma_{\alpha}|\bar{s}_{t},\bar{s}_{t}^{c},\mathcal{F}_{0})f(\bar{s}_{t}^{c}|\bar{s}_{t}\mathcal{F}_{0})d\Pi_{\alpha}d\Gamma_{\alpha}$
(22) $\displaystyle=$ $\displaystyle
f(\bar{y}_{t}|\bar{s}_{t},\mathcal{F}_{0})f(\bar{s}_{t}^{c}|\bar{s}_{t},\mathcal{F}_{0}).$

Thus, conditional on $\bar{s}_{t}$ and $\mathcal{F}_{0}$, the random vectors
$\bar{y}_{t}$ and $\bar{s}_{t}^{c}$ are independent. Consequently, it holds
that

$f(\bar{y}_{t},s|\mathcal{F}_{0})=f(\bar{y}_{t},\bar{s}_{t}|\mathcal{F}_{0})f(\bar{s}_{t}^{c}|\bar{s}_{t},\mathcal{F}_{0}).$
(23)

If we divide equation (21) by the above equation, then one obtains

$f\big{(}\Pi_{\hat{s}},\Gamma_{\hat{s}}\big{|}s,\mathcal{F}_{t}\big{)}=f\big{(}\Pi_{\alpha},\Gamma_{\alpha}\big{|}\bar{s}_{t},\mathcal{F}_{t}\big{)}f_{*}\big{(}\Pi_{\delta},\Gamma_{\delta}\big{|}\delta,\mathcal{F}_{0}\big{)}.$
(24)

Since
$f(\bar{y}_{t},\Pi_{\alpha},\Gamma_{\alpha}|s,\mathcal{F}_{0})=f(\bar{y}_{t},\Pi_{\alpha},\Gamma_{\alpha}|\bar{s}_{t},\mathcal{F}_{0})$
and conditional on $\bar{s}_{t}$ and $\mathcal{F}_{0}$, the random vectors
$\bar{y}_{t}$ and $\bar{s}_{t}^{c}$ are independent, we have
$f\big{(}\Pi_{\alpha},\Gamma_{\alpha}\big{|}s,\mathcal{F}_{t}\big{)}=f\big{(}\Pi_{\alpha},\Gamma_{\alpha}\big{|}\bar{s}_{t},\mathcal{F}_{t}\big{)}$.
Consequently, due to the conditional probability formula, we have that

$f\big{(}\Pi_{\hat{s}},\Gamma_{\hat{s}}\big{|}s,\mathcal{F}_{t}\big{)}=f\big{(}\Pi_{\alpha},\Gamma_{\alpha}\big{|}\bar{s}_{t},\mathcal{F}_{t}\big{)}f\big{(}\Pi_{\delta},\Gamma_{\delta}\big{|}\Pi_{\alpha},\Gamma_{\alpha},s,\mathcal{F}_{0}\big{)}.$
(25)

Thus, equating equation (24) with the above equation, we get equation (15). If
we integrate equation (24) by $(\Pi_{\epsilon},\Gamma_{\epsilon})$, then by
equation (7) and (9), we obtain equation (16). Integrating equation (13) by
$(\Pi_{\hat{s}},\Gamma_{\hat{s}})$, one obtains that

$\displaystyle
f(\bar{y}_{t},s,\mathsf{P}|\mathcal{F}_{0})=f\big{(}\bar{y}_{t},\bar{s}_{t}\big{|}\mathcal{F}_{0}\big{)}f(s,\mathsf{P}|\mathcal{F}_{0})\big{/}f(\bar{s}_{t}|\mathcal{F}_{0}).$
(26)

Consequently if we integrate equation (26) by $\bar{s}_{t}^{c}$ and divide
$f(\bar{y}_{t},\bar{s}_{t}|\mathcal{F}_{0})$, we get equation (17). Thus,
conditional on $\bar{s}_{t}$ and $\mathcal{F}_{0}$, $\bar{y}_{t}$ and
$\mathsf{P}$ are independent. Also, it follows from equation (26) that

$\displaystyle
f(\bar{s}_{t}^{c},\mathsf{P}|\bar{s}_{t},\mathcal{F}_{t})=f(s,\mathsf{P}|\mathcal{F}_{0})\big{/}f(\bar{s}_{t}|\mathcal{F}_{0}).$
(27)

By dividing the above equation by equation (17), we obtain equation (18). To
prove equation (19), let us consider a conditional density function
$f(\Pi_{\hat{s}},\Gamma_{\hat{s}},\mathsf{P}|s,\mathcal{F}_{t})$. Equation
(20) can be written by

$f(\bar{y}_{t},\Pi_{\hat{s}},\Gamma_{\hat{s}},s,\mathsf{P}|\mathcal{F}_{0})=f(\bar{y}_{t},\Pi_{\hat{s}},\Gamma_{\hat{s}}|s,\mathcal{F}_{0})f(s,\mathsf{P}|\mathcal{F}_{0}).$
(28)

On the other hand, by the conditional probability formula, it holds

$f(\bar{y}_{t},\Pi_{\hat{s}},\Gamma_{\hat{s}},s,\mathsf{P}|\mathcal{F}_{0})=f(\bar{y}_{t},\Pi_{\hat{s}},\Gamma_{\hat{s}}|s,\mathsf{P},\mathcal{F}_{0})f(s,\mathsf{P}|\mathcal{F}_{0}).$
(29)

Therefore, one conclude that conditional on $s$ and $\mathcal{F}_{0}$, the
random matrix $(\bar{y}_{t},\Pi_{\hat{s}},\Gamma_{\hat{s}})$ and transition
probability matrix $\mathsf{P}$ are independent. Thus, equation (19) holds.
That completes the proof. ∎

It should be noted that according to the Markov property (4), it follows from
equation (18) that the assumption for a Markov chain in the book of Hamilton94
always holds, namely,

$\displaystyle
f(s_{t+1}|\bar{s}_{t},\mathsf{P},\mathcal{F}_{t})=f(s_{t+1}|s_{t},\mathsf{P},\mathcal{F}_{0}).$
(30)

##  3 Special Case of Bayesian MS–VAR($p$) process

In this section, we consider a special case of the Bayesian MS–VAR($p$)
process. The Bayesian MS–VAR($p$) process can be written by the following
equation

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

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

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

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

$\xi_{t}~{}|~{}\Sigma_{s_{t}},s_{t},\mathcal{F}_{0}\sim\mathcal{N}\Big{(}0,\Sigma_{s_{t}}\Big{)},$
(33)

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

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

where $\pi_{0,s_{t}}$ is an $(nd\times 1)$ prior hyperparameter vector at
regime $s_{t}$ and $\Lambda_{0,s_{t}}$ is a symmetric positive definite
$(d\times d)$ prior hyperparameter matrix at regime $s_{t}$.

###  3.1 Distributions

For the regime vector $\bar{s}_{t}$ and regime $\alpha_{k}$, we define sets

$S_{t,\alpha_{k}}:=\big{\\{}u\in\\{1,\dots,t\\}\big{|}s_{u}=\alpha_{k},~{}u=1,\dots,t\big{\\}},~{}~{}~{}k=1,\dots,r_{\alpha}.$
(35)

For $k=1,\dots,r_{\alpha}$, the set $S_{t,\alpha_{k}}$ consists of indexes of
regimes in the regime vector $\bar{s}_{t}$ that equal the regime $\alpha_{k}$.
Let us suppose that $q_{t,\alpha_{k}}:=|S_{t,\alpha_{k}}|$ is a number of
regimes in the regime vector $\bar{s}_{t}$ that equal the regime $\alpha_{k}$
and elements of the set $S_{t,\alpha_{k}}$ are given by

$S_{t,\alpha_{k}}=\\{k_{t,1},\dots,k_{t,q_{t,\alpha_{k}}}\\},~{}~{}~{}k=1,\dots,r_{\alpha}.$
(36)

Further, we define indexes

$o_{t}:=\big{\\{}k\in\\{1,\dots,r_{\alpha}\\}\big{|}s_{t}=\alpha_{k},~{}k=1,\dots,r_{\alpha}\big{\\}},~{}~{}~{}t=1,\dots,t.$
(37)

The index $o_{t}$ represents a position of the regime $s_{t}$ in the regime
vector $\alpha$. Let
$\pi_{\alpha}:=\big{(}\pi_{\alpha_{1}}^{\prime},\dots,\pi_{\alpha_{r_{\alpha}}}^{\prime}\big{)}^{\prime}$
be an $([ndr_{\alpha}]\times 1)$ random coefficient vector, whose sub–vectors
are different and which corresponds to the regime vector $\bar{s}_{t}$,
$y_{t,\alpha_{k}}:=\Big{(}y_{k_{t,1}}^{\prime},\dots,y_{k_{t,q_{t,\alpha_{k}}}}^{\prime}\Big{)}^{\prime}$
be an $([nq_{t,\alpha_{k}}]\times 1)$ vector of endogenous variables,
corresponding to the regime $\alpha_{k}$, and
$\mathsf{Y}_{t,\alpha_{k}}^{\circ}:=\big{[}\mathsf{Y}_{k_{t,1}}:\dots:\mathsf{Y}_{k_{t,q_{t,\alpha_{k}}}}\big{]}$
be a $(d\times q_{t,\alpha_{k}})$ matrix of exogenous and endogenous
variables, corresponding to the regime $\alpha_{k}$. By using a $(t\times
r_{\alpha})$ matrix $D_{\alpha}:=[j_{o_{1}}:\dots:j_{o_{t}}]^{\prime}$ one can
revive the vector $\pi_{\bar{s}_{t}}:=\text{vec}(\Pi_{\bar{s}_{t}})$ from the
vector $\pi_{\alpha}$, that is, $\pi_{\bar{s}_{t}}=(D_{\alpha}\otimes
I_{nd})\pi_{\alpha}$, where $j_{o}$ is an $(r_{\alpha}\times 1)$ unit vector,
whose $o$–th element equals one and others zero.

####  3.1.1 Conditional Distributions

It follows from equations (33) and (34) that distributions of the random
vectors $\bar{\xi}_{t}$ and $\pi_{\alpha}$ are obtained by

$\bar{\xi}_{t}~{}|~{}\Sigma_{\alpha},\bar{s}_{t},\mathcal{F}_{0}\sim\mathcal{N}\big{(}0,\Sigma_{\bar{s}_{t}}\big{)},$
(38)

and

$\pi_{\alpha}~{}|~{}\Sigma_{\alpha},\alpha,\mathcal{F}_{0}\sim\mathcal{N}\Big{(}\pi_{0,\alpha},\Sigma_{\pi_{\alpha}}\Big{)},$
(39)

respectively, where
$\Sigma_{\bar{s}_{t}}:=\text{diag}\\{\Sigma_{s_{1}},\dots,\Sigma_{s_{t}}\\}$
is an ($[nt]\times[nt]$) block diagonal matrix, corresponding to the regime
vector $\bar{s}_{t}$ and
$\Sigma_{\alpha}:=\big{[}\Sigma_{\alpha_{1}}:\dots:\Sigma_{\alpha_{r_{\alpha}}}\big{]}^{\prime}$
is an $([nr_{\alpha}]\times n)$ matrix,
$\pi_{0,\alpha}:=\big{(}\pi_{0,\alpha_{1}}^{\prime},\dots,\pi_{0,\alpha_{r_{\alpha}}}^{\prime}\big{)}^{\prime}$
is an $([ndr_{\alpha}]\times 1)$ prior hyperparameter vector, and
$\Sigma_{\pi_{\alpha}}:=\text{diag}\big{\\{}\Lambda_{0,\alpha_{1}}\otimes\Sigma_{\alpha_{1}},\dots,\Lambda_{0,\alpha_{r_{\alpha}}}\otimes\Sigma_{\alpha_{r_{\alpha}}}\big{\\}}$
is an $([ndr_{\alpha}]\times[ndr_{\alpha}])$ block diagonal matrix, all of
which correspond to the duplication removed regime vector $\alpha$. A
connection between the random matrices $\Sigma_{\bar{s}_{t}}$ and
$\Sigma_{\alpha}$ is

$\Sigma_{\bar{s}_{t}}=\text{diag}\big{\\{}\big{(}(D_{\alpha}\otimes
I_{n})\Sigma_{\alpha}\big{)}_{1},\dots,\big{(}(D_{\alpha}\otimes
I_{n})\Sigma_{\alpha}\big{)}_{t}\big{\\}},$ (40)

where the matrix $\big{(}(D_{\alpha}\otimes I_{n})\Sigma_{\alpha}\big{)}_{j}$
equals $j$–th block matrix of the matrix $(D_{\alpha}\otimes
I_{n})\Sigma_{\alpha}$. On the other hand, by following Battulga24a, a
distribution of $([nT]\times 1)$ the random vector
$y=(y_{1}^{\prime},\dots,y_{T}^{\prime})^{\prime}$ is given by

$y~{}\big{|}~{}\pi_{\hat{s}},\Sigma_{\hat{s}},s,\mathcal{F}_{0}\sim\mathcal{N}\Big{(}\Psi_{s}^{-1}\varphi_{s},\Psi_{s}^{-1}\Sigma_{s}(\Psi_{s}^{-1})^{\prime}\Big{)},$
(41)

where the matrix $\Psi_{s}$ and the vector $\delta_{s}$ are

$\Psi_{s}:=\begin{bmatrix}I_{n}&0&\dots&0&\dots&0&0\\\
-A_{1,s_{2}}&I_{n}&\dots&0&\dots&0&0\\\
\vdots&\vdots&\dots&\vdots&\dots&\vdots&\vdots\\\
0&0&\dots&-A_{p-1,s_{T-1}}&\dots&I_{n}&0\\\
0&0&\dots&-A_{p,s_{T}}&\dots&-A_{1,s_{T}}&I_{n}\end{bmatrix}$ (42)

and

$\varphi_{s}:=\begin{bmatrix}A_{0,s_{1}}\psi_{1}+A_{1,s_{1}}y_{0}+\dots+A_{p,s_{1}}y_{1-p}\\\
A_{0,s_{2}}\psi_{2}+A_{2,s_{2}}y_{0}+\dots+A_{p,s_{2}}y_{2-p}\\\ \vdots\\\
A_{0,s_{T-1}}\psi_{T-1}\\\ A_{0,s_{T}}\psi_{T}\end{bmatrix},$ (43)

respectively. To price default–free options, Battulga24a used the conditional
distribution of the random vector $y$. For a generic vector
$o=(o_{1}^{\prime},\dots,o_{n}^{\prime})^{\prime}$ with ($m\times 1$) vector
$o_{i}$, we introduce an $(m\times n)$ matrix notation
$o^{\circ}:=[o_{1}:\dots:o_{n}]$. Then, the following Proposition holds.

######  Proposition 2.

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

$\displaystyle
f(\bar{y}_{t}|\bar{s}_{t},\mathcal{F}_{0})=\frac{1}{\pi^{nt/2}}\prod_{k=1}^{r_{\alpha}}\frac{|\Lambda_{0,\alpha_{k}}^{-1}|^{n/2}\Gamma_{n}\big{(}(\nu_{0,\alpha_{k}}+q_{t,\alpha_{k}})/2\big{)}|V_{0,\alpha_{k}}|^{\nu_{0,\alpha_{k}}/2}}{|\Lambda_{0,\alpha_{k}|t}|^{n/2}\Gamma_{n}(\nu_{0,\alpha_{k}}/2)\big{|}\bar{B}_{t,\alpha_{k}}+V_{0,\alpha_{k}}\big{|}^{(\nu_{0,\alpha_{k}}+q_{t,\alpha_{k}})/2}},$
(44)

where $\Gamma_{n}(\cdot)$ is the multivariate gamma function,
$\Lambda_{0,\alpha_{k}|t}:=\mathsf{Y}_{t,\alpha_{k}}^{\circ}(\mathsf{Y}_{t,\alpha_{k}}^{\circ})^{\prime}+\Lambda_{0,\alpha_{k}}^{-1}$
is a $(d\times d)$ matrix, and $\bar{B}_{t,\alpha_{k}}$ is an $(n\times n)$
positive semi–definite matrix and equals

$\displaystyle\bar{B}_{t,\alpha_{k}}$ $\displaystyle:=$ $\displaystyle
y_{t,\alpha_{k}}^{\circ}(y_{t,\alpha_{k}}^{\circ})^{\prime}+\pi_{0,\alpha_{k}}^{\circ}\Lambda_{0,\alpha_{k}}^{-1}(\pi_{0,\alpha_{k}}^{\circ})^{\prime}-C_{t,\alpha_{k}}\Lambda_{0,\alpha_{k}|t}C_{t,\alpha_{k}}^{\prime}$
$\displaystyle=$
$\displaystyle\big{(}y_{t,\alpha_{k}}^{\circ}-\pi_{0,\alpha_{k}}^{\circ}\mathsf{Y}_{t,\alpha_{k}}^{\circ}\big{)}\big{(}I_{q_{t,\alpha_{k}}}+(\mathsf{Y}_{t,\alpha_{k}}^{\circ})^{\prime}\Lambda_{0,\alpha_{k}}\mathsf{Y}_{t,\alpha_{k}}^{\circ}\big{)}^{-1}\big{(}y_{t,\alpha_{k}}^{\circ}-\pi_{0,\alpha_{k}}^{\circ}\mathsf{Y}_{t,\alpha_{k}}^{\circ}\big{)}^{\prime}$

with $(n\times d)$ matrix
$C_{t,\alpha_{k}}:=\big{(}y_{t,\alpha_{k}}^{\circ}(\mathsf{Y}_{t,\alpha_{k}}^{\circ})^{\prime}+\pi_{0,\alpha_{k}}^{\circ}\Lambda_{0,\alpha_{k}}^{-1}\big{)}\Lambda_{0,\alpha_{k}|t}^{-1}$.
Second, conditional on the random covariance matrix $\Sigma_{\alpha}$, regime
vector $\bar{s}_{t}$, and information $\mathcal{F}_{t}$, a joint density
function of the random coefficient vector $\pi_{\alpha}$ is given by

$\displaystyle f(\pi_{\alpha}|\Sigma_{\alpha},\bar{s}_{t},\mathcal{F}_{t})$
$\displaystyle=\frac{1}{(2\pi)^{ndr_{\alpha}/2}\prod_{k=1}^{r_{\alpha}}|A_{\alpha_{k}}|^{1/2}}\exp\bigg{\\{}-\frac{1}{2}\sum_{k=1}^{r_{\alpha}}\Big{(}\pi_{\alpha_{k}}-A_{\alpha_{k}}b_{\alpha_{k}}\Big{)}^{\prime}A_{\alpha_{k}}^{-1}\Big{(}\pi_{\alpha_{k}}-A_{\alpha_{k}}b_{\alpha_{k}}\Big{)}\bigg{\\}},$
(46)

where for $k=1,\dots,r_{\alpha}$,
$A_{\alpha_{k}}:=\big{(}\Lambda_{0,\alpha_{k}|t}^{-1}\otimes\Sigma_{\alpha_{k}}\big{)}$
is an $([nd]\times[nd])$ matrix and
$b_{\alpha_{k}}:=\big{(}\mathsf{Y}_{t,\alpha_{k}}^{\circ}\otimes\Sigma_{\alpha_{k}}^{-1}\big{)}y_{t,\alpha_{k}}+\big{(}\Lambda_{0,\alpha_{k}}^{-1}\otimes\Sigma_{\alpha_{k}}^{-1}\big{)}\pi_{0,\alpha_{k}}$
is an $([nd]\times 1)$ vector. Third, conditional on the regime vector
$\bar{s}_{t}$ and information $\mathcal{F}_{t}$, a joint density function of
the random coefficient matrix $\Sigma_{\alpha}$ is given by

$\displaystyle f(\Sigma_{\alpha}|\bar{s}_{t},\mathcal{F}_{t})$
$\displaystyle=$
$\displaystyle\prod_{k=1}^{r_{\alpha}}\frac{\big{|}\bar{B}_{t,\alpha_{k}}+V_{0,\alpha_{k}}\big{|}^{(\nu_{0,\alpha_{k}}+q_{t,\alpha_{k}})/2}}{\Gamma_{n}\big{(}(\nu_{0,\alpha_{k}}+q_{t,\alpha_{k}})/2\big{)}2^{n(\nu_{0,\alpha_{k}}+q_{t,\alpha_{k}})/2}}|\Sigma_{\alpha_{k}}|^{-(\nu_{0,\alpha_{k}}+q_{t,\alpha_{k}}+n+1)/2}$
(47) $\displaystyle\times$
$\displaystyle\exp\bigg{\\{}-\frac{1}{2}\sum_{k=1}^{r_{\alpha}}\mathrm{tr}\Big{(}\big{(}\bar{B}_{t,\alpha_{k}}+V_{0,\alpha_{k}}\big{)}\Sigma_{\alpha_{k}}^{-1}\Big{)}\bigg{\\}}.$

Finally, conditional on the regime vector $s$ and information
$\mathcal{F}_{t}$, a joint density function of the random coefficient matrix
$\pi_{\beta}^{\circ}$ is given by

$\displaystyle f(\pi_{\beta}^{\circ}|s,\mathcal{F}_{t})$ $\displaystyle=$
$\displaystyle\prod_{k=1}^{r_{\gamma}}\frac{|\Lambda_{0,\gamma_{k}|t}|^{n/2}\big{|}\bar{B}_{t,\gamma_{k}}+V_{0,\gamma_{k}}\big{|}^{-d/2}\Gamma_{n}\big{(}(\nu_{0,\gamma_{k}}+q_{t,\gamma_{k}}+d)/2\big{)}}{\pi^{nd/2}\Gamma_{n}((\nu_{0,\gamma_{k}}+q_{t,\gamma_{k}})/2)}$
$\displaystyle\times$
$\displaystyle\big{|}I_{n}+(\bar{B}_{t,\gamma_{k}}+V_{0,\gamma_{k}})^{-1}(\pi_{\gamma_{k}}^{\circ}-C_{t,\gamma_{k}})\Lambda_{0,\gamma_{k}|t}(\pi_{\gamma_{k}}^{\circ}-C_{t,\gamma_{k}})^{\prime}\big{|}^{-(\nu_{0,\gamma_{k}}+q_{t,\gamma_{k}}+d)/2}$
$\displaystyle\times$
$\displaystyle\prod_{\ell=1}^{r_{\delta}}\frac{|\Lambda_{0,\delta_{\ell}}|^{-n/2}|V_{0,\delta_{\ell}}|^{-d/2}\Gamma_{n}\big{(}(\nu_{0,\delta_{\ell}}+d)/2\big{)}}{\pi^{nd/2}\Gamma_{n}(\nu_{0,\delta_{\ell}}/2)}$
$\displaystyle\times$
$\displaystyle\big{|}I_{n}+V_{0,\delta_{\ell}}^{-1}(\pi_{\delta_{\ell}}^{\circ}-\pi_{0,\delta_{\ell}}^{\circ})\Lambda_{0,\delta_{\ell}}^{-1}(\pi_{\delta_{\ell}}^{\circ}-\pi_{0,\delta_{\ell}}^{\circ})^{\prime}\big{|}^{-(\nu_{0,\delta_{\ell}}+d)/2}.$

###### Proof.

First, since
$(\pi_{\alpha_{1}},\Sigma_{\alpha_{1}}),\dots,(\pi_{r_{\alpha}},\Sigma_{r_{\alpha}})$
are independent for given regime vector $\bar{s}_{t}$ and initial information
$\mathcal{F}_{0}$, observe that conditional density functions of the random
vectors $y$ and $\pi_{\hat{s}}$ and random matrix $\Sigma_{\hat{s}}$ are given
by

$\displaystyle
f(\bar{y}_{t}|\pi_{\alpha},\Sigma_{\alpha},\bar{s}_{t},\mathcal{F}_{0})=\frac{1}{(2\pi)^{nt/2}}\prod_{k=1}^{r_{\alpha}}|\Sigma_{\alpha_{k}}|^{-q_{t,\alpha_{k}}/2}$
(49)
$\displaystyle\times\exp\bigg{\\{}-\frac{1}{2}\sum_{k=1}^{r_{\alpha}}\Big{(}y_{t,\alpha_{k}}-\big{(}(\mathsf{Y}_{t,\alpha_{k}}^{\circ})^{\prime}\otimes
I_{n}\big{)}\pi_{\alpha_{k}}\Big{)}^{\prime}\big{(}I_{q_{t,\alpha_{k}}}\otimes\Sigma_{\alpha_{k}}^{-1}\big{)}\Big{(}y_{t,\alpha_{k}}-\big{(}(\mathsf{Y}_{t,\alpha_{k}}^{\circ})^{\prime}\otimes
I_{n}\big{)}\pi_{\alpha_{k}}\Big{)}\bigg{\\}},$ $\displaystyle
f(\pi_{\alpha}|\Sigma_{\alpha},\bar{s}_{t},\mathcal{F}_{0})$ $\displaystyle=$
$\displaystyle\frac{1}{(2\pi)^{nr_{\alpha}d/2}\prod_{k=1}^{r_{\alpha}}|\Lambda_{0,\alpha_{k}}|^{n/2}}\prod_{k=1}^{r_{\alpha}}|\Sigma_{\alpha_{k}}|^{-d/2}$
$\displaystyle\times$
$\displaystyle\exp\bigg{\\{}-\frac{1}{2}\sum_{k=1}^{r_{\alpha}}\big{(}\pi_{\alpha_{k}}-\pi_{0,\alpha_{k}}\big{)}^{\prime}\big{(}\Lambda_{0,\alpha_{k}}^{-1}\otimes\Sigma_{\alpha_{k}}^{-1}\big{)}\big{(}\pi_{\alpha_{k}}-\pi_{0,\alpha_{k}}\big{)}\bigg{\\}},$

and

$f(\Sigma_{\alpha}|\bar{s}_{t},\mathcal{F}_{0})=\prod_{k=1}^{r_{\alpha}}\frac{|V_{0,\alpha_{k}}|^{\nu_{0,\alpha_{k}}/2}}{\Gamma_{n}(\nu_{0,\alpha_{k}}/2)2^{n\nu_{0,\alpha_{k}}/2}}|\Sigma_{\alpha_{k}}|^{-(\nu_{0,\alpha_{k}}+n+1)/2}\exp\bigg{\\{}-\frac{1}{2}\text{tr}\Big{(}V_{0,\alpha_{k}}\Sigma_{\alpha_{k}}^{-1}\Big{)}\bigg{\\}},$
(51)

respectively. Consequently, by the completing square method, a joint
conditional density function of the random vectors $y$ and $\pi_{\hat{s}}$ is

$\displaystyle
f(\bar{y}_{t},\pi_{\alpha}|\Sigma_{\alpha},\bar{s}_{t},\mathcal{F}_{0})$
$\displaystyle=c_{1}\prod_{k=1}^{r_{\alpha}}|\Sigma_{\alpha_{k}}|^{-(q_{t,\alpha_{k}}+d)/2}\exp\bigg{\\{}-\frac{1}{2}\sum_{k=1}^{r_{\alpha}}\Big{(}\pi_{\alpha_{k}}-A_{\alpha_{k}}b_{\alpha_{k}}\Big{)}^{\prime}A_{\alpha_{k}}^{-1}\Big{(}\pi_{\alpha_{k}}-A_{\alpha_{k}}b_{\alpha_{k}}\Big{)}\bigg{\\}}$
(52)
$\displaystyle\times\exp\bigg{\\{}-\frac{1}{2}\sum_{k=1}^{r_{\alpha}}\Big{(}y_{t,\alpha_{k}}^{\prime}\big{(}I_{q_{t,\alpha_{k}}}\otimes\Sigma_{\alpha_{k}}^{-1}\big{)}y_{t,\alpha_{k}}+\pi_{0,\alpha_{k}}^{\prime}\big{(}\Lambda_{0,\alpha_{k}}^{-1}\otimes\Sigma_{\alpha_{k}}^{-1}\big{)}\pi_{0,\alpha_{k}}-b_{\alpha_{k}}^{\prime}A_{\alpha_{k}}b_{\alpha_{k}}\Big{)}\bigg{\\}},$

where normalizing coefficient equals

$c_{1}:=\frac{1}{(2\pi)^{n(t+r_{\alpha}d)/2}\prod_{k=1}^{r_{\alpha}}|\Lambda_{0,\alpha_{k}}|^{n/2}}.$
(53)

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

$\displaystyle
f(\bar{y}_{t}|\Sigma_{\alpha},\bar{s}_{t},\mathcal{F}_{0})=c_{2}\prod_{k=1}^{r_{\alpha}}|\Sigma_{\alpha_{k}}|^{-q_{t,\alpha_{k}}/2}$
$\displaystyle\times\exp\bigg{\\{}-\frac{1}{2}\sum_{k=1}^{r_{\alpha}}\Big{(}y_{t,\alpha_{k}}^{\prime}\big{(}I_{q_{t,\alpha_{k}}}\otimes\Sigma_{\alpha_{k}}^{-1}\big{)}y_{t,\alpha_{k}}+\pi_{0,\alpha_{k}}^{\prime}\big{(}\Lambda_{0,\alpha_{k}}^{-1}\otimes\Sigma_{\alpha_{k}}^{-1}\big{)}\pi_{0,\alpha_{k}}-b_{\alpha_{k}}^{\prime}A_{\alpha_{k}}b_{\alpha_{k}}\Big{)}\bigg{\\}},$
(54)

where the normalizing coefficient equals

$c_{2}:=\frac{1}{(2\pi)^{nt/2}\prod_{k=1}^{r_{\alpha}}|\Lambda_{0,\alpha_{k}}|^{n/2}|\Lambda_{0,\alpha_{k}|t}|^{n/2}}.$
(55)

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),$ (56)

we find that

$f(\bar{y}_{t}|\Sigma_{\alpha},\bar{s}_{t},\mathcal{F}_{0})=c_{2}\prod_{k=1}^{r_{\alpha}}|\Sigma_{\alpha_{k}}|^{-q_{t,\alpha_{k}}/2}\exp\bigg{\\{}-\frac{1}{2}\sum_{k=1}^{r_{\alpha}}\text{tr}\big{(}\bar{B}_{t,\alpha_{k}}\Sigma_{\alpha_{k}}^{-1}\big{)}\bigg{\\}}.$
(57)

Thus, it follows from equations (51) and (57) that a joint conditional density
of the random vector $y$ and random matrix $\Sigma_{\hat{s}}$ is

$f(\bar{y}_{t},\Sigma_{\alpha}|\bar{s}_{t},\mathcal{F}_{0})=c_{3}\prod_{k=1}^{r_{\alpha}}|\Sigma_{\alpha_{k}}|^{-(q_{t,\alpha_{k}}+\nu_{0,\alpha_{k}}+n+1)/2}\exp\bigg{\\{}-\frac{1}{2}\sum_{k=1}^{r_{\alpha}}\text{tr}\Big{(}\big{(}\bar{B}_{t,\alpha_{k}}+V_{0,\alpha_{k}}\big{)}\Sigma_{\alpha_{k}}^{-1}\Big{)}\bigg{\\}},$
(58)

where the normalizing coefficient equals

$c_{3}:=\frac{1}{(2\pi)^{nt/2}}\prod_{k=1}^{r_{\alpha}}\frac{|V_{0,\alpha_{k}}|^{\nu_{0,\alpha_{k}}/2}}{|\Lambda_{0,\alpha_{k}}|^{n/2}|\Lambda_{0,\alpha_{k}|t}|^{n/2}\Gamma_{n}(\nu_{0,\alpha_{k}}/2)2^{n\nu_{0,\alpha_{k}}/2}}.$
(59)

Consequently, a prior density of the random vector $y$ is given by

$\displaystyle f(\bar{y}_{t}|\bar{s}_{t},\mathcal{F}_{0})$ $\displaystyle=$
$\displaystyle\int_{\Sigma_{\hat{s}_{1}},\dots,\Sigma_{\hat{s}_{r_{\alpha}}}>0}f(\bar{y}_{t},\hat{\Sigma}_{\bar{s}_{t}}|\bar{s}_{t},\mathcal{F}_{0})d\Sigma_{\hat{s}_{1}}\dots
d\Sigma_{\hat{s}_{r_{\alpha}}}$ (60) $\displaystyle=$ $\displaystyle
c_{3}\prod_{k=1}^{r_{\alpha}}\frac{\Gamma_{n}\big{(}(\nu_{0,\alpha_{k}}+q_{t,\alpha_{k}})/2\big{)}2^{n(\nu_{0,\alpha_{k}}+q_{t,\alpha_{k}})/2}}{\big{|}\bar{B}_{t,\alpha_{k}}+V_{0,\alpha_{k}}\big{|}^{(\nu_{0,\alpha_{k}}+q_{t,\alpha_{k}})/2}}.$

If we divide equations (3.1.1) and (58) by equations (3.1.1) and (60),
respectively, then one obtains equations (2) and (47). According to equation
(16), we have that

$f(\pi_{\beta},\Sigma_{\beta}|s,\mathcal{F}_{t})=f(\pi_{\gamma},\Sigma_{\gamma}|\bar{s}_{t},\mathcal{F}_{t})f_{*}(\pi_{\delta},\Sigma_{\delta}|\delta,\mathcal{F}_{0}).$
(61)

We consider the first joint conditional density of the right–hand side of the
above equation. By integrating a product of density functions (2) and (47) by
$(\pi_{\epsilon},\Sigma_{\epsilon})$ and taking account that

$\sum_{k=1}^{r_{\gamma}}\Big{(}\pi_{\gamma_{k}}-A_{\gamma_{k}}b_{\gamma_{k}}\Big{)}^{\prime}A_{\gamma_{k}}^{-1}\Big{(}\pi_{\gamma_{k}}-A_{\gamma_{k}}b_{\gamma_{k}}\Big{)}=\sum_{k=1}^{r_{\gamma}}\mathrm{tr}\Big{(}(\pi_{\gamma_{k}}^{\circ}-C_{t,\gamma_{k}})\Lambda_{0,\gamma_{k}|t}(\pi_{\gamma_{k}}^{\circ}-C_{t,\gamma_{k}})^{\prime}\Sigma_{\gamma_{k}}^{-1}\Big{)}$
(62)

the joint conditional density function is

$\displaystyle
f(\pi_{\gamma},\Sigma_{\gamma}|\bar{s}_{t},\mathcal{F}_{t})=c_{4}\prod_{k=1}^{r_{\gamma}}|\Sigma_{\gamma_{k}}|^{-(\nu_{0,\gamma_{k}}+q_{t,\gamma_{k}}+d+n+1)/2}$
$\displaystyle\times\exp\bigg{\\{}-\frac{1}{2}\sum_{k=1}^{r_{\gamma}}\mathrm{tr}\Big{(}\big{(}\bar{B}_{t,\gamma_{k}}+V_{0,\gamma_{k}}+(\pi_{\gamma_{k}}^{\circ}-C_{t,\gamma_{k}})\Lambda_{0,\gamma_{k}|t}(\pi_{\gamma_{k}}^{\circ}-C_{t,\gamma_{k}})^{\prime}\big{)}\Sigma_{\alpha_{k}}^{-1}\Big{)}\bigg{\\}},$
(63)

where the normalizing coefficient equals

$c_{4}:=\frac{1}{(2\pi)^{ndr_{\gamma}/2}\prod_{k=1}^{r_{\gamma}}|\Lambda_{0,\gamma_{k}|t}|^{-n/2}}\prod_{k=1}^{r_{\gamma}}\frac{\big{|}\bar{B}_{t,\gamma_{k}}+V_{0,\gamma_{k}}\big{|}^{(\nu_{0,\gamma_{k}}+q_{t,\gamma_{k}})/2}}{\Gamma_{n}\big{(}(\nu_{0,\gamma_{k}}+q_{t,\gamma_{k}})/2\big{)}2^{n(\nu_{0,\gamma_{k}}+q_{t,\gamma_{k}})/2}}.$
(64)

If we integrate the above equation by $\Sigma_{\gamma}$, then one obtains that

$\displaystyle f(\pi_{\gamma}^{\circ}|\bar{s}_{t},\mathcal{F}_{t})$
$\displaystyle=$
$\displaystyle\prod_{k=1}^{r_{\gamma}}\frac{|\Lambda_{0,\gamma_{k}|t}|^{n/2}\big{|}\bar{B}_{t,\gamma_{k}}+V_{0,\gamma_{k}}\big{|}^{-d/2}\Gamma_{n}\big{(}(\nu_{0,\gamma_{k}}+q_{t,\gamma_{k}}+d)/2\big{)}}{\pi^{nd/2}\Gamma_{n}((\nu_{0,\gamma_{k}}+q_{t,\gamma_{k}})/2)}$
$\displaystyle\times$
$\displaystyle\big{|}I_{n}+(\bar{B}_{t,\gamma_{k}}+V_{0,\gamma_{k}})^{-1}(\pi_{\gamma_{k}}^{\circ}-C_{t,\gamma_{k}})\Lambda_{0,\gamma_{k}|t}(\pi_{\gamma_{k}}^{\circ}-C_{t,\gamma_{k}})^{\prime}\big{|}^{-(\nu_{0,\gamma_{k}}+q_{t,\gamma_{k}}+d)/2}.$

Similarly, if $r_{\delta}>0$, it can be shown that

$\displaystyle f(\pi_{\delta}^{\circ}|\delta,\mathcal{F}_{0})$
$\displaystyle=$
$\displaystyle\prod_{\ell=1}^{r_{\delta}}\frac{|\Lambda_{0,\delta_{\ell}}|^{-n/2}|V_{0,\delta_{\ell}}|^{-d/2}\Gamma_{n}\big{(}(\nu_{0,\delta_{\ell}}+d)/2\big{)}}{\pi^{nd/2}\Gamma_{n}(\nu_{0,\delta_{\ell}}/2)}$
(66) $\displaystyle\times$
$\displaystyle\big{|}I_{n}+V_{0,\delta_{\ell}}^{-1}(\pi_{\delta_{\ell}}^{\circ}-\pi_{0,\delta_{\ell}}^{\circ})\Lambda_{0,\delta_{\ell}}^{-1}(\pi_{\delta_{\ell}}^{\circ}-\pi_{0,\delta_{\ell}}^{\circ})^{\prime}\big{|}^{-(\nu_{0,\delta_{\ell}}+d)/2}.$

Therefore, equation (2) holds. By the completing square method, the matrix
$\bar{B}_{t,\alpha_{k}}$ can be written by

$\displaystyle\bar{B}_{t,\alpha_{k}}$ $\displaystyle=$
$\displaystyle\big{(}\tilde{y}_{t,\alpha_{k}}^{\circ}-\pi_{0,\alpha_{k}}^{\circ}\Lambda_{0,\alpha_{k}}^{-1}\Lambda_{0,\alpha_{k}|t}^{-1}\mathsf{Y}_{t,\alpha_{k}}^{\circ}\Phi_{t,\alpha_{k}}^{-1}\big{)}\Phi_{t,\alpha_{k}}^{-1}\big{(}\tilde{y}_{t,\alpha_{k}}^{\circ}-\pi_{0,\alpha_{k}}^{\circ}\Lambda_{0,\alpha_{k}}^{-1}\Lambda_{0,\alpha_{k}|t}^{-1}\mathsf{Y}_{t,\alpha_{k}}^{\circ}\Phi_{t,\alpha_{k}}^{-1}\big{)}^{\prime}$
(67) $\displaystyle-$
$\displaystyle\pi_{0,\alpha_{k}}^{\circ}\Lambda_{0,\alpha_{k}}^{-1}\Lambda_{0,\alpha_{k}|t}^{-1}\mathsf{Y}_{t,\alpha_{k}}^{\circ}\Phi_{t,\alpha_{k}}^{-1}(\mathsf{Y}_{t,\alpha_{k}}^{\circ})^{\prime}\Lambda_{0,\alpha_{k}|t}^{-1}\Lambda_{0,\alpha_{k}}^{-1}(\pi_{0,\alpha_{k}}^{\circ})^{\prime}$
$\displaystyle+$
$\displaystyle\pi_{0,\alpha_{k}}^{\circ}\Lambda_{0,\alpha_{k}}^{-1}(\pi_{0,\alpha_{k}}^{\circ})^{\prime}-\pi_{0,\alpha_{k}}^{\circ}\Lambda_{0,\alpha_{k}}^{-1}\Lambda_{0,\alpha_{k}|t}^{-1}\Lambda_{0,\alpha_{k}}^{-1}(\pi_{0,\alpha_{k}}^{\circ})^{\prime},$

where
$\Phi_{t,\alpha_{k}}:=I_{q_{t,\alpha_{k}}}-(\mathsf{Y}_{t,\alpha_{k}}^{\circ})^{\prime}\Lambda_{0,\alpha_{k}|t}^{-1}\mathsf{Y}_{t,\alpha_{k}}^{\circ}$
is a symmetric $(q_{t,\alpha_{k}}\times q_{t,\alpha_{k}})$ matrix. We consider
the following product

$I_{t,\alpha_{k}}:=\big{(}I_{q_{t,\alpha_{k}}}+(\mathsf{Y}_{t,\alpha_{k}}^{\circ})^{\prime}\Lambda_{0,\alpha_{k}}\mathsf{Y}_{t,\alpha_{k}}^{\circ}\big{)}\big{(}I_{q_{t,\alpha_{k}}}-(\mathsf{Y}_{t,\alpha_{k}}^{\circ})^{\prime}\Lambda_{0,\alpha_{k}|t}^{-1}\mathsf{Y}_{t,\alpha_{k}}^{\circ}\big{)}.$
(68)

It equals

$\displaystyle I_{t,\alpha_{k}}$ $\displaystyle=$ $\displaystyle
I_{q_{t,\alpha_{k}}}+(\mathsf{Y}_{t,\alpha_{k}}^{\circ})^{\prime}\Lambda_{0,\alpha_{k}}\mathsf{Y}_{t,\alpha_{k}}^{\circ}-(\mathsf{Y}_{t,\alpha_{k}}^{\circ})^{\prime}\Lambda_{0,\alpha_{k}|t}^{-1}\mathsf{Y}_{t,\alpha_{k}}^{\circ}$
(69) $\displaystyle-$
$\displaystyle(\mathsf{Y}_{t,\alpha_{k}}^{\circ})^{\prime}\Lambda_{0,\alpha_{k}}\mathsf{Y}_{t,\alpha_{k}}^{\circ}(\mathsf{Y}_{t,\alpha_{k}}^{\circ})^{\prime}\Lambda_{0,\alpha_{k}|t}^{-1}\mathsf{Y}_{t,\alpha_{k}}^{\circ}.$

If we add and subtract the matrix $\Lambda_{0,\alpha_{k}}^{-1}$ into the term
$\mathsf{Y}_{t,\alpha_{k}}^{\circ}(\mathsf{Y}_{t,\alpha_{k}}^{\circ})^{\prime}$
in the last line of the above equation, then the product matrix equals
$I_{t,\alpha_{k}}=I_{q_{t,\alpha_{k}}}$. Consequently, the matrix
$I_{q_{t,\alpha_{k}}}+(\mathsf{Y}_{t,\alpha_{k}}^{\circ})^{\prime}\Lambda_{0,\alpha_{k}}\mathsf{Y}_{t,\alpha_{k}}^{\circ}$
is an inverse matrix of the matrix $\Phi_{t,\alpha_{k}}$, that is,

$\Phi_{t,\alpha_{k}}^{-1}=I_{q_{t,\alpha_{k}}}+(\mathsf{Y}_{t,\alpha_{k}}^{\circ})^{\prime}\Lambda_{0,\alpha_{k}}\mathsf{Y}_{t,\alpha_{k}}^{\circ}.$
(70)

Since it is a positive definite matrix, the matrix $\bar{B}_{t,\alpha_{k}}$ is
a positive semi–definite matrix. Now, we consider the term
$\Lambda_{0,\alpha_{k}}^{-1}\Lambda_{0,\alpha_{k}|t}^{-1}\mathsf{Y}_{t,\alpha_{k}}^{\circ}\Phi_{t,\alpha_{k}}^{-1}$
in the first line in equation (67). Similarly as before, by adding and
subtracting $\Lambda_{0,\alpha_{k}}^{-1}$ into the term
$\mathsf{Y}_{t,\alpha_{k}}^{\circ}(\mathsf{Y}_{t,\alpha_{k}}^{\circ})^{\prime}$,
one obtains that

$\Lambda_{0,\alpha_{k}}^{-1}\Lambda_{0,\alpha_{k}|t}^{-1}\mathsf{Y}_{t,\alpha_{k}}^{\circ}\Phi_{t,\alpha_{k}}^{-1}=\mathsf{Y}_{t,\alpha_{k}}^{\circ}.$
(71)

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

$|\Phi_{t,\alpha_{k}}^{-1}|=\big{|}I_{q_{t,\alpha_{k}}}+\big{(}\Lambda_{0,s_{t}}^{1/2}\mathsf{Y}_{t,\alpha_{k}}^{\circ}\big{)}^{\prime}\Lambda_{0,s_{t}}^{1/2}\mathsf{Y}_{t,\alpha_{k}}^{\circ}\big{|}=\big{|}I_{d}+\Lambda_{0,s_{t}}^{1/2}\mathsf{Y}_{t,\alpha_{k}}^{\circ}\big{(}\mathsf{Y}_{t,\alpha_{k}}^{\circ}\big{)}^{\prime}\big{(}\Lambda_{0,s_{t}}^{1/2}\big{)}^{\prime}\big{|}.$
(72)

That completes the proof of the Proposition. ∎

It follows from equations (2) and (47) that sub coefficient vectors and sub
covariance matrices are conditional independent. Note that the conditional
independence is consistent with the assumption (10). From equations (2) and
(47) one deduces that for $k=1,\dots,r_{\alpha}$, the conditional density
functions of the coefficient vector $\pi_{\alpha_{k}}$ and the covariance
matrix $\Sigma_{\alpha_{k}}$ are given by

$\displaystyle
f(\pi_{\alpha_{k}}|\Sigma_{\alpha_{k}},\alpha_{k},\mathcal{F}_{t})=\frac{1}{(2\pi)^{nd/2}|A_{\alpha_{k}}|^{1/2}}\exp\bigg{\\{}-\frac{1}{2}\Big{(}\pi_{\alpha_{k}}-A_{\alpha_{k}}b_{\alpha_{k}}\Big{)}^{\prime}A_{\alpha_{k}}^{-1}\Big{(}\pi_{\alpha_{k}}-A_{\alpha_{k}}b_{\alpha_{k}}\Big{)}\bigg{\\}}$
(73)

and

$\displaystyle f(\Sigma_{\alpha_{k}}|\alpha_{k},\mathcal{F}_{t})$
$\displaystyle=$
$\displaystyle\frac{\big{|}\bar{B}_{t,\alpha_{k}}+V_{0,\alpha_{k}}\big{|}^{(\nu_{0,\alpha_{k}}+q_{t,\alpha_{k}})/2}}{\Gamma_{n}\big{(}(\nu_{0,\alpha_{k}}+q_{t,\alpha_{k}})/2\big{)}2^{n(\nu_{0,\alpha_{k}}+q_{t,\alpha_{k}})/2}}|\Sigma_{\alpha_{k}}|^{-(\nu_{0,\alpha_{k}}+q_{t,\alpha_{k}}+n+1)/2}$
(74) $\displaystyle\times$
$\displaystyle\exp\bigg{\\{}-\frac{1}{2}\mathrm{tr}\Big{(}\big{(}\bar{B}_{t,\alpha_{k}}+V_{0,\alpha_{k}}\big{)}\Sigma_{\alpha_{k}}^{-1}\Big{)}\bigg{\\}},$

respectively. Thus, the conditional distribution functions of the coefficient
vector $\pi_{\alpha_{k}}$ and the covariance matrix $\Sigma_{\alpha_{k}}$ are
multivariate normal and inverse Wishart, respectively. Also, it follows from
equation (2) that the conditional density function of the random coefficient
matrix $\pi_{\beta}^{\circ}$ equals products of matrix variate student $t$
density functions. The conditional density function of the random matrix
$\pi_{\beta}^{\circ}$ can be used to impulse response analysis. Because
marginal density functions of the random coefficient matrix
$\pi_{\beta}^{\circ}$ are the matrix variate student $t$, their means are
given by

$\mathbb{E}\big{[}\pi_{\gamma_{k}}^{\circ}\big{|}\gamma_{k},\mathcal{F}_{t}\big{]}=C_{t,\gamma_{k}},~{}~{}~{}k=1,\dots,r_{\gamma}$
(75)

and

$\mathbb{E}\big{[}\pi_{\delta_{k}}^{\circ}\big{|}\delta_{k},\mathcal{F}_{0}\big{]}=\pi_{0,\delta_{k}},~{}~{}~{}k=1,\dots,r_{\delta}.$
(76)

Because, according to equation (2), density function (44) has a form of the
matrix variate student $t$ distribution, we refer to the density function as a
conditional matrix variate student $t$ density function. Furthermore, it
follows from equation (9) and (44) that conditional on the regime vector $s$
and information $\mathcal{F}_{t}$, a density function of future values of the
vector of endogenous variables is given by the following equation

$\displaystyle f(\bar{y}_{t}^{c}|s,\mathcal{F}_{t})$ $\displaystyle=$
$\displaystyle\frac{1}{\pi^{n(T-t)/2}}\prod_{k=1}^{r_{\alpha}}\frac{|\Lambda_{0,\alpha_{k}|t}|^{n/2}\Gamma_{n}\big{(}(\nu_{0,\alpha_{k}}+q_{T,\alpha_{k}})/2\big{)}|\bar{B}_{t,\alpha_{k}}+V_{0,\alpha_{k}}|^{(\nu_{0,\alpha_{k}}+q_{t,\alpha_{k}})/2}}{|\Lambda_{0,\alpha_{k}|T}|^{n/2}\Gamma_{n}((\nu_{0,\alpha_{k}}+q_{t,\alpha_{k}}/2)\big{|}\bar{B}_{T,\alpha_{k}}+V_{0,\alpha_{k}}\big{|}^{(\nu_{0,\alpha_{k}}+q_{T,\alpha_{k}})/2}}$
(77) $\displaystyle\times$
$\displaystyle\prod_{\ell=1}^{r_{\delta}}\frac{|\Lambda_{0,\delta_{\ell}}^{-1}|^{n/2}\Gamma_{n}\big{(}(\nu_{0,\delta_{\ell}}+q_{T,\delta_{\ell}})/2\big{)}|V_{0,\delta_{\ell}}|^{\nu_{0,\delta_{\ell}}/2}}{|\Lambda_{0,\delta_{\ell}|T}|^{n/2}\Gamma_{n}(\nu_{0,\delta_{\ell}}/2)\big{|}\bar{B}_{T,\delta_{\ell}}+V_{0,\delta_{\ell}}\big{|}^{(\nu_{0,\delta_{\ell}}+q_{T,\delta_{\ell}})/2}}.$

Let us consider the following matrix

$\Theta_{t,\alpha_{k}}:=\Lambda_{0,\alpha_{k}}-\Lambda_{0,\alpha_{k}}\mathsf{Y}_{t,\alpha_{k}}^{\circ}\Phi_{t,\alpha_{k}}(\mathsf{Y}_{t,\alpha_{k}}^{\circ})^{\prime}\Lambda_{0,\alpha_{k}}.$
(78)

Since
$\Phi_{t,\alpha_{k}}=I_{q_{t,\alpha_{k}}}-(\mathsf{Y}_{t,\alpha_{k}}^{\circ})^{\prime}\Lambda_{0,\alpha_{k}|t}^{-1}\mathsf{Y}_{t,\alpha_{k}}^{\circ}$,
we have that

$\Theta_{t,\alpha_{k}}=\Lambda_{0,\alpha_{k}}-\Lambda_{0,\alpha_{k}}\mathsf{Y}_{t,\alpha_{k}}^{\circ}(\mathsf{Y}_{t,\alpha_{k}}^{\circ})^{\prime}\Lambda_{0,\alpha_{k}}+\Lambda_{0,\alpha_{k}}\mathsf{Y}_{t,\alpha_{k}}^{\circ}(\mathsf{Y}_{t,\alpha_{k}}^{\circ})^{\prime}\Lambda_{0,\alpha_{k}|t}^{-1}\mathsf{Y}_{t,\alpha_{k}}^{\circ}(\mathsf{Y}_{t,\alpha_{k}}^{\circ})^{\prime}\Lambda_{0,\alpha_{k}}.$
(79)

If we substitute the matrix
$\mathsf{Y}_{t,\alpha_{k}}^{\circ}(\mathsf{Y}_{t,\alpha_{k}}^{\circ})^{\prime}=\Lambda_{0,\alpha_{k}|t}-\Lambda_{0,\alpha_{k}}^{-1}$
into the above equation, one obtains
$\Theta_{t,\alpha_{k}}=\Lambda_{0,\alpha_{k}|t}^{-1}$. Consequently, by the
formula of partitioned matrix’s inverse (e.g., see Lutkepohl05) it can be
shown that the matrix $\bar{B}_{T,\alpha_{k}}+V_{0,\alpha_{k}}$ equals

$\displaystyle\bar{B}_{T,\alpha_{k}}+V_{0,\alpha_{k}}$ $\displaystyle=$
$\displaystyle\bar{B}_{t,\alpha_{k}}+V_{0,\alpha_{k}}$ (80) $\displaystyle+$
$\displaystyle\big{(}y_{t,\alpha_{k}}^{*}-\pi_{0,\alpha_{k}}^{\circ}\mathsf{Y}_{t,\alpha_{k}}^{*}-(y_{t,\alpha_{k}}^{\circ}-\pi_{0,\alpha_{k}}^{\circ}\mathsf{Y}_{t,\alpha_{k}}^{\circ})\Phi_{t,\alpha_{k}}(\mathsf{Y}_{t,\alpha_{k}}^{\circ})^{\prime}\Lambda_{0,\alpha_{k}}\mathsf{Y}_{t,\alpha_{k}}^{*}\big{)}$
$\displaystyle\times$
$\displaystyle(I_{q_{t,\alpha}^{*}}+\mathsf{Y}_{t,\alpha_{k}}^{*}\Lambda_{0,\alpha_{k}|t}^{-1}\mathsf{Y}_{t,\alpha_{k}}^{*})^{-1}$
$\displaystyle\times$
$\displaystyle\big{(}y_{t,\alpha_{k}}^{*}-\pi_{0,\alpha_{k}}^{\circ}\mathsf{Y}_{t,\alpha_{k}}^{*}-(y_{t,\alpha_{k}}^{\circ}-\pi_{0,\alpha_{k}}^{\circ}\mathsf{Y}_{t,\alpha_{k}}^{\circ})\Phi_{t,\alpha_{k}}(\mathsf{Y}_{t,\alpha_{k}}^{\circ})^{\prime}\Lambda_{0,\alpha_{k}}\mathsf{Y}_{t,\alpha_{k}}^{*}\big{)}^{\prime},$

where the matrices $y_{t,\alpha_{k}}^{*}$ and $\mathsf{Y}_{t,\alpha_{k}}^{*}$
are come from the matrices
$y_{T,\alpha_{k}}^{\circ}=[y_{t,\alpha_{k}}^{\circ}:y_{t,\alpha_{k}}^{*}]$ and
$\mathsf{Y}_{T,\alpha_{k}}^{\circ}=[\mathsf{Y}_{t,\alpha_{k}}^{\circ}:\mathsf{Y}_{t,\alpha_{k}}^{*}]$,
corresponding to the random vectors $\bar{y}_{t}$ and $\bar{y}_{t}^{c}$, and
$q_{t,\alpha_{k}}^{*}=q_{T,\alpha_{k}}-q_{t,\alpha_{k}}$. As a result, density
function (77) is represented by a product of the conditional matrix variate
student $t$ density functions. Using the idea of proof of Proposition 2, one
obtains conditional density function of the random vector $y_{t}$ for given
$s_{t}$ and $\mathcal{F}_{t}$

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

where

$\displaystyle
B_{t,s_{t}}:=\frac{1}{1+\mathsf{Y}_{t}^{\prime}\Lambda_{0,s_{t}}\mathsf{Y}_{t}}\big{(}y_{t}-\pi_{0,s_{t}}^{\circ}\mathsf{Y}_{t}\big{)}\big{(}y_{t}-\pi_{0,s_{t}}^{\circ}\mathsf{Y}_{t}\big{)}^{\prime}$
(82)

is a symmetric positive semi–definite $(n\times n)$ matrix. Note that the
conditional density function only depends on the regime $s_{t}$ and does not
depend on the other regimes. To calculate smoothed probabilities, the
conditional density function will be used, see below.

Let us assume that the prior density functions of each row of the transition
probability matrix $\mathsf{P}$ follow Dirichlet distribution and they are
mutually independent. Under the assumption, a joint density function of them
is given by

$f(\mathsf{P}|\mathcal{F}_{0})=\prod_{i=0}^{N}\frac{\Gamma\big{(}\sum_{j=1}^{N}\alpha_{ij}\big{)}}{\prod_{j=1}^{N}\Gamma(\alpha_{ij})}\prod_{j=1}^{N}p_{ij}^{\alpha_{ij}-1}$
(83)

where $\Gamma(\cdot)$ is the gamma function and the parameters of Dirichlet
distribution satisfy $\alpha_{ij}>0$ for $i=0,\dots,N$ and $j=1,\dots,N$. Let
us denote $i$–th row of the random transition probability matrix $\mathsf{P}$
by $\mathsf{P}_{i}$, corresponding prior hyperparameter by
$\alpha_{i}:=(\alpha_{i1},\dots,\alpha_{iN})^{\prime}$, and Dirichlet
distribution by $\text{Dir}(\alpha_{i})$. Then, the following Lemma holds.

######  Proposition 3.

Let for $i=0,\dots,N$, $\mathsf{P}_{i}\sim\mathrm{Dir}(\alpha_{i})$ and they
are mutually independent. Then, the followings are hold

  * (i)

for $t=1,\dots,T$, conditional on the information $\mathcal{F}_{0}$, a density
function of regime vector $\bar{s}_{t}$ is given by

$f(\bar{s}_{t}|\mathcal{F}_{0})=\prod_{i=0}^{N}\frac{\Gamma\big{(}\sum_{j=1}^{N}\alpha_{ij}\big{)}}{\prod_{j=1}^{N}\Gamma(\alpha_{ij})}\frac{\prod_{j=1}^{N}\Gamma(\alpha_{ij}+n_{ij}(\bar{s}_{t}))}{\Gamma\big{(}\sum_{j=1}^{N}(\alpha_{ij}+n_{ij}(\bar{s}_{t}))\big{)}}$
(84)

  * (ii)

and for $t=2,\dots,T$, conditional on the regime vector $\bar{s}_{t}$ and
information $\mathcal{F}_{0}$, a density function of regime $s_{t}$ is given
by

$f(s_{t}|\bar{s}_{t-1},\mathcal{F}_{0})=\frac{\alpha_{s_{t-1}s_{t}}+n_{s_{t-1}s_{t}}(\bar{s}_{t-1})}{\sum_{s_{t}=1}^{N}\big{(}\alpha_{s_{t-1}s_{t}}+n_{s_{t-1}s_{t}}(\bar{s}_{t-1})\big{)}},$
(85)

where the random variable $n_{ij}(\bar{s}_{t})$ equals

$n_{ij}(\bar{s}_{t}):=\\#\big{\\{}m\in\\{0,1,\dots,t-1\\}\big{|}s_{m-1}=i,s_{m}=j,~{}m=2,\dots,t\big{\\}}$
(86)

for $t=2,\dots,T$, $i=1,\dots,N$, and $j=1,\dots,N$ and

$n_{ij}(s_{1}):=\begin{cases}1&\mathrm{if}~{}~{}~{}i=0,~{}s_{1}=j\\\
0&\mathrm{if}~{}~{}~{}\mathrm{otherwise}\end{cases}$ (87)

for $i=0,\dots,N$ and $j=1,\dots,N$.

###### Proof.

Since for $t=2,\dots,T$, the random variable $n_{ij}(\bar{s}_{t})$ represents
a number of consequential elements, which equals $(i,j)$ of the regime vector
$\bar{s}_{t}$, we have that

$f(\bar{s}_{t}|\mathsf{P},\mathcal{F}_{0})=\prod_{m=1}^{t}p_{s_{m-1}s_{m}}=\prod_{i=0}^{N}\prod_{j=1}^{N}p_{ij}^{n_{ij}(\bar{s}_{t})}.$
(88)

Consequently, it follows from the joint density function of the random
transition probability matrix $\mathsf{P}$, which is given by equation (75)
that

$f(\bar{s}_{t}|\mathcal{F}_{0})=\int_{\mathsf{P}}f(\bar{s}_{t}|\mathsf{P},\mathcal{F}_{0})f(\mathsf{P}|\mathcal{F}_{0})d\mathsf{P}=\prod_{i=0}^{N}\frac{\Gamma\big{(}\sum_{j=1}^{N}\alpha_{ij}\big{)}}{\prod_{j=1}^{N}\Gamma(\alpha_{ij})}\frac{\prod_{j=1}^{N}\Gamma(\alpha_{ij}+n_{ij}(\bar{s}_{t}))}{\Gamma\big{(}\sum_{j=1}^{N}(\alpha_{ij}+n_{ij}(\bar{s}_{t}))\big{)}}.$
(89)

On the other hand, as
$f(s_{t}|\bar{s}_{t},\mathcal{F}_{0})=f(\bar{s}_{t}|\mathcal{F}_{0})/f(\bar{s}_{t-1}|\mathcal{F}_{0})$,
one obtains

$f(s_{t}|\bar{s}_{t-1},\mathcal{F}_{0})=\prod_{i=0}^{N}\frac{\prod_{j=1}^{N}\Gamma(\alpha_{ij}+n_{ij}(\bar{s}_{t}))}{\Gamma\big{(}\sum_{j=1}^{N}(\alpha_{ij}+n_{ij}(\bar{s}_{t}))\big{)}}\frac{\Gamma\big{(}\sum_{j=1}^{N}(\alpha_{ij}+n_{ij}(\bar{s}_{t-1}))\big{)}}{\prod_{j=1}^{N}\Gamma(\alpha_{ij}+n_{ij}(\bar{s}_{t-1}))}.$
(90)

Consequently, since
$n_{ij}(\bar{s}_{t})=n_{ij}(\bar{s}_{t-1})+\delta_{ij}(s_{t})$ for
$t=2,\dots,T$ and $\Gamma(x+1)=x\Gamma(x)$ for $x>0$, we find that

$f(s_{t}|\bar{s}_{t-1},\mathcal{F}_{0})=\frac{\alpha_{s_{t-1}s_{t}}+n_{s_{t-1}s_{t}}(\bar{s}_{t-1})}{\sum_{s_{t}=1}^{N}\big{(}\alpha_{s_{t-1}s_{t}}+n_{s_{t-1}s_{t}}(\bar{s}_{t-1})\big{)}},$
(91)

where the random variable $\delta_{ij}(s_{t})$ equals

$\delta_{ij}(s_{t})=\begin{cases}1&\mathrm{if}~{}~{}~{}s_{t-1}=i,~{}s_{t}=j,\\\
0&\mathrm{if}~{}~{}~{}\mathrm{otherwise}.\end{cases}$ (92)

That completes the proof. ∎

It is worth mentioning that according to equation (85) in the above
Proposition, conditional on $\mathcal{F}_{0}$, the regime–switching process
$s_{t}$ is not a Markov chain.

####  3.1.2 Characteristic Function

By equation (62), for $k=1,\dots,r_{\gamma}$, equation (73) can be written by

$\displaystyle
f(\pi_{\gamma_{k}}^{\circ}|\Sigma_{\gamma_{k}},\gamma_{k},\mathcal{F}_{t})$
$\displaystyle=$
$\displaystyle\frac{1}{(2\pi)^{nd/2}|\Lambda_{0,\gamma_{k}|t}|^{-n/2}|\Sigma_{\gamma_{k}}|^{d/2}}$
(93) $\displaystyle\times$
$\displaystyle\exp\bigg{\\{}-\frac{1}{2}\mathrm{tr}\Big{(}\pi_{\gamma_{k}}^{\circ}-C_{t,\gamma_{k}})\Lambda_{0,\gamma_{k}|t}(\pi_{\gamma_{k}}^{\circ}-C_{t,\gamma_{k}})^{\prime}\Sigma_{\gamma_{k}}^{-1}\Big{)}\bigg{\\}}.$

To obtain characteristic function of the random coefficient matrix
$\pi_{\gamma_{k}}^{\circ}$ for given regime $\gamma_{k}$ and information
$\mathcal{F}_{t}$, we use the matrix generalized inverse Gaussian (MGIG)
distribution. For a positive definite $(n\times n)$ matrix $\Sigma$, the
density function of the MGIG distribution is given by

$f(\Sigma)=\frac{2^{n\lambda}}{|\mathsf{A}|^{\lambda}\mathcal{B}_{\lambda}\big{(}\frac{1}{4}\mathsf{B}\mathsf{A}\big{)}}|\Sigma|^{\lambda-(n+1)/2}\exp\bigg{\\{}-\frac{1}{2}\mathrm{tr}\big{(}\mathsf{A}\Sigma^{-1}+\mathsf{B}\Sigma\big{)}\bigg{\\}},$
(94)

where $\mathcal{B}_{\lambda}(\cdot)$ is the matrix argument modified Bessel
function of the second kind with index $\lambda$, which is defined by

$\mathcal{B}_{\lambda}\bigg{(}\frac{1}{4}\mathsf{B}\mathsf{A}\bigg{)}:=\bigg{|}\frac{1}{2}\mathsf{B}\bigg{|}^{-\lambda}\int_{\Sigma>0}|\Sigma|^{-\lambda-(n+1)/2}\exp\bigg{\\{}-\frac{1}{2}\mathrm{tr}\big{(}\mathsf{A}\Sigma^{-1}+\mathsf{B}\Sigma\big{)}\bigg{\\}}d\Sigma$
(95)

and for $n\geq 2$, the index $\lambda\in\mathbb{R}$ and the $(n\times n)$
matrices $\mathsf{A}$ and $\mathsf{B}$ satisfy

$\begin{cases}\mathsf{A}\geq
0,~{}\mathsf{B}>0&\text{if}~{}~{}~{}\lambda\geq\frac{1}{2}\\\
\mathsf{A}>0,~{}\mathsf{B}>0&\text{if}~{}~{}~{}-\frac{1}{2}(n-1)\leq\lambda<\frac{1}{2}\\\
\mathsf{A}>0,~{}\mathsf{B}\geq
0&\text{if}~{}~{}~{}\lambda<-\frac{1}{2}(n-1)\\\ \end{cases},$ (96)

see Butler98. Its one dimensional version is called generalized inverse
Gaussian distribution and it is widely used to model returns of financial
assets, see McNeil05. It is the well–known fact that a characteristic function
of the random coefficient matrix $\pi_{\gamma_{k}}^{\circ}$ for given regime
$\gamma_{k}$, covariance matrix $\Sigma_{\gamma_{k}}$, and information
$\mathcal{F}_{t}$ is given by

$\displaystyle\varphi(Z_{\gamma_{k}}|\Sigma_{\gamma_{k}},\gamma_{k},\mathcal{F}_{t})$
$\displaystyle=$
$\displaystyle\mathbb{E}\big{[}\exp\big{\\{}\mathrm{tr}\big{(}\mathrm{i}Z_{\gamma_{k}}^{\prime}\pi_{\gamma_{k}}^{\circ}\big{)}\big{\\}}\big{|}\Sigma_{\gamma_{k}},\gamma_{k},\mathcal{F}_{t}\big{]}$
(97) $\displaystyle=$
$\displaystyle\exp\bigg{\\{}\mathrm{tr}\Big{(}\mathrm{i}Z_{\gamma_{k}}^{\prime}C_{t,\gamma_{k}}-\frac{1}{2}Z_{\gamma_{k}}^{\prime}\Lambda_{0,\gamma_{k}|t}^{-1}Z_{\gamma_{k}}\Sigma_{\gamma_{k}}\Big{)}\bigg{\\}},$

where $\mathrm{i}=\sqrt{-1}$ is the imaginary unit, $Z_{\gamma_{k}}$ is an
$(n\times d)$ matrix, corresponding to the regime $\gamma_{k}$. Consequently,
by the iterated expectation formula, conditional density function (74), and
the above characteristic function of the matrix normal distribution,
conditional on the regime $\gamma_{k}$ and information $\mathcal{F}_{t}$, a
characteristic function of the random coefficient matrix
$\pi_{\gamma_{k}}^{\circ}$ is obtained by

$\displaystyle\varphi(Z_{\gamma_{k}}|\gamma_{k},\mathcal{F}_{t})$
$\displaystyle=$
$\displaystyle\mathbb{E}\big{[}\exp\big{\\{}\mathrm{tr}\big{(}\mathrm{i}Z_{\gamma_{k}}^{\prime}\pi_{\gamma_{k}}^{\circ}\big{)}\big{\\}}\big{|}\gamma_{k},\mathcal{F}_{t}\big{]}$
(98) $\displaystyle=$
$\displaystyle\frac{\exp\big{\\{}\mathrm{tr}\big{(}\mathrm{i}Z_{\gamma_{k}}^{\prime}C_{t,\gamma_{k}}\big{)}\big{\\}}\big{|}\bar{B}_{t,\gamma_{k}}+V_{0,\gamma_{k}}\big{|}^{\nu_{0,\gamma_{k}}+q_{t,\gamma_{k}}}|Z_{\gamma_{k}}^{\prime}\Lambda_{0,\gamma_{k}|t}^{-1}Z_{\gamma_{k}}|^{\nu_{0,\gamma_{k}}+q_{t,\gamma_{k}}}}{\Gamma_{n}\big{(}(\nu_{0,\gamma_{k}}+q_{t,\gamma_{k}})/2\big{)}2^{n(\nu_{0,\gamma_{k}}+q_{t,\gamma_{k}})}}$
$\displaystyle\times$
$\displaystyle\mathcal{B}_{\nu_{0,\gamma_{k}}+q_{t,\gamma_{k}}}\big{(}Z_{\gamma_{k}}^{\prime}\Lambda_{0,\gamma_{k}|t}^{-1}Z_{\gamma_{k}}\big{(}\bar{B}_{t,\gamma_{k}}+V_{0,\gamma_{k}}\big{)}/4\big{)}.$

Similarly, it can be shown that

$\displaystyle\varphi(Z_{\delta_{k}}|\delta_{k},\mathcal{F}_{0})$
$\displaystyle=$
$\displaystyle\mathbb{E}\big{[}\exp\big{\\{}\mathrm{tr}\big{(}\mathrm{i}Z_{\delta_{k}}^{\prime}\pi_{\delta_{k}}^{\circ}\big{)}\big{\\}}\big{|}\delta_{k},\mathcal{F}_{0}\big{]}$
(99) $\displaystyle=$
$\displaystyle\frac{\exp\big{\\{}\mathrm{tr}\big{(}\mathrm{i}Z_{\delta_{k}}^{\prime}\pi_{0,\delta_{k}}^{\circ}\big{)}\big{\\}}\big{|}V_{0,\delta_{k}}\big{|}^{\nu_{0,\delta_{k}}}|Z_{\delta_{k}}^{\prime}\Lambda_{0,\delta_{k}}Z_{\delta_{k}}|^{\nu_{0,\delta_{k}}}}{\Gamma_{n}\big{(}(\nu_{0,\delta_{k}})/2\big{)}2^{n\nu_{0,\delta_{k}}}}$
$\displaystyle\times$
$\displaystyle\mathcal{B}_{\nu_{0,\delta_{k}}}\big{(}Z_{\delta_{k}}^{\prime}\Lambda_{0,\delta_{k}}Z_{\delta_{k}}V_{0,\delta_{k}}/4\big{)}.$

The above characteristic functions can be used to obtain raw moments of the
random coefficient matrix $\pi_{\beta}^{\circ}$ for given the regime vector
$s$ and information $\mathcal{F}_{t}$. For example, since conditional on $s$
and $\mathcal{F}_{t}$, for $k=1,\dots,r_{\gamma}$ and
$\ell=1,\dots,r_{\delta}$, coefficient matrices $\pi_{\beta_{k}}$ and
$\pi_{\delta_{\ell}}$ are independent, we have that

$\displaystyle\mathbb{E}\Big{[}\big{(}\pi_{\gamma_{1}}^{\circ}\big{)}_{i_{\gamma_{1}},j_{\gamma_{1}}}^{m_{\gamma_{1}}}\dots\big{(}\pi_{\gamma_{r_{\gamma}}}^{\circ}\big{)}_{i_{\gamma_{r_{\gamma}}},j_{\gamma_{r_{\gamma}}}}^{m_{\gamma_{r_{\gamma}}}}\big{(}\pi_{\delta_{1}}^{\circ}\big{)}_{i_{\delta_{1}},j_{\delta_{1}}}^{m_{\delta_{1}}}\dots\big{(}\pi_{\delta_{r_{\delta}}}^{\circ}\big{)}_{i_{\delta_{r_{\delta}}},j_{\delta_{r_{\delta}}}}^{m_{\delta_{r_{\delta}}}}\Big{|}s,\mathcal{F}_{t}\Big{]}$
$\displaystyle=\prod_{k=1}^{r_{\gamma}}\frac{1}{\mathrm{i}^{m_{\gamma_{k}}}}\frac{\partial^{m_{\gamma_{k}}}\varphi(Z_{\gamma_{k}}|\gamma_{k},\mathcal{F}_{t})}{\partial(Z_{\gamma_{k}})_{i_{\gamma_{k}},j_{\gamma_{k}}}^{m_{\gamma_{k}}}}\times\prod_{\ell=1}^{r_{\delta}}\frac{1}{\mathrm{i}^{m_{\delta_{\ell}}}}\frac{\partial^{m_{\delta_{\ell}}}\varphi(Z_{\delta_{\ell}}|\delta_{\ell},\mathcal{F}_{0})}{\partial(Z_{\gamma_{k}})_{i_{\delta_{\ell}},j_{\delta_{\ell}}}^{m_{\gamma_{k}}}},$
(100)

where for a generic $(n\times m)$ matrix $O$, $(O)_{i,j}$ denotes an
$(i,j)$–th element of the matrix $O$ for $i=1,\dots,n$ and $j=1,\dots,m$. The
partial derivatives can be calculated by the numerical methods. The raw
moments may be used to obtain forecast of the vector of endogenous variables.
In particular, conditional on $\bar{s}_{t+2}$, the optimal forecast, which
minimizes the mean squared errors for forecast horizon 2 at forecast origin
$t$ equals an expectation of the vector of endogenous variables at time
$(t+2)$ for given $\mathcal{F}_{t}$. Thus, the forecast is given by the
following equation

$\displaystyle\mathbb{E}\big{[}y_{t+2}\big{|}\bar{s}_{t+2},\mathcal{F}_{t}\big{]}$
$\displaystyle=$
$\displaystyle\mathbb{E}\big{[}A_{0,s_{t+2}}\big{|}\bar{s}_{t+2},\mathcal{F}_{t}\big{]}+\mathbb{E}\big{[}A_{1,s_{t+2}}A_{0,s_{t+1}}\big{|}\bar{s}_{t+2},\mathcal{F}_{t}\big{]}\psi_{t+1}$
$\displaystyle+$
$\displaystyle\sum_{k=1}^{p}\mathbb{E}\big{[}A_{1,s_{t+2}}A_{k,s_{t+1}}|\bar{s}_{t+2},\mathcal{F}_{t}\big{]}y_{t+1-k}+\sum_{k=2}^{p}\mathbb{E}\big{[}A_{k,s_{t+2}}|\bar{s}_{t+2},\mathcal{F}_{t}\big{]}y_{t+2-k},$

where the conditional expectations
$\mathbb{E}\big{[}A_{k,s_{t+2}}\big{|}\bar{s}_{t+2},\mathcal{F}_{t}\big{]}$
for $k=0,2,\dots,p$ are calculated by equations (75) and (76) and the
conditional expectations
$\mathbb{E}\big{[}A_{1,s_{t+2}}A_{k,s_{t+1}}\big{|}\bar{s}_{t+2},\mathcal{F}_{t}\big{]}$
for $k=0,\dots,p$ are calculated equation (3.1.2). To illustrative purpose, we
assume that
$s_{t+1},s_{t+2}\in\mathcal{A}_{\bar{s}_{t}}\cap\mathcal{A}_{\bar{s}_{t}^{c}}$
and positions of the regimes $s_{t+1}$ and $s_{t+2}$ in the regime vector
$\gamma$ are $k_{1*}$ and $k_{2*}$, respectively, that is,

$k_{i*}=\\{k\in\\{1,\dots,r_{\gamma}\\}|\gamma_{k}=s_{t+i},k=1,\dots,r_{\gamma}\\}$
(102)

for $i=1,2$. Then, we have that for $i=1,\dots,n$, $j=1,\dots,l$, and $k=0$,

$\Big{(}\mathbb{E}\big{[}A_{1,s_{t+2}}A_{k,s_{t+1}}|\bar{s}_{t+2},\mathcal{F}_{t}\big{]}\Big{)}_{i,j}=-\sum_{\ell=1}^{n}\frac{\partial\varphi(Z_{\gamma_{k_{2*}}}|\gamma_{k_{2*}},\mathcal{F}_{t})}{\partial(Z_{\gamma_{k_{2*}}})_{i,l+\ell}}\frac{\partial\varphi(Z_{\gamma_{k_{1*}}}|\gamma_{k_{1*}},\mathcal{F}_{t})}{\partial(Z_{\gamma_{k_{1*}}})_{\ell,j}},$
(103)

for $s_{t+1}=s_{t+2}$, $i=1,\dots,n$, and $k=1$,

$\displaystyle\Big{(}\mathbb{E}\big{[}A_{1,s_{t+2}}A_{k,s_{t+1}}|\bar{s}_{t+2},\mathcal{F}_{t}\big{]}\Big{)}_{i,i}$
$\displaystyle=-\sum_{\ell=1,\ell\neq
i}^{n}\frac{\partial\varphi(Z_{\gamma_{k_{1*}}}|\gamma_{k_{1*}},\mathcal{F}_{t})}{\partial(Z_{\gamma_{k_{1*}}})_{i,l+\ell}}\frac{\partial\varphi(Z_{\gamma_{k_{1*}}}|\gamma_{k_{1*}},\mathcal{F}_{t})}{\partial(Z_{\gamma_{k_{1*}}})_{l+\ell,i}}-\frac{\partial^{2}\varphi(Z_{\gamma_{k_{1*}}}|\gamma_{k_{1*}},\mathcal{F}_{t})}{\partial(Z_{\gamma_{k_{1*}}})_{i,l+i}^{2}},$
(104)

and for other cases ($i,j=1,\dots,n$ and $k=1,\dots,p$),

$\displaystyle\Big{(}\mathbb{E}\big{[}A_{1,s_{t+2}}A_{k,s_{t+1}}|\bar{s}_{t+2},\mathcal{F}_{t}\big{]}\Big{)}_{i,j}=-\sum_{\ell=1}^{n}\frac{\partial\varphi(Z_{\gamma_{k_{2*}}}|\gamma_{k_{2*}},\mathcal{F}_{t})}{\partial(Z_{\gamma_{k_{2*}}})_{i,l+\ell}}\frac{\partial\varphi(Z_{\gamma_{k_{1*}}}|\gamma_{k_{1*}},\mathcal{F}_{t})}{\partial(Z_{\gamma_{k_{1*}}})_{l+(k-1)p+\ell,j}}.$
(105)

Because the exact calculation of forecast of the process of endogenous
variables is complicated, we consider an approximation, which is used to
calculate the forecast of the endogenous variables in Banbura10. For
$u=t+1,\dots,T$, by the iterated expectation formula, conditional on the
regime vector $s$, the exact forecast is given by the following equation

$\mathbb{E}\big{[}y_{u}\big{|}s,\mathcal{F}_{t}\big{]}=\mathbb{E}\big{[}\Pi_{s_{u}}\mathsf{Y}_{u}\big{|}s,\mathcal{F}_{t}\big{]}=\mathbb{E}\big{[}\mathbb{E}\big{[}\Pi_{s_{u}}\big{|}s,\mathcal{F}_{u-1}\big{]}\mathsf{Y}_{u}\big{|}s,\mathcal{F}_{t}\big{]}.$
(106)

Banbura10 approximate the last expression by
$\mathbb{E}\big{[}\Pi_{s_{u}}\big{|}s,\mathcal{F}_{u-1}\big{]}\mathbb{E}\big{[}\mathsf{Y}_{u}\big{|}s,\mathcal{F}_{t}\big{]}.$
Consequently, the forecast is approximated by

$\mathbb{E}\big{[}y_{u}\big{|}s,\mathcal{F}_{t}\big{]}\approx\mathbb{E}\big{[}\Pi_{s_{u}}\big{|}s,\mathcal{F}_{u-1}\big{]}\mathbb{E}\big{[}\mathsf{Y}_{u}\big{|}s,\mathcal{F}_{t}\big{]}.$
(107)

For $t=u-1$, the approximation becomes exact, namely,

$\mathbb{E}\big{[}y_{u}\big{|}s,\mathcal{F}_{u-1}\big{]}=\mathbb{E}\big{[}\Pi_{s_{u}}\big{|}s,\mathcal{F}_{u-1}\big{]}\mathsf{Y}_{u}.$
(108)

However, for $u=t+2,\dots,T$, one should study the quality of the very simple
approximation.

####  3.1.3 Minnesota Prior

In practice, one usually adopts 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,s_{t}}$ and all elements
of other matrices $A_{0,s_{t}},A_{2,s_{t}},\dots,A_{p,s_{t}}$ 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}$. For the prior, it is assumed that conditional on
$\Sigma_{s_{t}}$, $s_{t}$, 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,s_{t}}$ $(\ell=0,\dots,p)$, it holds

$\mathbb{E}\big{(}(A_{\ell,s_{t}})_{i,j}\big{|}\Sigma_{s_{t}},s_{t},\mathcal{F}_{0}\big{)}=\begin{cases}\phi_{i}&\text{if}~{}~{}~{}i=j,~{}\ell=1,\\\
0&\text{if}~{}~{}~{}\text{otherwise}\end{cases},$ (109)
$\text{Var}\big{(}(A_{0,s_{t}})_{i,j}\big{|}\Sigma_{s_{t}},s_{t},\mathcal{F}_{0}\big{)}=(\sigma_{i,s_{t}}/\varepsilon_{s_{t}})^{2},$
(110)

and

$\text{Var}\big{(}(A_{\ell,s_{t}})_{ij}\big{|}\Sigma_{s_{t}},s_{t},\mathcal{F}_{0}\big{)}=\begin{cases}\displaystyle\bigg{(}\frac{\sigma_{i,s_{t}}}{\ell^{\lambda_{2,s_{t}}}\lambda_{1,s_{t}}\tau_{i,s_{t}}}\bigg{)}^{2}&\text{if}~{}~{}~{}i=j,\\\
\displaystyle\bigg{(}\frac{\sigma_{i,s_{t}}}{\ell^{\lambda_{2,s_{t}}}\lambda_{1,s_{t}}\tau_{j,s_{t}}}\bigg{)}^{2}&\text{if}~{}~{}~{}\text{otherwise}\end{cases}~{}~{}~{}\text{for}~{}\ell=1,\dots,p.$
(111)

The parameter $\varepsilon_{s_{t}}^{2}$ is a small number and it corresponds
to an uninformative diffuse prior for $(A_{0,s_{t}})_{i,j}$, the parameter
$\lambda_{1,s_{t}}$ controls the overall tightness of the prior distribution,
the parameter $\lambda_{2,s_{t}}$ controls amount of information prior
information at higher lags, and $\tau_{i,s_{t}}$ is a scaling parameter, see
Miranda18. Thus, the factor $1/\ell^{2\lambda_{2,s_{t}}}$ 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 (109), (110), and (111)

$\hat{y}_{s_{t}}^{\circ}=\Pi_{s_{t}}\hat{\mathsf{Y}}_{s_{t}}^{\circ}+\hat{\xi}^{\circ},$
(112)

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

$\hat{y}_{s_{t}}^{\circ}:=\big{[}0_{[n\times
l]}:\lambda_{1,s_{t}}\text{diag}\\{\phi_{1}\tau_{1,s_{t}},\dots,\phi_{n}\tau_{n,s_{t}}\\}:0_{[n\times(n-1)p]}\big{]}$
(113)

and

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

with
$J_{s_{t}}:=\text{diag}\\{1^{\lambda_{2,s_{t}}},\dots,p^{\lambda_{2,s_{t}}}\\}$,
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_{s_{t}}$ to the matrices of dummy variables,
see Miranda18. It is worth mentioning that the matrices of dummy variables
$\hat{y}_{s_{t}}$ and $\hat{\mathsf{Y}}_{s_{t}}$ should not depend on the
covariance matrix $\Sigma_{s_{t}}$. If the dummy variables depend on the
covariance matrix, an OLS estimator $\hat{\pi}_{s_{t}}$, and matrix
$\Lambda_{0,s_{t}}$ depend on the covariance matrix $\Sigma_{s_{t}}$, see
below. Consequently, in this case, one can not use the results of Proposition
2. For this reason, we choose the prior condition (111). Equation (112), can
be written by

$\hat{y}_{s_{t}}=\big{(}(\hat{\mathsf{Y}}_{s_{t}}^{\circ})^{\prime}\otimes
I_{n}\big{)}\pi_{s_{t}}+\hat{\xi},$ (115)

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

$\pi_{s_{t}}\overset{d}{=}\big{(}((\hat{\mathsf{Y}}_{s_{t}}^{\circ}(\hat{\mathsf{Y}}_{s_{t}}^{\circ})^{\prime})^{-1}\hat{\mathsf{Y}}_{s_{t}}^{\circ})\otimes
I_{n}\big{)}\hat{y}_{s_{t}}+\big{(}((\hat{\mathsf{Y}}_{s_{t}}^{\circ}(\hat{\mathsf{Y}}_{s_{t}}^{\circ})^{\prime})^{-1}\hat{\mathsf{Y}}_{s_{t}}^{\circ})\otimes
I_{n}\big{)}\hat{\xi},$ (116)

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_{s_{t}}$, namely,
$\pi_{0,s_{t}}:=\text{vec}\big{(}\hat{\Pi}_{s_{t}}\big{)}=\text{vec}\big{(}\hat{y}_{s_{t}}^{\circ}(\hat{\mathsf{Y}}_{s_{t}}^{\circ})^{\prime}(\hat{\mathsf{Y}}_{s_{t}}^{\circ}(\hat{\mathsf{Y}}_{s_{t}}^{\circ})^{\prime})^{-1}\big{)}$.
Consequently, conditional on $\Sigma_{s_{t}}$, $s_{t}$, and $\mathcal{F}_{0}$,
a distribution of the coefficient vector $\pi_{s_{t}}$ is given by

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

where
$\Lambda_{0,s_{t}}:=(\hat{\mathsf{Y}}_{s_{t}}(\hat{\mathsf{Y}}_{s_{t}}^{\circ})^{\prime})^{-1}$
is a $(d\times d)$ diagonal matrix. Consequently, conditional on
$\Sigma_{s_{t}}$, $s_{t}$, and $\mathcal{F}_{0}$, columns of the random
coefficient matrix $\Pi_{s_{t}}$ are independent. For the moment conditions
(109), (110), and (111), we can use the results of Proposition 2.

###  3.2 Simulation Methods

By applying Propositions 2 and 3, one can obtain exact density function of the
random vector $\bar{y}_{t}$, namely,

$\displaystyle
f(\bar{y}_{t}|\mathcal{F}_{0})=\sum_{\bar{s}_{t}}f(\bar{y}_{t}|\bar{s}_{t},\mathcal{F}_{0})\prod_{i=0}^{N}\frac{\Gamma\big{(}\sum_{j=1}^{N}\alpha_{ij}\big{)}}{\prod_{j=1}^{N}\Gamma(\alpha_{ij})}\frac{\prod_{j=1}^{N}\Gamma(\alpha_{ij}+n_{ij}(\bar{s}_{t}))}{\Gamma\big{(}\sum_{j=1}^{N}(\alpha_{ij}+n_{ij}(\bar{s}_{t}))\big{)}}.$
(118)

However, this exact density function has the following two main disadvantages:

  * $(i)$

it is difficult to obtain characteristics (such as mean, quantile, marginal
densities, distribution, and so on) of the mixture density function
$f(\bar{y}_{t}|\mathcal{F}_{0})$,

  * $(ii)$

and it is difficult to calculate the sum with respect to $\bar{s}_{t}$. For
example, if the length of the regime vector $\bar{s}_{t}$ equals 30 and the
regime number equals 3, then we have to calculate $3^{30}\approx 2.06\times
10^{14}$ summands.

If the dimensions increase, the disadvantages are seriously worsen. Therefore,
from a practical point of view, we need to develop Monte–Carlo simulation
method.

####  3.2.1 General Method

According to the conditional probability formula and Proposition 1, we have
that

$\displaystyle
f(\bar{y}_{t}^{c},\pi_{\hat{s}},\Sigma_{\hat{s}},s,\mathsf{P}|\mathcal{F}_{t})$
$\displaystyle=$ $\displaystyle
f(\bar{y}_{t}^{c}|\pi_{\beta},\Sigma_{\beta},s,\mathcal{F}_{t})f_{*}(\pi_{\delta}|\Sigma_{\delta},\delta,\mathcal{F}_{0})f_{*}(\Sigma_{\delta}|\delta,\mathcal{F}_{0})$
(119) $\displaystyle\times$ $\displaystyle
f(\pi_{\alpha}|\Sigma_{\alpha},\bar{s}_{t},\mathcal{F}_{t})f(\Sigma_{\alpha}|\bar{s}_{t},\mathcal{F}_{t})f(\bar{s}_{t}^{c}|s_{t},\mathsf{P},\mathcal{F}_{0})f(\bar{s}_{t},\mathsf{P}|\mathcal{F}_{t}).$

The above equation tells us that how to generate random sample from
$(\bar{y}_{t}^{c},\pi_{\hat{s}_{t}},\Sigma_{\hat{s}_{t}},s,\mathsf{P})$ for
given information $\mathcal{F}_{t}$. The direction of our simulation method
move toward from right to left for the above equation. To generate random
samples, first, generate the regime vector and transition probability matrix
$(\bar{s}_{t},\mathsf{P})$ from the posterior density function
$f(\bar{s}_{t},\mathsf{P}|\mathcal{F}_{t})$. Next, using the regime vector
$\bar{s}_{t}$ and transition probability matrix $\mathsf{P}$, generate regime
vector $\bar{s}_{t}^{c}$ from the conditional density function
$f(\bar{s}_{t}^{c}|\bar{s}_{t},\mathsf{P},\mathcal{F}_{0})$, so on.

First, we consider a simulation method that generate the regime vector and
transition probability matrix $(\bar{s}_{t},\mathsf{P})$ from the posterior
density function $f(\bar{s}_{t},\mathsf{P}|\mathcal{F}_{t})$. As mentioned
above, for given $\mathcal{F}_{0}$, the regime–switching process $s_{t}$ is
not a Markov chain. Therefore, we develop the Gibbs sampling method to
generate $(\bar{s}_{t},\mathsf{P})$. In the Bayesian statistics, the Gibbs
sampling is often used when the joint distribution is not known explicitly or
is difficult to sample from directly, but the conditional distribution of each
variable is known and is easy to sample from. Constructing the Gibbs sampler
to approximate the joint posterior distribution
$f(\bar{s}_{t},\mathsf{P}|\mathcal{F}_{t})$ is straightforward: New values
$\big{(}\bar{s}_{t}(\ell),\mathsf{P}(\ell)\big{)}$, $\ell=1,\dots,\mathcal{L}$
can be generated by

  * •

generate $\mathsf{P}(\ell)$ from
$f\left(\mathsf{P}|\bar{s}_{t}(\ell),\mathcal{F}_{t}\right)$,

  * •

generate $\bar{s}_{t}(\ell)$ from
$f\big{(}\bar{s}_{t}\big{|}\mathsf{P}(\ell),\mathcal{F}_{t}\big{)}$.

To generate $\mathsf{P}(\ell)$ from
$f(\mathsf{P}|\bar{s}_{t}(\ell),\mathcal{F}_{t})$, we apply equation (17) in
Proposition 1, density function (83), equation (84) in Proposition 3, and the
Bayesian formula. Then, we have that

$f(\mathsf{P}|\bar{s}_{t}(\ell),\mathcal{F}_{t})=\prod_{i=0}^{N}\frac{\Gamma\big{(}\sum_{j=1}^{N}(\alpha_{ij}+n_{ij}(\bar{s}_{t}(\ell)))\big{)}}{\prod_{j=1}^{N}\Gamma(\alpha_{ij}+n_{ij}(\bar{s}_{t}(\ell)))}\prod_{j=1}^{N}p_{ij}^{\alpha_{ij}+n_{ij}(\bar{s}_{t}(\ell))-1}.$
(120)

Thus, one can deduce that conditional on $\bar{s}_{t}(\ell)$ and
$\mathcal{F}_{t}$, for $i=0,\dots,N$, each row of the transition probability
matrix $\mathsf{P}$ are independent and has Dirichlet distribution with
parameter
$\alpha_{i}(\bar{s}_{t}(\ell)):=(\alpha_{i1}+n_{i1}(\bar{s}_{t}(\ell)),\dots,\alpha_{iN}+n_{iN}(\bar{s}_{t}(\ell)))^{\prime}$.
Consequently, the simulation method that generate $\mathsf{P}(\ell)$ for given
$\bar{s}_{t}(\ell)$ and $\mathcal{F}_{t}$ as follows:

  * •

generate $\mathsf{P}_{i}(\ell)$ from
$\text{Dir}(\alpha_{i}(\bar{s}_{t}(\ell)))$ for $i=0,\dots,N.$

Collect $\mathsf{P}_{i}(\ell)$ for $i=0,\dots,N$ into an $([N+1]\times N)$
matrix $\mathsf{P}(\ell)$, that is,
$\mathsf{P}(\ell):=[\mathsf{P}_{0}(\ell)^{\prime}:\dots:\mathsf{P}_{N}(\ell)^{\prime}]^{\prime}.$
Let
$\hat{\mathsf{P}}(\ell):=[\mathsf{P}_{1}(\ell)^{\prime}:\dots:\mathsf{P}_{N}(\ell)^{\prime}]^{\prime}$
be an $(N\times N)$ transition probability matrix, which omits first row of
the matrix $\mathsf{P}(\ell)$.

Now we consider a sampling method that generate $\bar{s}_{t}(\ell)$ from
$f\big{(}\bar{s}_{t}\big{|}\mathsf{P}(\ell),\mathcal{F}_{t}\big{)}$. Here we
follow the book of Hamilton94, see also Battulga22b and Battulga23a. If we
assume that the regime–switching process in regime $j$ at time $u$, then
according to equation (81), the conditional density function of the random
vector $y_{u}$ is given by the following equation

$\displaystyle\eta_{u,j}$ $\displaystyle:=$ $\displaystyle
f(y_{u}|s_{u}=j,\mathcal{F}_{u-1})$ (121) $\displaystyle=$
$\displaystyle\frac{1}{\pi^{n/2}}\frac{|\Lambda_{0,j}|^{-n/2}\Gamma_{n}\big{(}(\nu_{0,j}+1)/2\big{)}|V_{0,j}|^{\nu_{0,j}/2}}{\big{(}1+\mathsf{Y}_{u}^{\prime}\Lambda_{0,j}\mathsf{Y}_{u}\big{)}^{n/2}\Gamma_{n}(\nu_{0,j}/2)\big{|}B_{u,j}+V_{0,j}\big{|}^{(\nu_{0,j}+1)/2}},$

for $u=1,\dots,t$ and $j=1,\dots,N$, where $B_{u,j}$ is given by equation
(82). For all $u=1,\dots,t$, we collect the conditional density functions of
$y_{u}$ into an $(n\times 1)$ vector $\eta_{t}$, that is,
$\eta_{u}:=(\eta_{u,1},\dots,\eta_{u,N})^{\prime}$. Let us denote a
probabilistic inference about the value of the regime–switching process
$s_{u}$ equals to $j$, based on the information $\mathcal{F}_{u}$ and
transition probability matrix $\mathsf{P}(\ell)$ by
$\mathbb{P}(s_{u}=j|\mathsf{P}(\ell),\mathcal{F}_{u})$. Collect these
conditional probabilities
$\mathbb{P}(s_{u}=j|\mathsf{P}(\ell),\mathcal{F}_{u})$ for $j=1,\dots,N$ into
an $(N\times 1)$ vector $z_{u|u}(\ell)$, that is,
$z_{u|u}(\ell):=\big{(}\mathbb{P}(s_{u}=1|\mathsf{P}(\ell),\mathcal{F}_{u}),\dots,\mathbb{P}(s_{u}=N|\mathsf{P}(\ell),\mathcal{F}_{u})\big{)}^{\prime}$.
Also, we need a probabilistic forecast about the value of the regime–switching
process at time $u+1$ equals $j$ conditional on data up to and including time
$u$ and transition probability matrix $\mathsf{P}(\ell)$. Collect these
forecasts into an $(N\times 1)$ vector $z_{u+1|u}(\ell)$, that is,
$z_{u+1|u}(\ell):=\big{(}\mathbb{P}(s_{u+1}=1|\mathsf{P}(\ell),\mathcal{F}_{u}),\dots,\mathbb{P}(s_{u+1}=N|\mathsf{P}(\ell),\mathcal{F}_{u})\big{)}^{\prime}$.

The probabilistic inference and forecast for each time $u=1,\dots,t$ can be
found by iterating on the following pair of equations:

$z_{u|u}(\ell)=\frac{(z_{u|u-1}(\ell)\odot\eta_{u})}{i_{N}^{\prime}(z_{u|u-1}(\ell)\odot\eta_{u})}~{}~{}~{}\text{and}~{}~{}~{}z_{u+1|u}(\ell)=\hat{\mathsf{P}}(\ell)^{\prime}z_{u|u}(\ell),~{}~{}~{}u=1,\dots,t,$
(122)

where $\odot$ is the Hadamard product of two vectors, $\eta_{u}$ is the
$(N\times 1)$ vector, whose $j$-th element is given by equation (121),
$\hat{\mathsf{P}}(\ell)$ is the $(N\times N)$ transition probability matrix,
and $i_{N}$ is an $(N\times 1)$ vector, whose elements equal 1. Given a
starting value $z_{1|0}(\ell):=\mathsf{P}_{0}(\ell)^{\prime}$ one can iterate
on (122) for $u=1,\dots,t$ to calculate the values of $z_{u|u}(\ell)$ and
$z_{u+1|u}(\ell)$. To obtain marginal distributions of the regime vector
$\bar{s}_{t}$ conditional on the transition probability matrix
$\mathsf{P}(\ell)$ and information $\mathcal{F}_{t}$, let us introduce
$(N\times 1)$ smoothed inference vector
$z_{u|t}(\ell):=\big{(}\mathbb{P}(s_{u}=1|\mathsf{P}(\ell),\mathcal{F}_{t}),\dots,\mathbb{P}(s_{u}=N|\mathsf{P}(\ell),\mathcal{F}_{t})\big{)}^{\prime}$
for $u=1,\dots,t$. The smoothed inference vectors can be obtained by using the
Kim94’s smoothing algorithm:

$z_{u|t}(\ell)=z_{u|u}(\ell)\odot\big{\\{}\hat{\mathsf{P}}(\ell)^{\prime}(z_{u+1|t}(\ell)\oslash
z_{u+1|u}(\ell))\big{\\}},~{}~{}~{}u=t-1,\dots,1,$ (123)

where $\oslash$ is an element–wise division of two vectors. The smoothed
probabilities $z_{u|t}(\ell)$ are found by iterating on (123) backward for
$u=t-1,\dots,1$. This iteration is started with $z_{t|t}(\ell)$, which is
obtained from (122) for $u=t$. Thus, the simulation method that generate
$\bar{s}_{t}(\ell)$ for given $\mathsf{P}(\ell)$ and $\mathcal{F}_{t}$ as
follows:

  * •

generate $s_{u}(\ell)$ from $z_{u|t}(\ell)$ for $u=1,\dots,t.$

Collect $s_{u}(\ell)$ for $u=1,\dots,t$ into $(t\times 1)$ vector
$\bar{s}_{t}(\ell)$, namely,
$\bar{s}_{t}(\ell):=(s_{1}(\ell),\dots,s_{t}(\ell))^{\prime}$.

Second, we consider a simulation method that generate the regime vector
$\bar{s}_{t}^{c}(\ell)$ from the density function
$f(\bar{s}_{t}^{c}|s_{t}(\ell),\mathsf{P}(\ell),\mathcal{F}_{0})$. Note that
conditional on the transition probability matrix $\mathsf{P}(\ell)$ and
initial information $\mathcal{F}_{0}$, the regime–switching process $s_{t}$ is
a Markov chain. Thus, to sample the regime vector $\bar{s}_{t}^{c}(\ell)$ from
the density function
$f(\bar{s}_{t}^{c}|s_{t}(\ell),\mathsf{P}(\ell),\mathcal{F}_{0})$, we can use
the Markov property. That is, we have that

$f(\bar{s}_{t}^{c}|s_{t}(\ell),\mathsf{P}(\ell),\mathcal{F}_{0})=f(s_{T}|s_{T-1},\mathsf{P}(\ell),\mathcal{F}_{0})\dots
f(s_{t+1}|s_{t}(\ell),\mathsf{P}(\ell),\mathcal{F}_{0}).$ (124)

Thus, a simulation method that generate $\bar{s}_{t}^{c}(\ell)$ for given
$s_{t}(\ell)$, $\mathsf{P}(\ell)$, and $\mathcal{F}_{0}$ as follows:

  * •

generate $s_{t+1}(\ell)$ from
$f(s_{t+1}|s_{t}(\ell),\mathsf{P}(\ell),\mathcal{F}_{0})$,

  * •

generate $s_{t+2}(\ell)$ from
$f(s_{t+2}|s_{t+1}(\ell),\mathsf{P}(\ell),\mathcal{F}_{0})$,

  * $\dots$

  * •

generate $s_{T}(\ell)$ from
$f(s_{T}|s_{T-1}(\ell),\mathsf{P}(\ell),\mathcal{F}_{0})$.

Collect $s_{u}(\ell)$ for $u=t+1,\dots,T$ into $([T-t]\times 1)$ vector
$\bar{s}_{t}^{c}(\ell)$, namely,
$\bar{s}_{t}^{c}(\ell):=(s_{t+1}(\ell),\dots,s_{T}(\ell))^{\prime}$.

Third, we consider a simulation method that generate the coefficient vector
$\pi_{\bar{s}_{t}^{c}(\ell)}(\ell)$ and covariance matrix
$\Sigma_{\bar{s}_{t}^{c}(\ell)}(\ell)$ from the density function
$f\big{(}\pi_{\bar{s}_{t}^{c}(\ell)},\Sigma_{\bar{s}_{t}^{c}(\ell)}\big{|}s(\ell),\mathcal{F}_{t}\big{)}$.
Let
$\mathcal{A}_{{\bar{s}}_{t}(\ell)}=\big{\\{}\alpha_{1}(\ell),\dots,\alpha_{r_{\alpha(\ell)}}(\ell)\big{\\}}$
and
$\mathcal{A}_{{\bar{s}}_{t}^{c}(\ell)}=\big{\\{}\beta_{1}(\ell),\dots,\beta_{r_{\beta(\ell)}}(\ell)\big{\\}}$
be the duplication removed sets, corresponding the regime vectors
$\bar{s}_{t}(\ell)$ and $\bar{s}_{t}^{c}(\ell)$, respectively. To eliminate
unnecessary simulations, instead of the regimes in the set
$\mathcal{A}_{{\bar{s}}_{t}(\ell)}$, one should consider regimes in
intersection set of the sets $\mathcal{A}_{{\bar{s}}_{t}(\ell)}$ and
$\mathcal{A}_{{\bar{s}}_{t}^{c}(\ell)}$. Let us assume that the intersection
set and a difference set of the sets are given by
$\mathcal{A}_{{\bar{s}}_{t}(\ell)}\cap\mathcal{A}_{{\bar{s}}_{t}^{c}(\ell)}=\big{\\{}\gamma_{1}(\ell),\dots,\gamma_{r_{\gamma(\ell)}}(\ell)\big{\\}}$
and
$\mathcal{A}_{{\bar{s}}_{t}^{c}(\ell)}\backslash\mathcal{A}_{{\bar{s}}_{t}(\ell)}=\big{\\{}\delta_{1}(\ell),\dots,\delta_{k}(\ell)\big{\\}}$,
respectively. Then, according to equations (10), (73) and (74), a simulation
method that generates
$\big{(}\pi_{\gamma_{k}(\ell)}(\ell),\Sigma_{\gamma_{k}(\ell)}(\ell)\big{)}$
for $\ell=1,\dots,r_{\gamma(\ell)}$ as follows: for
$k=1,\dots,r_{\gamma(\ell)}$,

  * •

generate $\Sigma_{\gamma_{k}(\ell)}(\ell)$ from
$\mathcal{IW}\Big{(}\nu_{0,\gamma_{k}(\ell)}+q_{t,\gamma_{k}(\ell)},\bar{B}_{t,\gamma_{k}(\ell)}+V_{0,\gamma_{k}(\ell)}\Big{)},$

  * •

generate $\pi_{\gamma_{k}(\ell)}(\ell)$ from
$\mathcal{N}\Big{(}A_{\gamma_{k}(\ell)}b_{\gamma_{k}(\ell)},A_{\gamma_{k}(\ell)}\Big{)}$.

On the other hand, according to equations (10) and (15) and marginal density
functions of equations (3.1.1) and (51), a simulation method that generates
$\big{(}\pi_{\delta_{k}(\ell)}(\ell),\Sigma_{\delta_{k}(\ell)}(\ell)\big{)}$
for $\ell=1,\dots,r_{\delta(\ell)}$ as follows: if $r_{\delta(\ell)}>0$, then
for $k=1,\dots,r_{\delta(\ell)}$,

  * •

generate $\Sigma_{\delta_{k}(\ell)}(\ell)$ from
$\mathcal{IW}\Big{(}\nu_{0,\delta_{k}(\ell)},V_{0,\delta_{k}(\ell)}\Big{)}$,

  * •

generate $\pi_{\delta_{k}(\ell)}(\ell)$ from
$\mathcal{N}\Big{(}\pi_{0,\delta_{k}(\ell)},\Lambda_{0,\delta_{k}(\ell)}\otimes\Sigma_{\delta_{k}(\ell)}(\ell)\Big{)}$.

Let
$\hat{\beta}(\ell):=\big{(}\gamma_{1}(\ell),\dots,\gamma_{r_{\gamma}(\ell)}(\ell),\delta_{1}(\ell),\dots,\delta_{r_{\delta(\ell)}}(\ell)\big{)}^{\prime}$
be an $(r_{\beta(\ell)}\times 1)$ regime vector. The regime vector
$\hat{\beta}(\ell)$ has same elements as the duplication removed regime vector
$\beta(\ell)$, but positions of the elements are different for the two regime
vectors. Collect the realizations for $k=1,\dots,r_{\gamma(\ell)}$,
$\Sigma_{\gamma_{k}(\ell)}(\ell)$ and $\pi_{\gamma_{k}(\ell)}(\ell)$ and for
$k=1,\dots,r_{\delta(\ell)}$, $\Sigma_{\delta_{k}(\ell)}(\ell)$ and
$\pi_{\delta_{k}(\ell)}(\ell)$ into $([ndr_{\beta(\ell)}]\times 1)$ vector
$\pi_{\hat{\beta}(\ell)}(\ell)$ and $(n\times[ndr_{\beta(\ell)}])$ matrix
$\Sigma_{\hat{\beta}(\ell)}(\ell)$, namely,

$\pi_{\hat{\beta}(\ell)}(\ell):=\Big{(}\pi_{\gamma_{1}(\ell)}(\ell)^{\prime},\dots,\pi_{\gamma_{r_{\gamma(\ell)}}(\ell)}(\ell)^{\prime},\pi_{\delta_{1}(\ell)}(\ell)^{\prime},\dots,\pi_{\delta_{r_{\delta(\ell)}}(\ell)}(\ell)^{\prime}\Big{)}^{\prime}$
(125)

and

$\Sigma_{\hat{\beta}(\ell)}(\ell):=\Big{[}\Sigma_{\gamma_{1}(\ell)}(\ell):\dots:\Sigma_{\gamma_{r_{\gamma(\ell)}}(\ell)}(\ell),\Sigma_{\delta_{1}(\ell)}(\ell):\dots:\Sigma_{\delta_{r_{\delta(\ell)}}(\ell)}(\ell)\Big{]}.$
(126)

Similar to equation (37), for $u=t+1,\dots,T$, we denote a position of the
regime $s_{u}(\ell)$ in the regime vector $\hat{\beta}(\ell)$ by
$o_{\ell}(\ell)$. Let us define a matrix
$D_{\hat{\beta}(\ell)}:=\big{[}j_{o_{t+1}(\ell)}:\dots:j_{o_{T}(\ell)}\big{]}^{\prime}$,
where $j_{o}(\ell)$ is an ($r_{\beta(\ell)}\times 1$) unit vector, whose
$o$–th element 1 and others 0. Then, one revives the vector
$\pi_{\bar{s}_{t}^{c}(\ell)}=\big{(}D_{\hat{\beta}(\ell)}\otimes
I_{nd}\big{)}\pi_{\hat{\beta}(\ell)}$ and matrix

$\Sigma_{\bar{s}_{t}^{c}(\ell)}=\text{diag}\big{\\{}\big{(}(D_{\hat{\beta}(\ell)}\otimes
I_{n})\Sigma_{\hat{\beta}(\ell)}\big{)}_{1},\dots,\big{(}(D_{\hat{\beta}(\ell)}\otimes
I_{n})\Sigma_{\hat{\beta}(\ell)}\big{)}_{T}\big{\\}}.$ (127)

Fourth, we consider a simulation method that generate the regime vector
$\bar{y}_{t}^{c}(\ell)$ from the density function
$f(\bar{y}_{t}^{c}|\pi_{\beta},\Sigma_{\beta},s,\mathcal{F}_{t})$. Let us
assume that

$\varphi_{s}=\begin{bmatrix}\varphi_{1,\bar{s}_{t}}\\\
\varphi_{2,\bar{s}_{t}^{c}}\end{bmatrix}~{}~{}~{}\text{and}~{}~{}~{}\Psi_{s}=\begin{bmatrix}\Psi_{11,\bar{s}_{t}}&0\\\
\Psi_{21,\bar{s}_{t}^{c}}&\Psi_{22,\bar{s}_{t}^{c}}\end{bmatrix}$ (128)

are partitions of the vector $\varphi_{s}$ and matrix $\Psi_{s}$,
corresponding to random sub vectors $\bar{y}_{t}$ and $\bar{y}_{t}^{c}$ of the
random vector $y$. Then, due to Battulga24a, a distribution of the random
vector $\bar{y}_{t}^{c}$ is given by

$\displaystyle\bar{y}_{t}^{c}~{}|~{}\pi_{\beta},\Sigma_{\beta},s,\mathcal{F}_{t}$
$\displaystyle\sim$
$\displaystyle\mathcal{N}\Big{(}\Psi_{22,\bar{s}_{t}^{c}}^{-1}\big{(}\varphi_{2,\bar{s}_{t}^{c}}-\Psi_{21,\bar{s}_{t}^{c}}\bar{y}_{t}\big{)},\Psi_{22,\bar{s}_{t}^{c}}^{-1}\Sigma_{\bar{s}_{t}^{c}}(\Psi_{22,\bar{s}_{t}^{c}}^{-1})^{\prime}\Big{)},$
(129)

where
$\Sigma_{\bar{s}_{t}^{c}}=\text{diag}\\{\Sigma_{s_{t+1}},\dots,\Sigma_{s_{T}}\\}$
is an $([n(T-t)]\times[n(T-t)])$ block diagonal matrix, corresponding to the
regime vector $\bar{s}_{t}^{c}$. Thus, a simulation method that generate the
vector of endogenous variables $\bar{y}_{t}^{c}(\ell)$ for given
$\pi_{\bar{s}_{t}(\ell)}(\ell)$, $\Sigma_{\bar{s}_{t}(\ell)}(\ell)$,
$s(\ell)$, and $\mathcal{F}_{t}$ as follows:

  * •

generate $\bar{y}_{t}^{c}(\ell)$ from
$\mathcal{N}\Big{(}\Psi_{22,\bar{s}_{t}^{c}(\ell)}^{-1}\big{(}\varphi_{2,\bar{s}_{t}^{c}(\ell)}-\Psi_{21,\bar{s}_{t}^{c}(\ell)}\bar{y}_{t}\big{)},\Psi_{22,\bar{s}_{t}^{c}(\ell)}^{-1}\Sigma_{\bar{s}_{t}^{c}(\ell)}(\ell)\big{(}\Psi_{22,\bar{s}_{t}^{c}(\ell)}^{-1}\big{)}^{\prime}\Big{)},$

where the matrix $\Psi_{22,\bar{s}_{t}^{c}(\ell)}$ and vector
$\varphi_{\bar{s}_{t}^{c}(\ell)}-\Psi_{21,\bar{s}_{t}^{c}(\ell)}\bar{y}_{t}$
are given by

$\Psi_{22,\bar{s}_{t}^{c}(\ell)}:=\begin{bmatrix}I_{n}&0&\dots&0&\dots&0&0\\\
-A_{1,s_{t+2}(\ell)}&I_{n}&\dots&0&\dots&0&0\\\
\vdots&\vdots&\dots&\vdots&\dots&\vdots&\vdots\\\
0&0&\dots&-A_{p-1,s_{T-1}(\ell)}&\dots&I_{n}&0\\\
0&0&\dots&-A_{p,s_{T}(\ell)}&\dots&-A_{1,s_{T}(\ell)}&I_{n}\end{bmatrix}$
(130)

and

$\varphi_{2,\bar{s}_{t}^{c}(\ell)}-\Psi_{21,\bar{s}_{t}^{c}(\ell)}\bar{y}_{t}:=\begin{bmatrix}A_{0,s_{t+1}(\ell)}\psi_{t+1}+A_{1,s_{t+1}(\ell)}y_{t}+\dots+A_{p,s_{t+1}(\ell)}y_{t+1-p}\\\
A_{0,s_{t+2}(\ell)}\psi_{t+2}+A_{2,s_{t+2}(\ell)}y_{t}+\dots+A_{p,s_{t+2}(\ell)}y_{t+2-p}\\\
\vdots\\\ A_{0,s_{T-1}(\ell)}\psi_{T-1}\\\
A_{0,s_{T}(\ell)}\psi_{T}\end{bmatrix}$ (131)

and they are obtained from the vector $\pi_{\bar{s}_{t}^{c}(\ell)}(\ell)$. It
should be noted that traditional methods that generate the vector
$\bar{y}_{t}^{c}$ are based on an iterative method for $y_{t+1},\dots,y_{T}$
by generating $\xi_{t+1},\dots,\xi_{T}$, see Karlsson13. As a result, if $T-t$
is large, the simulation method reduces the computational burden that
generates the random vector $\bar{y}_{t}^{c}$ as compared to the traditional
algorithms.

####  3.2.2 Importance Sampling Method

Now, we consider the importance sampling method for the Bayesian MS–VAR
process. We estimate a probability of a rare event, corresponding to the
endogenous variables by the important sampling method. In the importance
sampling method, one changes the real probability measure $\mathbb{P}$. The
new probability measure $\tilde{\mathbb{P}}$ must be chosen that the rare
event more frequently comes from than the real probability measure
$\mathbb{P}$. Let
$f(y,\pi_{\hat{s}},\Sigma_{\hat{s}},s,\mathsf{P}|\mathcal{F}_{0})$ and
$\tilde{f}(y,\pi_{\hat{s}},\Sigma_{\hat{s}},s,\mathsf{P}|\mathcal{F}_{0})$ be
joint density functions under the real probability measure $\mathbb{P}$ and
new probability measure $\tilde{\mathbb{P}}$, respectively, for given initial
information $\mathcal{F}_{0}$ and

$L=L(y,\pi_{\hat{s}},\Sigma_{\hat{s}},s,\mathsf{P}|\mathcal{F}_{0})=f(y,\pi_{\hat{s}},\Sigma_{\hat{s}},s,\mathsf{P}|\mathcal{F}_{0})/\tilde{f}(y,\pi_{\hat{s}},\Sigma_{\hat{s}},s,\mathsf{P}|\mathcal{F}_{0})$
(132)

be the likelihood ratio. Let us choose density functions that corresponds to
the new probability measure $\mathbb{\tilde{P}}$ by

$\tilde{f}(\bar{y}_{t},\Sigma_{\hat{s}},s,\mathsf{P}|\mathcal{F}_{0})=f(\bar{y}_{t},\Sigma_{\hat{s}},s,\mathsf{P}|\mathcal{F}_{0}),$
(133)

if $s_{u}\in\gamma$, then

$\displaystyle\tilde{f}(\pi_{\hat{s}}|\Sigma_{\hat{s}},s,\mathsf{P},\mathcal{F}_{t})=f(\pi_{\hat{s}}|\Sigma_{\hat{s}},s,\mathsf{P},\mathcal{F}_{t})=f(\pi_{\hat{s}}|\Sigma_{\hat{s}},s,\mathcal{F}_{t})$
$\displaystyle\Big{/}\frac{1}{(2\pi)^{nd}|A_{s_{u}}|^{1/2}}\exp\bigg{\\{}-\frac{1}{2}\Big{(}\pi_{s_{u}}-A_{s_{u}}b_{s_{u}}\Big{)}^{\prime}A_{s_{u}}^{-1}\Big{(}\pi_{s_{u}}-A_{s_{u}}b_{s_{u}}\Big{)}\bigg{\\}}$
$\displaystyle\times\frac{1}{(2\pi)^{nd}|A_{s_{u}}|^{1/2}}\exp\bigg{\\{}-\frac{1}{2}\Big{(}\pi_{s_{u}}-A_{s_{u}}b_{s_{u}}-\theta_{s_{u}}A_{s_{u}}(\mathsf{Y}_{u}\otimes
I_{n})z_{u}\Big{)}^{\prime}$ (134) $\displaystyle\times
A_{s_{u}}^{-1}\Big{(}\pi_{s_{u}}-A_{s_{u}}b_{s_{u}}-\theta_{s_{u}}A_{s_{u}}(\mathsf{Y}_{u}\otimes
I_{n})z_{u}\Big{)}\bigg{\\}}$

if $s_{u}\in\delta$, then

$\displaystyle\tilde{f}(\pi_{\hat{s}}|\Sigma_{\hat{s}},s,\mathsf{P},\mathcal{F}_{t})=f(\pi_{\hat{s}}|\Sigma_{\hat{s}},s,\mathsf{P},\mathcal{F}_{t})=f(\pi_{\hat{s}}|\Sigma_{\hat{s}},s,\mathcal{F}_{t})$
$\displaystyle\Big{/}\frac{1}{(2\pi)^{nd}|\Lambda_{0,s_{u}}|^{n/2}|\Sigma_{s_{u}}|^{d/2}}\exp\bigg{\\{}-\frac{1}{2}\Big{(}\pi_{s_{u}}-\pi_{0,s_{u}}\Big{)}^{\prime}\big{(}\Lambda_{0,s_{u}}^{-1}\otimes\Sigma_{s_{u}}^{-1}\big{)}\Big{(}\pi_{s_{u}}-\pi_{0,s_{u}}\Big{)}\bigg{\\}}$
$\displaystyle\times\frac{1}{(2\pi)^{nd}|\Lambda_{0,s_{u}}|^{n/2}|\Sigma_{s_{u}}|^{d/2}}\exp\bigg{\\{}-\frac{1}{2}\Big{(}\pi_{s_{u}}-\pi_{0,s_{u}}-\theta_{s_{u}}(\Lambda_{0,s_{u}}\mathsf{Y}_{u}\otimes\Sigma_{s_{u}})z_{u}\Big{)}^{\prime}$
$\displaystyle\times\big{(}\Lambda_{0,s_{u}}^{-1}\otimes\Sigma_{s_{u}}^{-1}\big{)}\Big{(}\pi_{s_{u}}-\pi_{0,s_{u}}-\theta_{s_{u}}(\Lambda_{0,s_{u}}\mathsf{Y}_{u}\otimes\Sigma_{s_{u}})z_{u}\Big{)}\bigg{\\}},$

and

$\displaystyle\tilde{f}(\bar{y}_{t}^{c}|\pi_{\hat{s}},\Sigma_{\hat{s}},s,\mathsf{P},\mathcal{F}_{t})=f(\bar{y}_{t}^{c}|\pi_{\hat{s}},\Sigma_{\hat{s}},s,\mathsf{P},\mathcal{F}_{t})=f(\bar{y}_{t}^{c}|\pi_{\hat{s}},\Sigma_{\hat{s}},s,\mathcal{F}_{t})$
$\displaystyle\Big{/}\frac{1}{(2\pi)^{n}|\Sigma_{s_{u}}|^{1/2}}\exp\bigg{\\{}-\frac{1}{2}\Big{(}y_{u}-\Pi_{s_{u}}\mathsf{Y}_{u}\Big{)}^{\prime}\Sigma_{s_{u}}^{-1}\Big{(}y_{u}-\Pi_{s_{u}}\mathsf{Y}_{u}\Big{)}\bigg{\\}}$
(136)
$\displaystyle\times\frac{1}{(2\pi)^{n}|\Sigma_{s_{u}}|^{1/2}}\exp\bigg{\\{}-\frac{1}{2}\Big{(}y_{u}-\Pi_{s_{u}}\mathsf{Y}_{u}-\theta_{s_{u}}\Sigma_{s_{u}}z_{u}\Big{)}^{\prime}\Sigma_{s_{u}}^{-1}\Big{(}y_{u}-\Pi_{s_{u}}\mathsf{Y}_{u}-\theta_{s_{u}}\Sigma_{s_{u}}z_{u}\Big{)}\bigg{\\}}$

for $u=t+1,\dots,T$, where $\theta_{s_{u}}$ is a positive constant, depending
on the random coefficient vector $\pi_{s_{u}}$, random covariance matrix
$\Sigma_{s_{t}}$, and regime $s_{u}$ and $z_{u}$ is an $(n\times 1)$ vector,
whose elements are known. Note that if $\theta_{s_{u}}=0$, then the new
probability measure $\tilde{\mathbb{P}}$ equals the real probability measure
$\mathbb{P}$. If we compare the density functions
$\tilde{f}(\bar{y}_{t}^{c}|\pi_{\hat{s}},\Sigma_{\hat{s}},s,\mathsf{P},\mathcal{F}_{t})$
and
$f(\bar{y}_{t}^{c}|\pi_{\hat{s}},\Sigma_{\hat{s}},s,\mathsf{P},\mathcal{F}_{t})$,
one can conclude that conditional distribution of $y_{u}$ changes from
$y_{u}|\pi_{s_{u}},s_{u},\mathcal{F}_{u-1}\sim\mathcal{N}(\Pi_{s_{u}}\mathsf{Y}_{u},\Sigma_{s_{u}})$
to
$y_{u}|\pi_{s_{u}},s_{u},\mathcal{F}_{u-1}\sim\mathcal{N}(\Pi_{s_{u}}\mathsf{Y}_{u}+\theta_{s_{u}}\Sigma_{s_{u}}z_{u},\Sigma_{s_{u}})$
and for each $v=t+1,\dots,T$ $(v\neq u)$, the conditional distribution of
other sub random vector $y_{v}$ of the random vector $\bar{y}_{t}^{c}$ does
not change. The same explanation holds for the density functions
$\tilde{f}(\pi_{\hat{s}}|\Sigma_{\hat{s}},s,\mathsf{P},\mathcal{F}_{t})$ and
$f(\pi_{\hat{s}}|\Sigma_{\hat{s}},s,\mathsf{P},\mathcal{F}_{t})$. Then, it can
be shown that for $u=t+1,\dots,T$, the likelihood ratio is given by

$L_{u}=\exp\big{\\{}-\theta_{s_{u}}X_{u}+\psi(\theta_{s_{u}})\big{\\}}.$ (137)

where the random variable $X_{u}$ is given by $X_{u}:=z_{u}^{\prime}y_{u}$ and
the quadratic function $\psi(\theta_{s_{u}})$ for $\theta_{s_{u}}$ is given by

$\psi(\theta_{s_{u}}):=\begin{cases}\theta_{u}z_{u}^{\prime}C_{u,s_{u}}\mathsf{Y}_{u}+\frac{1}{2}\theta_{u}^{2}\big{(}1+\mathsf{Y}_{u}^{\prime}\Lambda_{0,s_{u}|u}^{-1}\mathsf{Y}_{u}\big{)}z_{u}^{\prime}\Sigma_{s_{u}}z_{u}&\text{if}~{}~{}~{}s_{u}\in\gamma,\\\
\theta_{u}z_{u}^{\prime}\Pi_{s_{u}}\mathsf{Y}_{u}+\frac{1}{2}\theta_{u}^{2}\big{(}1+\mathsf{Y}_{u}^{\prime}\Lambda_{0,s_{u}}\mathsf{Y}_{u}\big{)}z_{u}^{\prime}\Sigma_{s_{u}}z_{u}&\text{if}~{}~{}~{}s_{u}\in\delta.\end{cases}$
(138)

For each $u=t+1,\dots,T$, by choosing $z_{u}$ by the unit vector, one can
extract elements of the vector of endogenous variables $y_{u}$. If the process
$y_{t}$ consists of returns of financial assets, then by choosing $z_{u}$ by
weight vector, one obtains portfolio return.

Now we consider a conditional probability
$\mathbb{P}(X_{u}>x_{u}|\pi_{\beta},\Sigma_{\beta},s,\mathcal{F}_{u-1})$ for
large $x_{u}\in\mathbb{R}$ and $u=t+1,\dots,T$. Since $\theta_{s_{u}}$ is the
positive constant, for the conditional probability, by equation (137), the
following inequality holds

$\mathbb{P}(X_{u}>x_{u}|\pi_{\beta},\Sigma_{\beta},s,\mathcal{F}_{u-1})=\mathbb{\tilde{E}}\big{[}1_{\\{X_{u}>x_{u}\\}}L_{u}\big{|}\pi_{\beta},\Sigma_{\beta},s,\mathcal{F}_{u-1}\big{]}\leq\exp\big{\\{}-\theta_{s_{u}}x_{u}+\psi(\theta_{s_{u}})\big{\\}}$
(139)

for $u=t+1,\dots,T$, where for a generic event $A\in\mathcal{H}_{T}$, $1_{A}$
is the indicator function of the event $A$, see Glasserman00. Also, for the
second order moment, it holds

$m_{2}(x_{u},\theta_{s_{u}}):=\mathbb{\tilde{E}}\big{[}1_{\\{X_{u}>x_{u}\\}}L_{u}^{2}\big{|}\pi_{\beta},\Sigma_{\beta},s,\mathcal{F}_{u-1}\big{]}\leq\exp\big{\\{}-2\theta_{s_{u}}x_{u}+2\psi(\theta_{s_{u}})\big{\\}}.$
(140)

To reduce a variance of the importance sampling, we need to keep the right
hand side of the above inequality as low as possible. To minimize the
right–hand side of the above equation, we minimize the exponent using the
parameter $\theta_{s_{t}}$. The minimizer of the right–hand side of the above
equation is given by

$\theta_{s_{u}}(x_{u}):=\begin{cases}\displaystyle\frac{x_{u}-z_{u}^{\prime}C_{u,s_{u}}\mathsf{Y}_{u}}{\big{(}1+\mathsf{Y}_{u}^{\prime}\Lambda_{0,s_{u}|u}^{-1}\mathsf{Y}_{u}\big{)}z_{u}^{\prime}\Sigma_{s_{u}}z_{u}}&\text{if}~{}~{}~{}s_{u}\in\gamma,\\\
\displaystyle\frac{x_{u}-z_{u}^{\prime}\Pi_{s_{u}}\mathsf{Y}_{u}}{\big{(}1+\mathsf{Y}_{u}^{\prime}\Lambda_{0,s_{u}}\mathsf{Y}_{u}\big{)}z_{u}^{\prime}\Sigma_{s_{u}}z_{u}}&\text{if}~{}~{}~{}s_{u}\in\delta\end{cases}$
(141)

for $u=t+1,\dots,T$. Therefore, an importance sampling method that estimates
the conditional probabilities $\mathbb{P}(X_{u}>x_{u}|\mathcal{F}_{t})$ for
$u=t+1,\dots,T$ as follows: for $\ell=1,\dots,\mathcal{L}$,

  * •

generate
$\big{(}\bar{y}_{t}^{c}(\ell),\pi_{\bar{s}_{t}^{c}(\ell)}(\ell),\Sigma_{\bar{s}_{t}^{c}(\ell)}(\ell),\bar{s}_{t}^{c}(\ell)\big{)}$
using the general simulation method,

  * •

for $u=t+1,\dots,T$,

    * –

calculate $\theta_{s_{u}(\ell)}^{*}(x_{u},\ell)$ using equation (141),

    * –

if $s_{u}(\ell)\in\gamma(\ell)$, generate $\pi_{s_{u}(\ell)}^{*}(\ell)$ from

$\mathcal{N}\Big{(}A_{s_{u}(\ell)}(\ell)b_{s_{u}(\ell)}(\ell)+\theta_{s_{u}(\ell)}^{*}(x_{u},\ell)(\mathsf{Y}_{u}(\ell)\otimes
I_{n})z_{u},A_{s_{u}}(\ell)\Big{)},$ (142)

    * –

if $s_{u}(\ell)\in\delta(\ell)$, generate $\pi_{s_{u}(\ell)}^{*}(\ell)$ from

$\mathcal{N}\Big{(}\pi_{0,s_{u}(\ell)}+\theta_{s_{u}(\ell)}^{*}(x_{u},\ell)(\Lambda_{0,s_{u}(\ell)}\mathsf{Y}_{u}(\ell)\otimes\Sigma_{s_{u}(\ell)}(\ell))z_{u},\Lambda_{0,s_{u}(\ell)}\otimes\Sigma_{s_{u}(\ell)}(\ell)\Big{)},$
(143)

    * –

generate $y_{u}^{*}(\ell)$ from

$\mathcal{N}\Big{(}\Pi_{s_{u}}^{*}\mathsf{Y}_{u}(\ell)+\theta_{s_{u}(\ell)}^{*}(x_{u},\ell)\Sigma_{s_{u}(\ell)}z_{u},\Sigma_{s_{u}(\ell)}(\ell)\Big{)},$
(144)

where $\pi_{s_{u}(\ell)}^{*}(\ell)=\text{vec}(\Pi_{s_{u}(\ell)}^{*}(\ell))$,

    * –

calculate $L_{u}^{*}(\ell)$ using equation (137), where
$X_{u}^{*}(\ell)=z_{u}^{\prime}y_{u}^{*}(\ell)$,

  * •

and for $u=t+1,\dots,T$, estimate the probabilities
$\mathbb{P}(X_{u}>x_{u}|\mathcal{F}_{t})$ by

$\hat{\mathbb{P}}(X_{u}>x_{u}|\mathcal{F}_{t})=\frac{1}{\mathcal{L}}\sum_{\ell=1}^{\mathcal{L}}1_{\\{X_{u}^{*}(\ell)>x_{u}\\}}L_{u}^{*}(\ell),$
(145)

where calculations of the matrix $A_{s_{u}(\ell)}(\ell)$ and vectors
$b_{s_{u}(\ell)}(\ell)$ and $\mathsf{Y}_{u}(\ell)$ are based on the vectors
$\bar{y}_{t}$ and $\bar{y}_{t}^{c}(\ell)$.

##  4 Numerical Results

Now we consider numerical results. For means of illustration, we have chosen
three companies, listed in the S&P 500 index from different sectors. In order
to increase the number of price and dividend observation points, we take
quarterly data instead of yearly data. Our data covers a period from Q1, 1985
to Q4, 2023. That leads to $t=156$ observations for Johnson & Johnson (J&J),
PepsiCo, and Walmart. All quarterly price and dividend data have been
collected from Yahoo Finance.

###  4.1 Dividend Discount Model

Before we move to the numerical results, we provide a brief review of the
dividend discount model (DDM), which is useful for the numerical results. Let
us assume that there are $m$ companies, and the companies will not default in
the future. For DDM with default risk, we refer to Battulga22a. As mentioned
before, the basic idea of all DDMs is that the market price of a stock equals
the sum of the stock’s next period price and dividend discounted at the
required rate of return. Therefore, for successive prices of $i$–th company,
the following relation holds

$P_{i,t}=(1+k_{i,t})P_{i,t-1}-d_{i,t},~{}~{}~{}i=1,\dots,m,~{}t=1,2,\dots,$
(146)

where $k_{i,t}$ is the required rate of return on stock, $P_{i,t}$ is the
stock price, and $d_{i,t}$ is the dividend, respectively, at time $t$ of
$i$–th company. In vector form, the above equation is written by

$P_{t}=(i_{m}+k_{t})\odot P_{t-1}-d_{t},~{}~{}~{}t=1,2,\dots,$ (147)

where $k_{s_{t}}:=(k_{1,t},\dots,k_{m,t})^{\prime}$ is an $(m\times 1)$ vector
of the required rate of returns on stocks at time $t$,
$P_{t}:=(P_{1,t},\dots,P_{m,t})^{\prime}$ is an $(m\times 1)$ price vector at
time $t$, and $d_{t}:=(d_{1,t},\dots,d_{m,t})^{\prime}$ is an $(m\times 1)$
dividend vector at time $t$ of the companies. On the other hand, we assume
that a dividend of $i$–th company is proportional to stock price of the
company. To price dividend paying option, Merton74 used the assumption.
Consequently, dividends of the companies are modeled by

$d_{t}=\alpha_{t}\odot P_{t-1},~{}~{}~{}t=1,2,\dots,$ (148)

where $\alpha_{t}=(\alpha_{1,t},\dots,\alpha_{m,t})^{\prime}$ is an $(m\times
1)$ vector of dividend–to–price ratios at time $t$ of the companies.
Consequently, by equation (147) and (148), the price vector at time $t$ is
written by

$P_{t}=(i_{m}+k_{t}-\alpha_{t})P_{t-1}.$ (149)

To the vector of dividend–to–price ratios take positive values, instead of
equation (148), we work a vector of log dividend–to–price ratios, which is
given by the following equation

$\tilde{d}_{t}=\tilde{\alpha}_{t}-\tilde{P}_{t-1},$ (150)

where $\tilde{d}_{t}:=\ln(d_{t})$ is an $(m\times 1)$ log dividend vector at
time $t$, $\tilde{P}_{t-1}:=\ln(P_{t-1})$ is an $(m\times 1)$ log price vector
at time $t-1$, $\tilde{\alpha}_{t}$ is the $(m\times 1)$ vector of log
dividend–to–price ratios at time $t$. As a result, equation (149) becomes

$P_{t}=\big{(}i_{m}+k_{t}-\exp\\{\tilde{\alpha}_{t}\\}\big{)}P_{t-1}.$ (151)

It follows from equations (147) and (150) that for $t=1,2,\dots$, vectors of
the required rate of returns and log dividend–to–price ratios of the companies
are obtained by

$k_{t}=(P_{t}+d_{t})\oslash P_{t-1}-i_{m}$ (152)

and

$\tilde{\alpha}_{t}=\ln(d_{t}\oslash P_{t-1}),$ (153)

respectively. By repeating equation (151), for $t=0,1,\dots$ and
$r=0,1,\dots,$ we obtain the price vectors at time $t+r$ of the companies

$P_{t+r}=\prod_{j=1}^{r}\big{(}i_{m}+k_{t+j}-\exp\\{\tilde{\alpha}_{t+j}\\}\big{)}\odot
P_{t},$ (154)

where for $q=1,2,\dots$ and generic ($m\times 1$) vectors $o_{1},\dots,o_{q}$,
$\prod_{j=1}^{q}o_{j}=o_{1}\odot\dots\odot o_{q}$ is an element–wise product
of the vectors $o_{1},\dots,o_{q}$ and with convention
$\prod_{j=q}^{q-1}o_{j}=i_{m}$.

###  4.2 Maximum Likelihood Estimation

It is not difficult to show that the required rate of returns of the selected
companies are AR(0) processes, see Figure 1. Thus, for each firm, we model the
required rate of returns by the following AR(0) process

$k_{i,t}=c_{i,s_{t}}+\xi_{i,t},~{}~{}~{}i=1,2,3,$ (155)

where the regime–switching process $s_{t}$ takes one of the values
$\\{1,2,3\\}$ and $c_{i,s_{t}}$ is a constant, corresponding to $i$–th company
and the regime $s_{t}$. We suppose that a variance of the residual process
$\xi_{i,t}$ in the regime $s_{t}$ of the $i$–th company equals
$\sigma_{i,s_{t}}^{2}$.

We present maximum likelihood estimations of the parameters for the selected
companies in Table 1. Maximum likelihood estimation method can be found in
Hamilton90,Hamilton94, see also Battulga22b,Battulga23a. The 2–9th rows of
Table 1 correspond to that the required rate of returns of the companies are
modeled by the regime–switching process with three regimes and the 10–13th
rows of the same Table correspond to that the required rate of returns of the
companies take constant values (the regime–switching process takes one
regime).

Table 1: Maximum Likelihood Estimations of Parameters of Selected Companies

Row | Parameters | Johnson & Johnson | Pepsi | Walmart
---|---|---|---|---
2 | $c_{i,j}$ | 15.64% | 4.26% | –6.41% | 6.29% | 2.73% | –13.90% | 15.09% | 1.74% | –5.60%
3 | $\sigma_{i,j}$ | 0.047 | 0.052 | 0.054 | 0.098 | 0.055 | 0.022 | 0.096 | 0.066 | 0.083
4 | $\hat{\mathsf{P}}_{i}$ | 0.222 | 0.384 | 0.395 | 0.908 | 0.092 | 0.000 | 0.460 | 0.000 | 0.540
5 | 0.099 | 0.645 | 0.256 | 0.000 | 0.933 | 0.067 | 0.000 | 0.932 | 0.068
6 | 0.232 | 0.768 | 0.000 | 0.890 | 0.000 | 0.110 | 0.683 | 0.164 | 0.153
7 | $\tau_{i,j}$ | 1.285 | 2.814 | 1.000 | 10.842 | 14.865 | 1.124 | 1.851 | 14.789 | 1.181
8 | $\pi_{i,j}$ | 0.146 | 0.634 | 0.220 | 0.404 | 0.554 | 0.042 | 0.269 | 0.517 | 0.213
9 | $k_{i,\infty}$ | 3.58% | 3.47% | 3.77%
10 | $c_{i}$ | 3.58% | 3.54% | 4.02%
11 | $k_{i,L}$ | 2.24% | 2.18% | 2.25%
12 | $k_{i,U}$ | 4.92% | 4.90% | 5.80%
13 | $\sigma$ | 0.084 | 0.086 | 0.112

In order to obtain estimations of the parameters, which correspond to the
2–9th rows of Table 1 we assume that the regime–switching process $s_{t}$
follows a Markov chain with three regimes, namely, up regime (regime 1),
normal regime (regime 2), and down regime (regime 3). For the normal regime,
most of the time of the regime–switching process $s_{t}$ lasts in regime 2.
Since explanations are comparable for the other companies, we give
explanations only for J&J. In the 2nd row of Table 1, for each company
$i=1,2,3$, we provide estimations of the parameters $c_{i,1},c_{i,2},c_{i,3}$.
For J&J, in regimes 1, 2, and 3, estimations of the required rate of return
are 15.64%, 4.26%, and –6.41%, respectively. For example, in the normal
regime, the required rate of return of J&J could be 4.26% on average.

In the 3rd row of Table 1, we present parameter estimations of standard
deviations of the error random variables $\xi_{i,t}$ for the selected
companies. For J&J, the parameter estimations of the standard deviations equal
$0.047$ (regime 1), $0.052$ (regime 2), and $0.054$ (regime 3). The 13th row
of Table 1 corresponds to the parameter estimations of standard deviations, in
which the required rate of returns of the companies are modeled by
regime–switching process with one regime. For J&J, the parameter estimation
equals 0.084. As we compare the 9th row and 13th row of the Table, we can see
that the estimations that correspond to the regime–switching process with
three regimes are significantly lower than the ones that correspond to the
regime–switching process with one regime, except the up regime of PepsiCo.

The 4–6th rows of Table 1 correspond to the transition probability matrix
$\hat{\mathsf{P}}$. For the selected companies, their transition probability
matrices $\hat{\mathsf{P}}$s are ergodic, where ergodic means that one of the
eigenvalues of $\hat{\mathsf{P}}$ is unity and that all other eigenvalues of
$\hat{\mathsf{P}}$ are inside the unit circle, see Hamilton94. From the 4th
row of Table 1, one can deduce that if the required rate of return of J&J is
in the up regime, then in the next period, it will switch to the up regime
with a probability of 0.22, the normal regime with a probability of 0.384, or
the down regime with a probability of 0.395. If the required rate of return of
J&J in the normal regime, corresponding to row 5 of the Table, then in the
next period, it will switch to the up regime with a probability of 0.099, the
normal regime with a probability of 0.645, or the down regime with a
probability of 0.256. Finally, if the required rate of return of J&J is in the
down regime, then in the next period, it will switch to the up regime with a
probability of 0.232 or the normal regime with a probability of 0.768 because
of the down regime’s zero probability, see 6th row of the same Table.

We provide the average persistence times of the regimes in the 7th row of
Table 1. The average persistence time of the regime $s_{t}$ is defined by
$\tau_{s_{t}}:=1/(1-p_{s_{t}s_{t}})$ for $s_{t}\in\\{1,2,3\\}.$ From Table 1,
one can conclude that up, normal, and down regimes of J&J’s required rate of
return will persist on average for 1.285, 2.814, and 1.000 quarters,
respectively.

In the 8th row of Table 1, we give ergodic probabilities $\pi$ of the selected
companies. Ergodic probability vector $\pi$ of an ergodic Markov chain is
obtained from an equation $\hat{\mathsf{P}}\pi=\pi$. The ergodic probability
vector represents long–run probabilities, which do not depend on the initial
probability vector $z_{1|0}$. After sufficiently long periods, the required
rate of return of J&J will be in the up regime with a probability of 0.146,
the normal regime with a probability of 0.634, or the down regime with a
probability of 0.220, which are irrelevant to initial regimes.

The 9th row of Table 1 is devoted to long–run expectations of the required
rate of returns of the selected companies. The long–run expectation of the
required rate of return is defined by
$k_{\infty}:=\lim_{t\to\infty}\mathbb{E}(k(s_{t})).$ For J&J, it equals 3.58%.
So that, after long periods, the average required rate of return of J&J
converges to 3.58%.

Finally, for the one regime case, parameter estimations of the required rate
of returns at time Q4, 2023 of the firms are presented in row 10 of the Table,
while the corresponding 95% confidence intervals are included in rows 11 and
12 below. Here we use formulas in Battulga23a to calculate the parameter
estimation and confidence bands. Note that since the required rate of return
estimation expresses the average quarterly return of the companies, we can
convert them yearly using a formula $(1+k)^{4}-1$. The Table further
illustrates average returns (3.58% for J&J) and return variability, as the
return is supposed to lie within the (2.24%, 4.92%) interval with a 95%
probability.

For the selected firms, plotting the smoothed probabilistic inferences with a
return series will be interesting. For each period $u=1,\dots,t$, and each
firm, the smoothed inferences are calculated by equation (123), and the return
series are calculated by equation (152). In Figure 1, we plotted the resulting
series as a function of period $u$. In Figure 1, the left axis corresponds to
the return series, while the right axis corresponds to the smoothed inference
series for each company. From the Figure and the 9th and 13th rows of Table 1,
we can deduce that the regime–switching processes with three regimes are more
suited to explain the required rate of return series as compared to the
regime–switching processes with one regime.

Figure 1: Returns VS Regime Probabilities of Selected Companies

###  4.3 Bayesian Estimation

Now, we aim to obtain future theoretical prices of the selected companies
using equation (154). To get the future prices of the firms, we assume that
the vectors of the required rate of returns and the log dividend–to–price
ratios of the firms are placed on the Bayesian MS–VAR process $y_{t}$. We also
assume that the process $y_{t}$ is modeled by the Bayesian MS–VAR process of
order 4, which is given by the following equation

$y_{t}=a_{0,s_{t}}+A_{1,s_{t}}y_{t-1}+\dots+A_{4,s_{t}}y_{t-4}+\xi_{t},~{}t=1,\dots,T,$
(156)

where $y_{t}=(k_{t}^{\prime},\tilde{\alpha}_{t}^{\prime})^{\prime}$ is a
$(6\times 1)$ vector of endogenous variables with
$k_{t}:=(k_{1,t},k_{2,t},k_{3,t})^{\prime}$ and
$\tilde{\alpha}_{t}:=(\tilde{\alpha}_{1,t},\tilde{\alpha}_{2,t},\tilde{\alpha}_{3,t})^{\prime}$,
$\xi_{t}=(\xi_{1,t},\dots,\xi_{6,t})^{\prime}$ is a $(6\times 1)$ residual
process, $a_{0,s_{t}}$ is a $(6\times 1)$ random coefficient vector at regime
$s_{t}$, and for $i=1,\dots,4$, $A_{i,s_{t}}$ are $(6\times 6)$ random
coefficient matrices at regime $s_{t}$.

For the vector of log dividend–to–price ratios $\tilde{\alpha}_{t}$ of the
selected companies, it can be shown that the ratios follow the unit root
process. On the other hand, the vector of the required rate of returns $k_{t}$
of the selected companies is a stationary process, see Figure 1. Consequently,
conditional expectations of diagonal elements of the random matrix
$A_{1,s_{t}}$ are chosen by $\phi:=(0,0,0,1,1,1)^{\prime}$.

We suppose that the scale matrix $V_{0,s_{t}}$ is a diagonal matrix and its
diagonal elements, corresponding to the process of log dividend–to–price
ratios are estimated by sample variances of the processes, modeled by
univariate AR(4). For other diagonal elements, corresponding to the required
rate of returns, we use the maximum likelihood estimations of parameters,
which are given in Table 1. Therefore, for each regime $s_{t}=1,2,3$, we
choose the scale matrix $V_{0,s_{t}}$ by

$\displaystyle V_{0,1}$ $\displaystyle:=$
$\displaystyle\text{diag}\big{\\{}0.047^{2},0.098^{2},0.096^{2},0.092^{2},0.094^{2},0.133^{2}\big{\\}}$
$\displaystyle V_{0,2}$ $\displaystyle:=$
$\displaystyle\text{diag}\big{\\{}0.052^{2},0.055^{2},0.066^{2},0.092^{2},0.094^{2},0.133^{2}\big{\\}}$
(157) $\displaystyle V_{0,3}$ $\displaystyle:=$
$\displaystyle\text{diag}\big{\\{}0.054^{2},0.022^{2},0.083^{2},0.092^{2},0.094^{2},0.133^{2}\big{\\}}.$

Also, based on the same Table, for each regime, we choose conditional
expectation of the random vector $a_{0,s_{t}}$ by
$c_{1}:=\mathbb{E}\big{[}a_{0,s_{t}}\big{|}\Sigma_{s_{t}},s_{t}=1,\mathcal{F}_{0}\big{]}=(0.15,0.06,0.15,0,0,0)^{\prime}$,
$c_{2}:=\mathbb{E}\big{[}a_{0,s_{t}}\big{|}\Sigma_{s_{t}},s_{t}=2,\mathcal{F}_{0}\big{]}=(0.04,0.03,0.02,0,0,0)^{\prime}$,
and
$c_{3}:=\mathbb{E}\big{[}a_{0,s_{t}}\big{|}\Sigma_{s_{t}},s_{t}=3,\mathcal{F}_{0}\big{]}=(-0.07,-0.14,-0.05,0,0,0)^{\prime}$.
For the other hyperparameters, we choose $\varepsilon_{j}=\lambda_{1,j}=20$
and $\lambda_{2,i}=\tau_{i,j}=1$ for each regime $j=1,2,3$ and each endogenous
variable $i=1,\dots,6$. As a result, the matrices of the dummy variables
$\hat{y}_{s_{t}}^{\circ}$ and $\hat{\mathsf{Y}}_{s_{t}}^{\circ}$, which are
given in equations (113) and (114) are represented by

$\hat{y}_{s_{t}}^{\circ}:=20\big{[}c_{s_{t}}:\text{diag}\\{0,0,0,1,1,1\\}:0_{[6\times
20]}\big{]}$ (158)

and

$\hat{\mathsf{Y}}_{s_{t}}^{\circ}:=20\begin{bmatrix}1&0_{[1\times 24]}\\\
0_{[24\times
1]}&\big{(}\text{diag}\\{1,2,3,4\\}\otimes\text{diag}\\{1,1,1,1,1,1\\}\big{)}\end{bmatrix},$
(159)

where dimensions of the matrices of the dummy variables
$\hat{y}_{s_{t}}^{\circ}$ and $\hat{\mathsf{Y}}_{s_{t}}^{\circ}$ are $(6\times
25)$ and $(25\times 25)$, respectively.

In order to obtain future eight quarters’ theoretical prices and their 95%
confidence bands of the selected companies, we use equation (154) and make a
simulation of $\mathcal{L}=10,000$. The following Table presents the results
of the theoretical prices and confidence bands.

Table 2: Future Theoretical Prices and 95% Confidence Bands of Selected
Companies

Companies | $P_{t}$ | $P_{t+1}$ | $P_{t+2}$ | $P_{t+3}$ | $P_{t+4}$ | $P_{t+5}$ | $P_{t+6}$ | $P_{t+7}$ | $P_{t+8}$
---|---|---|---|---|---|---|---|---|---
J&J | T/Price | 158.90 | 163.44 | 167.00 | 170.74 | 174.29 | 178.03 | 181.79 | 185.65 | 189.84
L/Bound | 158.90 | 143.04 | 137.33 | 133.97 | 131.30 | 129.46 | 127.21 | 125.87 | 124.01
U/Bound | 158.90 | 184.10 | 199.14 | 213.59 | 225.62 | 238.55 | 250.06 | 261.65 | 273.48
PepsiCo | T/Price | 168.53 | 171.08 | 173.25 | 175.43 | 177.33 | 179.30 | 181.38 | 183.39 | 185.45
L/Bound | 168.53 | 150.54 | 144.19 | 140.09 | 136.47 | 133.36 | 130.78 | 128.37 | 125.59
U/Bound | 168.53 | 191.43 | 203.46 | 213.40 | 222.56 | 231.38 | 239.09 | 248.31 | 255.65
Walmart | T/Price | 55.08 | 57.00 | 58.75 | 60.70 | 62.48 | 64.47 | 66.40 | 68.40 | 67.38
L/Bound | 55.08 | 47.94 | 45.42 | 44.31 | 43.28 | 42.58 | 41.76 | 41.05 | 40.62
U/Bound | 55.08 | 66.19 | 72.91 | 79.20 | 84.91 | 91.24 | 97.23 | 102.37 | 109.35

For Table 2, T/Price, L/Bound, and U/Bound mean theoretical price, lower
bound, and upper bound, respectively. The column $P_{t}$ represents real
prices at time Q4, 2023 of the selected companies. For each $u=1,\dots,8$, the
column $P_{t+u}$ expresses future prices at time $t+u$ (after $u$ quarters
from Q4, 2023) of the companies. In particular, the column $P_{t+4}$
represents future prices in Q4, 2024. The same column further explains that
the theoretical price in Q4, 2024 of J&J equals 174.29, and the price in Q4,
2024 of the firm will lie within the (131.30,225.62) interval with a 95%
probability. The same explanations hold for other columns and other companies.
From the Table, one can deduce that histograms of the future prices are right
skewed, and the skewness increases as time increases.

##  5 Conclusion

In this paper, for the general Bayesian MS–VAR process, we obtain some useful
density functions for Monte–Carlo simulations. The density functions tell us
that conditional on the regime vectors and initial information, the vector of
endogenous variables is independent of model’s some random components. Thus,
one only needs the prior distributions to calculate the density functions, and
the results have yet to be explored before.

In a particular case of the general Bayesian MS–VAR process, we also get
closed–form density functions of the random components of the model. In
particular, we find that joint distributions of future values of the random
coefficient matrix are matrix variate student distribution, see equation (2).
Hence, one can analyze impulse response by directly generating the coefficient
matrix from the density. Also, we provide a new density function, which has
yet to be introduced before of future values of the endogenous variables; see
equations (44) and (77). Thus, future studies may concentrate on marginal
density functions and direct simulation methods for the density function.
Further, we obtain a characteristic function of the random coefficient matrix,
which can be used to calculate the forecast of the endogenous variables.

In the paper, we developed Monte–Carlo simulation algorithms. The simulation
method’s novelty is that it removes the regime vector duplication. As a
result, our proposed Monte–Carlo simulation method departs from the previous
simulation methods with regime switching. We also provide importance sampling
method to estimate probability of a rare event, corresponding to the future
endogenous variables. Since the method can be used to calculate quantiles, in
this case, the quantiles of the future endogenous variables become more
reliable than navy simulation methods.

## References