Bayesian Markov-Switching Vector Autoregressive Process

Bayesian Markov--Switching Vector Autoregressive Process

Keywords: Bayesian MS--VAR process, Monte--Carlo simulation methods.

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 \citeA{Sims80}, there are some other important works on multiple time series modeling, see, e.g., \citeA{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 \citeA{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, \citeA{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 \citeA{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 \citeA{Hamilton89,Hamilton90} (see also books of \citeA{Hamilton94} and \citeA{Krolzig97}), and the model is hidden Markov model with dependencies; see \citeA{Zucchini16}. However, Markov regime--switching models have been introduced before Hamilton (1989), see, \citeA{Goldfeld73}, \citeA{Quandt58}, and \citeA{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 \citeA{Battulga23a}, \citeA{Battulga24f}, and \citeA{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, \citeA{Battulga24a} derived pricing formulas. Also, the author used the Bayesian MS--VAR process to price equity--linked life insurance products and rainbow options, see \citeA{Battulga23a} and \citeA{Battulga24f} .

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 \citeA{Albert93}, and its multidimensional extension is given by \citeA{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. \citeA{Glasserman00} used the importance sampling method to model portfolio loss random variable by using approximation. Also, \citeA{Glasserman05b} study a loss random variable of credit portfolio applying the method, see also \citeA{McNeil05}.

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. Finally, Section 4 concludes the study.

Bayesian MS--VAR 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

where $y_t=(y_{1,t},\dots,y_{n,t})'$ is an $(n\times 1)$ vector of endogenous variables, $\psi_t=(1,\psi_{2,t},\dots,\psi_{l,t})'$ is an $(l\times 1)$ vector of exogenous variables, $\xi_t=(\xi_{1,t},\dots,\xi_{n,t})'$ 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 (ref) can be written by

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',y_{t-1}',\dots,y_{t-p}')'$ 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 \citeA{Lutkepohl05} and \citeA{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

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. Here we assume that initial values of the random covariance matrix are $\Sigma_{s_{1-j}}=\Sigma_{1-j}$ for $j=1,\dots,q_*$.

Let us introduce stacked vectors and matrices: $y:=(y_1',\dots,y_T')'$, $s:=(s_1,\dots,s_T)'$, $\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}',\dots,y_0',\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$ are the initial values of the random matrix process $\Sigma_{s_t}$, $\psi_1,\dots,\psi_T$ are the 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',\dots,o_T')'$, 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',\dots,o_t')'$ and $\bar{o}_t^c:=(o_{t+1}',\dots,o_T')'$. 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,

for $t=0,\dots,T-1$, where $p_{s_1}:=p_{0s_1}=\mathbb{P}[s_1=s_1|\mathsf{P},\mathcal{F}_0]$ is the initial probability. A distribution of a residual random vector $\xi:=(\xi_1',\dots,\xi_T')'$ is given by

where $\Sigma_s:=\text{diag}\{\Sigma_{s_1},\dots,\Sigma_{s_T}\}$ is an $([nT]\times [nT])$ 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)'$, we define sets

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

and

and intersection sets of the sets of right hand sides of equations (ref) and (ref), and (ref) are empty sets. We introduce the following regime vectors: $\hat{s}:=(\hat{s}_1,\dots,\hat{s}_{r_{\hat{s}}})'$ is an $(r_{\hat{s}}\times 1)$ vector, $\alpha:=(\alpha_1,\dots,\alpha_{r_\alpha})'$ is an $(r_\alpha\times 1)$ vector, $\beta=(\beta_1,\dots,\beta_{r_\beta})'$ is an $(r_\beta\times 1)$ vector, $\gamma=(\gamma_1,\dots,\gamma_{r_\gamma})'$ is an $(r_\gamma\times 1)$ vector, $\delta=(\delta_1,\dots,\delta_{r_\delta})'$ is an $(r_\delta\times 1)$ vector, and $\epsilon=(\epsilon_1,\dots,\epsilon_{r_\epsilon})'$ is an $(r_\epsilon\times 1)$ vector. For the regime vector $a=(a_1,\dots,a_{r_a})' \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, conditional on $\hat{s}$ and $\mathcal{F}_0$, a joint density function of the random coefficient matrix $(\Pi_{\hat{s}},\Gamma_{\hat{s}})$ is represented by

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$. Using the regime vectors $\alpha$ and $\delta$, the above joint density function can be written by

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

Then, the following Proposition, which is useful for Monte--Carlo simulation holds, see below. The proofs of this one and other Propositions are given in Appendix.

It follows from the Proposition that (i) conditional on $\bar{s}_t$, $\mathsf{P}$, and $\mathcal{F}_0$, $\bar{s}_t^c$ and $(\bar{y}_t,\Pi_\alpha,\Gamma_\alpha)$ are independent, (ii) conditional on $\delta$ and $\mathcal{F}_0$, $(\Pi_{\delta},\Gamma_{\delta})$ and $(\bar{y}_t,\Pi_{\alpha},\Gamma_{\alpha},\alpha,\mathsf{P})$ are independent, (iii) conditional on $\bar{s}_t$ and $\mathcal{F}_0$, $\mathsf{P}$ and $(\bar{y}_t,\Pi_\alpha,\Gamma_\alpha)$ are independent, and (iv) conditional on $\bar{s}_t$ and $\mathcal{F}_t$, $(\Pi_\alpha,\Gamma_\alpha)$ and $\mathsf{P}$ are independent. It should be noted that according to the Markov property (ref), it follows from equation (ref) that the assumption for a Markov chain in the book of \citeA{Hamilton94} always holds, namely,

Special Case of Bayesian MS--VAR 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

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

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

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

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$.

Distributions

For the regime vector $\bar{s}_t$ and regime $\alpha_k$, we define sets

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

Further, we define indexes

The index $o_t$ represents a position of the regime $s_t$ in the regime vector $\alpha$. Let $\pi_{\alpha}:=\big(\pi_{\alpha_1}',\dots,\pi_{\alpha_{r_{\alpha}}}'\big)'$ be an $([ndr_{\alpha}]\times 1)$ duplication removed 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}}',\dots,y_{k_{t,q_{t,\alpha_k}}}'\Big)'$ 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}]'$, 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.

Conditional Densities

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

and

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]'$ is an $([nr_{\alpha}]\times n)$ matrix, $\pi_{0,\alpha}:=\big(\pi_{0,\alpha_1}',\dots,\pi_{0,\alpha_{r_{\alpha}}}'\big)'$ 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

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 \citeA{Battulga24a}, a distribution of the $([nT]\times 1)$ random vector $y=(y_1',\dots,y_T')'$ is given by

where the matrix $\Psi_s$ and the vector $\varphi_s$ are

and

respectively. To price default--free options, \citeA{Battulga24a} used the conditional distribution of the random vector $y$. For a generic vector $o=(o_1',\dots,o_n')'$ 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.

It follows from equations (ref) and (ref) that sub coefficient vectors and sub covariance matrices are conditional independent. Note that the conditional independence is consistent with the assumption (ref). From equations (ref) and (ref) 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

and

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 (ref) 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

and

Because, according to equation (ref), density function (ref) 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 (ref) and (ref) 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

Consequently, due to equation (ref), the above density function is represented by a product of the conditional matrix variate student $t$ density functions. Note that if $r_\delta=0$, we must eliminate the second line of the above equation.

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

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})'$, and Dirichlet distribution by $\text{Dir}(\alpha_i)$. Then, the following Lemma holds.

Consequently, it is worth mentioning that according to equation (ref) in the above Proposition, conditional on $\mathcal{F}_0$, the regime--switching process $s_t$ is not a Markov chain because the conditional density function depends on the regime vector $\bar{s}_{t-1}$.

Characteristic Function

For $k=1,\dots,r_\gamma$, equation (ref) can be written by

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

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

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

see \citeA{Butler98}. Its one dimensional version is called generalized inverse Gaussian distribution and it is widely used to model returns of financial assets, see \citeA{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

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 (ref), 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

Similarly, it can be shown that

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

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

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 (ref) and (ref) 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 (ref). 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,

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

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

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

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 \citeA{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

\citeA{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

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

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

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 \citeA{Litterman79}. Also, \citeA{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 \citeA{Miranda18}. The idea of Minnesota prior is that it shrinks diagonal elements of the matrix $A_{1,s_t}$ 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,s_t},A_{1,s_t},\dots,A_{p,s_t}$ are jointly normally distributed, and for $(i,j)$--th element of the matrix $A_{\ell,s_t}$ $(\ell=0,\dots,p)$, it holds

and

where $\sigma_{i,s_t}^2$ is an $(i,i)$--th element of the random covariance matrix $\Sigma_{s_t}$. 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 \citeA{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 \citeA{Banbura10}, it can be shown that the following equation satisfies the prior conditions (ref), (ref), and (ref)

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

and

with $J_{s_t}:=\text{diag}\{1^{\lambda_{2,s_t}},\dots,p^{\lambda_{2,s_t}}\}$, respectively, and $\hat{\xi}_{s_t}^\circ:=[\xi_{1,s_t}:\dots:\xi_{d,s_t}]$ is an $(n\times d)$ matrix, whose columns are conditional independent for the given covariance matrix $\Sigma_{s_t}$, regime $s_t$, and initial information $\mathcal{F}_0$ and for $i=1,\dots,d$, each column has a distribution $\xi_{i,s_t}|\Sigma_{s_t},s_t,\mathcal{F}_0\sim \mathcal{N}(0,\Sigma_{s_t})$. Note that one can add constraints for elements of the coefficient matrix $\Pi_{s_t}$ to the matrices of dummy variables. 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 (ref). For this reason, we choose the prior condition (ref), c.f. \citeA{Banbura10}. Equation (ref), can be written by

where $\hat{y}_{s_t}$ and $\hat{\xi}_{s_t}$ 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}_{s_t}^\circ$, respectively, i.e., $\hat{y}_{s_t}:=\text{vec}(\hat{y}_{s_t}^\circ)$ and $\hat{\xi}_{s_t}:=\text{vec}(\hat{\xi}_{s_t}^\circ)$. It follows from equation (ref) that

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)'(\hat{\mathsf{Y}}_{s_t}^\circ(\hat{\mathsf{Y}}_{s_t}^\circ)')^{-1}\big)$. Consequently, since $\hat{\xi}_{s_t}|\Sigma_{s_t},s_t,\mathcal{F}_0\sim \mathcal{N}\big(0,I_d\otimes \Sigma_{s_t}\big)$, conditional on $\Sigma_{s_t}$, $s_t$, and $\mathcal{F}_0$, a distribution of the random coefficient vector $\pi_{s_t}$ is given by

where $\Lambda_{0,s_t}:=(\hat{\mathsf{Y}}_{s_t}^\circ(\hat{\mathsf{Y}}_{s_t}^\circ)')^{-1}$ is a $(d\times d)$ diagonal matrix. Therefore, 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 (ref), (ref), and (ref), we can use the results of Proposition (ref).

Simulation Methods

By applying Propositions (ref) and (ref), one can obtain exact density function of the random vector $\bar{y}_t$, namely,

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

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

General Method

According to the conditional probability formula and Proposition (ref), we have that

The above equation tells us that how to generate random samples $(\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 the random samples, firstly, we generate the random coefficient vector, random covariance matrix, regime vector, and transition probability matrix $(\pi_\alpha,\Sigma_\alpha,\bar{s}_t,\mathsf{P})$ from the posterior density function $f(\pi_\alpha,\Sigma_\alpha,\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|s_t,\mathsf{P},\mathcal{F}_0)$, so on.

First, we consider a simulation method that generate the random coefficient vector, random covariance matrix, regime vector, and transition probability matrix $(\pi_\alpha,\Sigma_\alpha,\bar{s}_t,\mathsf{P})$ from the posterior density function $f(\pi_\alpha,\Sigma_\alpha,\bar{s}_t,\mathsf{P}|\mathcal{F}_t)$. We develop the Gibbs sampling method to generate $(\pi_\alpha,\Sigma_\alpha,\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. According to Proposition (ref), constructing the Gibbs sampler to approximate the joint posterior distribution $f(\bar{s}_t,\mathsf{P}|\mathcal{F}_t)$ is straightforward: set initial values $\big(\pi_{\alpha(0)}(0),\Sigma_{\alpha(0)}(0),\bar{s}_t(0),\mathsf{P}(0)\big)$ and new values $\big(\pi_{\alpha(\ell)}(\ell),\Sigma_{\alpha(\ell)}(\ell),\bar{s}_t(\ell),\mathsf{P}(\ell)\big)$, $\ell=1,\dots,\mathcal{L}$ can be generated by

where $\alpha(\ell+1)$ is the duplication removed regime vector, corresponding to the regime vector $\bar{s}_t(\ell+1)$.

Now, we consider a sampling method that generate $\bar{s}_t(\ell+1)$ from $f\big(\bar{s}_t\big|\lambda_t(\ell),\mathcal{F}_t\big)$, where the parameter vector at iteration $\ell$, $\lambda_t(\ell):=\big(\pi_{\alpha(\ell)}(\ell),\Sigma_{\alpha(\ell)}(\ell),\mathsf{P}(\ell)\big)$ corresponds to the regime vector $\bar{s}_t$. Here we follow some results of the book of \citeA{Hamilton94}, see also \citeA{Battulga23b} and \citeA{Battulga22b}. If we assume that the regime--switching process in regime $j$ at time $u$, then the conditional density function of the random vector $y_u$ is given by the following equation

for $u=1,\dots,t$ and $j=1,\dots,N$. 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(\ell):=(\eta_{u,1}(\ell),\dots,\eta_{u,N}(\ell))'$. Let us denote a probabilistic inference about the value of the regime--switching process $s_u$ is equal to $j$, based on the random coefficient vector $\pi_{\alpha(\ell)}(\ell)$, random covariance matrix $\Sigma_{\alpha(\ell)}(\ell)$, information $\mathcal{F}_u$, and transition probability matrix $\mathsf{P}(\ell)$ by $\mathbb{P}(s_u=j|\lambda_t(\ell),\mathcal{F}_u)$. Collect these conditional probabilities $\mathbb{P}(s_u=j|\lambda_t(\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|\lambda_t(\ell),\mathcal{F}_u),\dots,\mathbb{P}(s_u=N|\lambda_t(\ell),\mathcal{F}_u)\big)'$. Also, we need a probabilistic forecast about the value of the regime--switching process at time $u+1$ is equal to $j$ conditional on the random coefficient vector $\pi_{\alpha(\ell)}(\ell)$, random covariance matrix $\Sigma_{\alpha(\ell)}(\ell)$, transition probability matrix $\mathsf{P}(\ell)$, and data up to and including time $u$, $\mathcal{F}_u$. 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|\lambda_t(\ell),\mathcal{F}_u),\dots,\mathbb{P}(s_{u+1}=N|\lambda_t(\ell),\mathcal{F}_u)\big)'$.

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

where $\odot$ is the Hadamard product of two vectors, $\eta_u(\ell)$ is the $(N\times 1)$ vector, whose $j$--th element is given by equation (ref), $\hat{\mathsf{P}}(\ell)$ is the $(N\times N)$ transition probability matrix, which omits the first row of the transition probability matrix $\mathsf{P}(\ell)$, 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)'$ one can iterate on (ref) 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|\lambda_t(\ell),\mathcal{F}_t),\dots,\mathbb{P}(s_u=N|\lambda_t(\ell),\mathcal{F}_t)\big)'$ for $u=1,\dots,t$. In practice, a popular method to calculate smoothed probability inference is the \citeA{Kim94}'s smoothing algorithm, which is based on approximation. In this paper, we introduce a new and simple smoothing method, which is also based on an approximation and is given the following Proposition.

The smoothed probabilities $z_{u|t}(\ell)$ are found by iterating on (ref) and (ref) backward for $u=t-1,\dots,1$. This iteration uses probabilistic inferences and forecasts $z_{u|u}(\ell)$ and $z_{u+1|u}(\ell)$ for $u=1,\dots,t-1$, which are obtained from (ref). Thus, the simulation method that generate $\bar{s}_t(\ell)$ for given $\mathsf{P}(\ell)$ and $\mathcal{F}_t$ as follows:

Collect $s_u(\ell+1)$ for $u=1,\dots,t$ into $(t\times 1)$ vector $\bar{s}_t(\ell+1)$, namely, $\bar{s}_t(\ell+1):=(s_1(\ell+1),\dots,s_t(\ell+1))'$. The smoothing method is not only used to Bayesian estimation but also maximum likelihood estimation (ML) of parameters of a model with regime switch. Also, to estimate parameters of the model, one can use the joint density function. It is worth mentioning that if vector of the endogenous variables $y_t$ does not depend $(\pi_{\alpha},\Sigma_\alpha)$, then the smoothing method and joint density function become exact. Therefore, when they are used to the ML estimation of the model's parameters, the ML estimation becomes exact.

Let us consider a simulation method that generate the random coefficient vector $\pi_{\alpha(\ell+1)}(\ell+1)$ and random covariance matrix $\Sigma_{\alpha(\ell+1)}(\ell+1)$ from the density function $f\big(\pi_{\alpha(\ell+1)},\Sigma_{\alpha(\ell+1)}\big|\alpha(\ell+1),\mathcal{F}_t\big)$. Let $\mathcal{A}_{{\bar{s}}_t(\ell+1)}=\big\{\alpha_1(\ell+1),\dots,\alpha_{r_{\alpha(\ell+1)}}(\ell+1)\big\}$ be the duplication removed set, corresponding the regime vector $\bar{s}_t(\ell+1)$. Then, according to equations (ref), (ref) and (ref), a simulation method that generates $\big(\pi_{\gamma_k(\ell+1)}(\ell),\Sigma_{\gamma_k(\ell+1)}(\ell+1)\big)$ as follows: for $k=1,\dots,r_{\alpha(\ell+1)}$,

where $\nu_{0,\alpha_k(\ell+1)|t}:=\nu_{0,\alpha_k(\ell+1)}+q_{t,\alpha_k(\ell+1)}$ and $q_{t,\alpha_k(\ell+1)}$ is a number of regimes in the regime vector $\bar{s}_t(\ell+1)$ that equal the regime $\alpha_k(\ell+1)$.

To generate $\mathsf{P}(\ell+1)$ from $f(\mathsf{P}|\bar{s}_t(\ell+1),\mathcal{F}_0)$, we apply equation (ref) in Proposition (ref), density function (ref), equation (ref) in Proposition (ref), and the Bayesian formula. Then, we have that

Thus, one can deduce that conditional on $\bar{s}_t(\ell+1)$ and $\mathcal{F}_0$, 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+1)):=(\alpha_{i1}+n_{i1}(\bar{s}_t(\ell+1)),\dots,\alpha_{iN}+n_{iN}(\bar{s}_t(\ell+1)))'$. Consequently, a simulation method that generate $\mathsf{P}(\ell+1)$ for given $\bar{s}_t(\ell+1)$ and $\mathcal{F}_0$ as follows:

Collect $\mathsf{P}_i(\ell+1)$ for $i=0,\dots,N$ into an $([N+1]\times N)$ matrix $\mathsf{P}(\ell+1)$, that is, $\mathsf{P}(\ell+1):=[\mathsf{P}_0(\ell+1)':\dots:\mathsf{P}_N(\ell+1)']'.$

Second, for $\ell=1,\dots,\mathcal{L}$, we consider a simulation method that generates 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

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:

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))'$.

Third, for $\ell=1,\dots,\mathcal{L}$, we consider a simulation method that generates the coefficient vector $\pi_{\delta(\ell)}(\ell)$ and covariance matrix $\Sigma_{\delta(\ell)}(\ell)$ from the density function $f\big(\pi_{\delta(\ell)},\Sigma_{\delta(\ell)}\big|\delta(\ell),\mathcal{F}_0\big)$. Let $\mathcal{A}_{{\bar{s}}_t^c(\ell)}=\big\{\beta_1(\ell),\dots,\beta_{r_{\beta(\ell)}}(\ell)\big\}$ be the duplication removed set, corresponding the regime vector $\bar{s}_t^c(\ell)$. To eliminate unnecessary simulations, instead of the regimes in the set $\mathcal{A}_{{\bar{s}}_t^c(\ell)}$, one should consider regimes in difference set of the sets $\mathcal{A}_{{\bar{s}}_t(\ell)}$ and $\mathcal{A}_{{\bar{s}}_t^c(\ell)}$. Let us assume that the difference set of the sets is given by $\mathcal{A}_{{\bar{s}}_t^c(\ell)}\backslash \mathcal{A}_{{\bar{s}}_t(\ell)}=\big\{\delta_1(\ell),\dots,\delta_k(\ell)\big\}$. Then, according to equations (ref), (ref), (ref), and (ref), a simulation method that generates $\big(\pi_{\delta_k(\ell)}(\ell),\Sigma_{\delta_k(\ell)}(\ell)\big)$ as follows: if $r_{\delta(\ell)}>0$, then for $k=1,\dots,r_{\delta(\ell)}$,

Note that if $r_{\delta(\ell)}=0$, we do not need to generate the coefficient vector $\pi_{\delta(\ell)}(\ell)$ and covariance matrix $\Sigma_{\delta(\ell)}(\ell)$ from the density function $f\big(\pi_{\delta(\ell)},\Sigma_{\delta(\ell)}\big|\delta(\ell),\mathcal{F}_0\big)$. Let

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_{\alpha(\ell)}$, $\Sigma_{\alpha_k(\ell)}(\ell)$ and $\pi_{\alpha_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 $([nr_{\beta(\ell)}]\times n)$ matrix $\Sigma_{\hat{\beta}(\ell)}(\ell)$, namely,

and

Similar to equation (ref), 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]'$, 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

Fourth, for $\ell=1,\dots,\mathcal{L}$, we consider a simulation method that generates the regime vector $\bar{y}_t^c(\ell)$ from the density function $f(\bar{y}_t^c|\pi_{\bar{s}_t^c(\ell)}(\ell),\Sigma_{\bar{s}_t^c(\ell)}(\ell),\bar{s}_t^c(\ell),\mathcal{F}_t)$. Let us assume that

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 \citeA{Battulga24a}, a distribution of the random vector $\bar{y}_t^c$ is given by

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:

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

and

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 \citeA{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.

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

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

if $s_u\in\mathcal{A}_{\bar{s}_t}\cap\mathcal{A}_{\bar{s}_t^c}$, then

if $s_u\in\mathcal{A}_{\bar{s}_t^c}\backslash\mathcal{A}_{\bar{s}_t}$, then

and

for $u=t+1,\dots,T$, where $A_{s_u}:=(\Lambda_{0,s_u}\otimes \Sigma_{s_u})$ is an $([nd]\times [nd])$ matrix, $\theta_{s_u}$ is a positive constant, depending on the 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

where the random variable $X_u$ is given by $X_u:=z_u'y_u$ and the quadratic function $\psi(\theta_{s_u})$ for $\theta_{s_u}$ is given by

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 (ref), the following inequality holds

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 \citeA{Glasserman00}. Also, for the second order moment, it holds

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

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}$,

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 distribution of future values of the random coefficient matrix is a product of matrix variate student distributions, see equation (ref). Hence, one can analyze impulse response by directly generating the coefficient matrices from the distribution functions. Also, we provide a new density function, which has yet to be introduced before of future values of the endogenous variables; see equations (ref) and (ref). 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 develop Monte--Carlo simulation algorithms. The simulation method's novelty is that it removes duplication in a regime vector. 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. To obtain smoothed probability inference, \citeA{Kim94}'s smoothing method is widely used in practice. But this method is based on approximation. For this reason, we introduce a very simple, exact (in some special cases), and new smoothing method.