Forecasting with a Bivariate Hysteretic Time Series Model Incorporating Asymmetric Volatility and Dynamic Correlations

1. Introduction

Shocks to a time series can cause persistent effects, whereby the influence of disturbances spreads and persists over time. This phenomenon, referred to as the hysteresis effect, reflects a form of path dependence in which system dynamics respond asymmetrically to past shocks. To address issues related to excessive or spurious regime shifts, a range of univariate hysteretic time series models have been developed by authors of ([6,7,19,20,24,25,26]).

In financial econometrics, it is well-established that asset returns tend to exhibit co-movement. Understanding and forecasting the temporal dependence in the second-order moments of these returns is a key concern in finance. Multivariate models provide a useful framework for capturing complex features such as volatility clustering across multiple assets, time-varying correlations, and joint downside tail risks across industries. These considerations have led researchers to extend univariate volatility models into the multivariate setting. For instance, authors of [5] introduced the VECH and BEKK models, while authors of [15] proposed the Dynamic Conditional Correlation (DCC) model, which allows for a time-varying conditional correlation matrix. In contrast, authors of [22] developed a model that captures correlation dynamics through a weighted average of past correlation matrices, reflecting the persistence of conditional correlations. Authors of [3] develop an asymmetric Dynamic Conditional Correlation (AG-DCC) model to examine the presence of asymmetric responses in conditional volatility and correlation between financial asset returns, particularly allowing for asymmetries in the correlations. A comprehensive discussion on generalized univariate volatility models can be found in [23]. Authors of [10] suggest an extension of [22] using a Bayesian Markov chain Monte Carlo (MCMC) technique to accommodate heavy-tailed distributions. Nonetheless, these models do not consider regime-switching behavior, which is potentially essential for modeling structural shifts and regime-dependent dynamics in financial markets.

In the multivariate context, authors of [8] propose the Hysteretic Vector Autoregressive (HVAR) model, which incorporates delayed regime switching based on a hysteresis variable. Specifically, transitions between regimes occur only when this variable exits a predefined hysteresis zone. Aauthors of [9] introduce a bivariate HAR model incorporating GARCH errors and time-varying correlations. This model integrates features of dynamic correlation, asymmetric effects on correlation and volatility, and heavy-tailed distribution within the multivariate HAR framework previously developed by [8]. However, the asymmetry in [9] is introduced only through the regime-switching behavior of a hysteresis variable within the system. This approach overlooks the leverage effects associated with individual asset returns, which have been emphasized in earlier studies in [3]. In univariate framework, authors of [21] and [17] examine the intricate dynamics between financial returns and volatility, emphasizing the asymmetric effects of shocks. Authors of [17] modify the GARCH model to account for seasonal volatility patterns, differential impacts of positive and negative return innovations, and the influence of nominal interest rates on conditional variance. Similarly, of [21] generalizes the ARCH framework by modeling conditional variance as a quadratic function of past innovations, allowing for a nuanced capture of volatility patterns, including asymmetries and leverage effects. Both studies underscore the importance of accommodating asymmetries in volatility modeling to better understand and predict financial market behaviors.

As a result, the volatility specification in [9] leaves room for improvement in modeling asymmetric effects at the level of individual return series. In this paper, we develop an extension of the multivariate hysteretic autoregressive (MHAR) model with GARCH errors and dynamic correlations (see [9]) to accommodate asymmetries in volatility dynamics. Specifically, we incorporate two well-known asymmetric volatility specifications: the GJR-GARCH, as defined in [17], and the QGARCH proposed in [21]. These extensions result in two model variants, namely the MHAR–GJR–GARCH and the MHAR–QGARCH models.

To the best of our knowledge, this is the first study to explore asymmetric volatility structures within the MHAR–GARCH framework. By introducing these asymmetric components, the proposed models offer greater flexibility in capturing the heterogeneity and nonlinear behavior commonly observed in financial asset returns. Such flexibility is particularly important in modeling risk dynamics, especially during periods of market turbulence where asymmetries in volatility play a crucial role. Based on the proposed models, we employ an adaptive multivariate t-distribution to account for heavy-tailed errors, and utilize the SMN representation (see [2]) to flexibly model marginal error distributions with varying degrees of freedom, improving the model’s fit to the target time series.

In finance, systemic risk refers to the possibility that problems in one financial institution or a group of them could spread throughout the financial system due to the strong connections between institutions. Such a chain reaction can lead to serious disruptions or even cause the entire market to collapse. Following [1], the Marginal Expected Shortfall (MES) is employed to empirically evaluate the extent to which this risk measure addresses practical concerns related to systemic risk, using a large sample of major U.S. banks. In this study, we consider two widely used risk measures: Value at Risk (VaR) and MES, which MES plays a more prominent role in capturing tail risk and systemic vulnerability. Additionally, we implement two backtesting procedures to assess the accuracy of out-of-sample VaR forecasts.

A major limitation of the proposed models lies in their increasing complexity, particularly due to the large number of parameters that must be estimated and the challenges involved in modeling nonlinear multivariate structures. As the nonlinearity and asymmetric structures of the proposed models, traditional estimation methods become inefficient or impractica. To overcome these difficulties, we adopt a Bayesian framework using Markov Chain Monte Carlo (MCMC) techniques, which allows for simultaneous inference of all unknown parameters.

The remainder of this paper is divided into the following sections: The multivariate hysteretic autoregressive model with time-varying correlations and asymmetry structures in volatility is presented in Section 2. Bayesian inference for model parameters is presented in Section 3. Forecasting VaR and the marginal expected shortfall(MES) are mentioned in Section 4. Section 5 examines simulation. The empirical study is demonstrated in Section 6 and marks are covered in Section 7.

2. Multivariate Hysteretic Autoregressive Model with Asymmetry Structures in Volatility and Time-Varying Correlation

Consider the MHAR-GARCH model, which is a multivariate hysteretic autoregressive model with GARCH errors:

$$\begin{array}{ccc}\hfill {\mathit{y}}_t& =& {\Phi }_0^{\left(J_t\right)}+\sum _{l=1}^c{\Phi }_l^{\left(J_t\right)}{\mathit{y}}_{t-l}+{\mathit{a}}_t,\hfill \\ \hfill {\mathit{a}}_t& =& \mathrm{diag}\left(\sqrt{h_{1,t}},\dots ,\sqrt{h_{k,t}}\right){\mathit{ϵ}}_{t}where{\mathit{ϵ}}_t\sim D(\mathbf{0},{\Sigma }_t),\hfill \\ \hfill h_{i,t}& =& {\omega }_{i0}^{\left(J_t\right)}+\sum _{l=1}^{p_i}{\alpha }_{il}^{\left(J_t\right)}{a}_{i,t-l}^2+\sum _{l=1}^{q_i}{\beta }_{il}^{\left(J_t\right)}{h}_{i,t-l},i=1,\dots ,k,\hfill \\ \hfill {\Sigma }_t& =& \left(1-{\theta }_1^{\left(J_t\right)}-{\theta }_2^{\left(J_t\right)}\right){\Sigma }^{\left(J_t\right)}+{\theta }_1^{\left(J_t\right)}{\Sigma }_{t-1}+{\theta }_2^{\left(J_t\right)}{\Psi }_{t-1},\hfill \end{array}$$

(1)

$$\begin{array}{ccc}\hfill J_t& =& \left\{\begin{array}{cc}1\hfill & \mathrm{if}{z}_t\le r_L\hfill \\ 2\hfill & \mathrm{if}{z}_t>r_U\hfill \\ J_{t-1}\hfill & \mathrm{otherwise},r_L<r_U,\hfill \end{array}\right.\hfill \end{array}$$

(2)

and the $(u,v)$th element of ${\mathbf{\Psi }}_{t-1}$ is formulated as:

$${\psi }_{uv,t-1}={\displaystyle \frac{{\displaystyle \sum _{s=1}^S}{ϵ}_{u,t-s}{ϵ}_{v,t-s}}{\sqrt{\left({\displaystyle \sum _{s=1}^S}{ϵ}_{u,t-s}^2\right)\left({\displaystyle \sum _{s=1}^S}{ϵ}_{v,t-s}^2\right)}}},1\le u<v\le S^{★}\le S.$$

(3)

where ${\mathit{y}}_t={(y_{1,t},\dots ,y_{k,t})}^{\prime }$ is a vector of k assets at time t, $J_t$ is a regime indicator, ${\Phi }_0^{\left(J_t\right)}$ is a k-dimensional vector, ${\Phi }_l^{\left(J_t\right)}$ is a $k\times k$ matrix, ${\Sigma }^{\left(J_t\right)}$ is a $k\times k$ positive-definite matrix with diagonal elements, scalar parameters are satisfied ${\theta }_1^{\left(J_t\right)},{\theta }_2^{\left(J_t\right)}>0$ and $0<{\theta }_1^{\left(J_t\right)}+{\theta }_2^{\left(J_t\right)}<1$, and ${\Psi }_{t-1}$ is a $k\times k$ sample correlation matrix shocks from $t-S,\dots ,t-1$ for a pre-specified S. Moreover, $z_t$ is a hysteresis variable. In this study, we investigate two distinct forms of asymmetric volatility within the framework of a multivariate hysteretic autoregressive (MHAR) model. The first approach incorporates the asymmetric volatility structure proposed by [17] into the MHAR-GARCH framework, resulting in the MHAR-GJR-GARCH model. The second approach introduces the quadratic GARCH specification, as developed by authors of [21], leading to the formulation of the MHAR-QGARCH model. We also derive the volatility dynamics of the MHAR-GJR-GARCH model:

$$\begin{array}{ccc}\hfill h_{i,t}& =& {\omega }_{i0}^{\left(J_t\right)}+\sum _{l=1}^{p_i}{\alpha }_{il}^{\left(J_t\right)}{a}_{i,t-l}^2+\sum _{l=1}^{q_i}{\beta }_{il}^{\left(J_t\right)}{h}_{i,t-l}+\sum _{l=1}^{r_i}{\delta }_{il}^{\left(J_t\right)}{I}_{i,t-l}\{{\mathit{ϵ}}_{t-l}<\mathbf{0}\}{a}_{i,t-l}^2,i=1,\dot{...},k,\hfill \end{array}$$

where $I_i(.)$ is a $k\times 1$ indicator function that returns the value of 1 when the argument is true or 0 otherwise. The volatility of the MHAR - QGARCH model is as follows:

$$\begin{array}{ccc}\hfill h_{i,t}& =& {\omega }_{i0}^{\left(J_t\right)}+\sum _{l=1}^{p_i}{\alpha }_{il}^{\left(J_t\right)}{a}_{i,t-l}^2+\sum _{l=1}^{q_i}{\beta }_{il}^{\left(J_t\right)}{h}_{i,t-l}+\sum _{l=1}^{r_i}{\delta }_{il}^{\left(J_t\right)}{a}_{i,t-l},i=1,\dot{...},k,\hfill \end{array}$$

We now consider the basic cases of two models: the bivariate HAR(1) – GJR – GARCH(1,1) model and the bivariate HAR(1) – QGARCH(1,1) model. We assume that innovations in Equation (1) follow the modified bivariate Student-t distribution (see [9]). In this case, we apply the scale SMN representation (see [2]) to the adapted bivariate Student-t distribution, ${\mathcal{T}}_2^{*}(\mathbf{0},{\Sigma }_t,\mathbf{\nu })$ and we choose $z_t=y_{1,t-d}$. Then, the bivariate HAR(1) - GJR - GARCH(1,1) model is described as follows:

$$\begin{array}{ccc}\hfill {\mathit{y}}_t& =& {\Phi }_0^{\left(J_t\right)}+{\Phi }_1^{\left(J_t\right)}{\mathit{y}}_{t-1}+{\mathit{a}}_t,\hfill \\ \hfill {\mathit{a}}_t& =& \mathrm{diag}(\sqrt{h_{1,t}},\sqrt{h_{2,t}}){\mathit{ϵ}}_t,\hfill \\ \hfill {\mathit{ϵ}}_t|{\Lambda }_t& \sim & N_2(\mathbf{0},{\Lambda }_t^{-1/2}{\Sigma }_t{\Lambda }_t^{-1/2}),\hfill \\ \hfill {\Lambda }_t^{1/2}& =& \mathrm{diag}(\sqrt{{\lambda }_{1,t}},\sqrt{{\lambda }_{2,t}}),{\lambda }_{i,t}\sim Ga\left({\displaystyle \frac{{\nu }_i}{2}},{\displaystyle \frac{{\nu }_i}{2}}\right),i=1,2,\hfill \\ \hfill {\Sigma }_t& =& \left(1-{\theta }_1^{\left(J_t\right)}-{\theta }_2^{\left(J_t\right)}\right){\Sigma }^{\left(J_t\right)}+{\theta }_1^{\left(J_t\right)}{\Sigma }_{t-1}+{\theta }_2^{\left(J_t\right)}{\Psi }_{t-1},\hfill \end{array}$$

(4)

$$\begin{array}{ccc}\hfill J_t& =& \left\{\begin{array}{cc}1\hfill & \mathrm{if}{y}_{1,t-d}\le r_L,\hfill \\ 2\hfill & \mathrm{if}{y}_{1,t-d}>r_U,\hfill \\ J_{t-1}\hfill & \mathrm{otherwise},\hfill \end{array}\right.\hfill \end{array}$$

(5)

with the $(u,v)$th element of ${\Psi }_{t-1}$ is described in (3) and the conditional volitilities as follows:

$$\begin{array}{ccc}\hfill h_{i,t}& =& {\omega }_{i0}^{\left(J_t\right)}+{\alpha }_{i1}^{\left(J_t\right)}{a}_{i,t-1}^2+{\beta }_{i1}^{\left(J_t\right)}{h}_{i,t-1}+{\delta }_{i1}^{\left(J_t\right)}{I}_{i,t-1}\{{\mathit{ϵ}}_{t-1}<\mathbf{0}\}{a}_{i,t-1}^2,i=1,2,\hfill \end{array}$$

where $I_i(.)$ is a $2\times 1$ indicator function that returns the value of 1 when the argument is true or 0 otherwise. The positivity and stationarity conditions for volatility are given as follows:

$$\begin{array}{c}\hfill {\omega }_{i0}^{\left(J_t\right)}>0,{\alpha }_{i1}^{\left(J_t\right)}\ge 0,{\beta }_{i1}^{\left(J_t\right)}\ge 0,{\delta }_{i1}^{\left(J_t\right)}\ge 0\mathrm{and}{\alpha }_{i1}^{\left(J_t\right)}+{\beta }_{i1}^{\left(J_t\right)}+{\delta }_{i1}^{\left(J_t\right)}\le 1.\end{array}$$

(6)

The bivariate HAR(1) - QGARCH(1,1) model modifies the conditional volatilities as follows:

$$\begin{array}{ccc}\hfill h_{i,t}& =& {\omega }_{i0}^{\left(J_t\right)}+{\alpha }_{i1}^{\left(J_t\right)}{a}_{i,t-1}^2+{\beta }_{i1}^{\left(J_t\right)}{h}_{i,t-1}+{\delta }_{i1}^{\left(J_t\right)}{a}_{i,t-1},i=1,2,\hfill \end{array}$$

where the positivity and stationarity conditions for volatility are given as follows:

$$\begin{array}{c}\hfill {\omega }_{i0}^{\left(J_t\right)}>0,{\alpha }_{i1}^{\left(J_t\right)}\ge 0,{\beta }_{i1}^{\left(J_t\right)}\ge 0,{\delta }_{i1}^{\left(J_t\right)}\ge 0,{\left({\delta }_{i1}^{\left(J_t\right)}\right)}^2\le 4(1-{\alpha }_{i1}^{\left(J_t\right)}-{\beta }_{i1}^{\left(J_t\right)}),\\ \hfill \mathrm{and}{\alpha }_{i1}^{\left(J_t\right)}+{\beta }_{i1}^{\left(J_t\right)}\le 1,\end{array}$$

(7)

and we specify the unconditional correlation matrix ${\Sigma }^{\left(J_t\right)}$:

$$\begin{array}{ccc}\hfill {\Sigma }^{\left(J_t\right)}& =& \left[\begin{array}{cc}1& {\rho }^{\left(J_t\right)}\\ {\rho }^{\left(J_t\right)}& 1\end{array}\right].\hfill \end{array}$$

(8)

3. Bayesian Inference

To estimate the unknown parameters of the proposed models in a Bayesian framework, for example, we create groups of the unknown parameters: (i) ${\mathit{\varphi }}_i^{\left(J_t\right)}={({\varphi }_{i0}^{\left(J_t\right)},{\varphi }_{i1}^{\left(J_t\right)},{\varphi }_{i2}^{\left(J_t\right)})}^{\prime },i,J_t=1,2$; (ii) $\mathbf{r}={(r_L,r_U)}^{\prime }$; (iii) $\mathbf{\nu }={({\nu }_1,{\nu }_2)}^{\prime }$; (iv) $\mathbf{\rho }={({\rho }^{\left(1\right)},{\rho }^{\left(2\right)})}^{\prime }$; (v) ${\mathbf{\gamma }}_i^{\left(J_t\right)}={({\omega }_{i0}^{\left(J_t\right)},{\alpha }_{i1}^{\left(J_t\right)},{\beta }_{i1}^{\left(J_t\right)},{\delta }_{i1}^{\left(J_t\right)})}^{\prime }$, $i,J_t=1,2$; (vi) ${\mathbf{\eta }}^{\left(J_t\right)}={({\theta }_1^{\left(J_t\right)},{\theta }_2^{\left(J_t\right)})}^{\prime },$ and (vii) d. We define $\mathit{\theta }$ as a vector of all unknown parameters of the proposed model. Following that, the bivariate HAR(1) - GJR - GARCH(1,1) and bivariate HAR(1) - QGARCH(1,1) models’ conditional likelihood functions are given by:

$$\begin{array}{ccc}\hfill ln\left(\mathcal{L}\left(\mathcal{y}|\mathit{\theta }\right)\right)& \propto & \sum _t\left\{\sum _{J_t=1}^2\right[-0.5ln\left({\displaystyle \frac{h_{1,t}{h}_{2,t}(1-{{\rho }^{\left(J_t\right)}}^2)}{{\lambda }_{1,t}{\lambda }_{2,t}}}\right)\hfill \\ \hfill & & -{\displaystyle \frac{1}{2(1-{{\rho }^{\left(J_t\right)}}^2)}}\left({\displaystyle \frac{{\lambda }_{1t}{a}_{1t}^2}{h_{1,t}}}+{\displaystyle \frac{{\lambda }_{2t}{a}_{2t}^2}{h_{2,t}}}-{\displaystyle \frac{2{\rho }^{\left(J_t\right)}{a}_{1t}{a}_{2t}\sqrt{{\lambda }_{1,t}{\lambda }_{2,t}}}{\sqrt{h_{1,t}}\sqrt{h_{2,t}}}}\right)\left]\right\},\hfill \end{array}$$

(9)

where ${\mathit{a}}_t={\mathit{y}}_t-{\mathbf{\Phi }}_0^{\left(J_t\right)}-{\mathbf{\Phi }}_1^{\left(J_t\right)}{\mathit{y}}_{t-1}$.

We set up prior distributions for the unknown parameters. Assume that ${\mathit{\varphi }}_i^{\left(J_t\right)}\sim N_3({\mu }_{0i},{\Sigma }_{0i}^{-1}),$$i,J_t=1,2$; for parameter threshold $r_L\sim Unif(l_1,l_2)$, where $l_1$ and $l_2$ are the pth and $(100-2p)$th percentiles of observed time series, respectively, for $0<p<33.$ Furthermore, $r_U|r_L\sim Unif(u_1,u_2),$ where $u_2$ is the $(100-p)$th percentile and $u_1=r_L+c^{*}$ for $c^{*}$ is a selected number that ensures $r_L+c^{*}\le r_U$ and at least $p\%$ of observations are in the range $(r_L,r_U)$. For degrees of freedom, we assume the scale mixture variables ${\lambda }_{i,t}\sim Ga\left({\nu }_i/2,{\nu }_i/2\right)$ and ${\nu }_i\sim \mathrm{Unif}(2,60),i=1,2$, and ${\rho }^{\left(J_t\right)}\sim \mathrm{Unif}(-1,1)$. For lag d, we choose the discrete uniform prior $p\left(d\right)=1/d_0$ with maximum delay $d_0$. In terms of volatility parameters, ${\mathbf{\gamma }}_i^{\left(J_t\right)}$ follows a uniform distribution such that ${\mathbf{\gamma }}_i^{\left(J_t\right)}$ is proportional to $I\left(S_1\right)$ or $I\left(S_2\right)$, where $S_1$ and $S_2$ are the sets that satisfy (6) and (7), respectively.

The conditional posterior distribution for each group is proportional to the conditional likelihood function multiplied by the prior density for that group, as shown below:

$$P\left({\mathit{\theta }}_l|\mathit{y},{\mathit{\theta }}_{-l}\right)\propto \mathcal{L}\left(\mathit{y}|\mathit{\theta }\right)P\left({\mathit{\theta }}_l|{\mathit{\theta }}_{-l}\right),$$

where ${\mathit{\theta }}_l$ is each parameter group, $P\left({\mathit{\theta }}_l\right)$ is its prior density, and ${\mathit{\theta }}_{-l}$ is the vector of all parameters, except ${\mathit{\theta }}_l$. The conditional posterior distribution of delay lag d follows a multinomial distribution with a probability:

$$\begin{array}{c}\hfill Pr(d=i|\mathit{y},{\mathit{\theta }}_{-d})={\displaystyle \frac{p\left(\mathit{y}\right|d=i,{\mathit{\theta }}_{-d})Pr(d=i)}{{\sum }_{j=1}^{d_0}p\left(\mathit{y}\right|d=j,{\mathit{\theta }}_{-d})Pr(d=j)}},i=1,\dots ,d_0.\end{array}$$

In this study, with the exception of the lag parameter d, the conditional posterior distributions of the remaining parameter groups exhibit non-standard forms. To make statistical inference, we employ an adaptive Markov Chain Monte Carlo (MCMC) method for selected parameter groups, complemented by a random walk Metropolis algorithm. Specifically, we assume that the innovation term in Equation (4) follows a Gaussian distribution, which serves as the kernel for sampling ${\mathit{\varphi }}_i^{\left(J_t\right)}$. For the parameter groups ${\mathit{\eta }}^{\left(J_t\right)}$ and ${\mathbf{\gamma }}_i^{\left(J_t\right)}$, where $i,J_t=1,2$, an adaptive MCMC approach is utilized to draw samples, whereas the remaining parameters are updated using the random walk Metropolis algorithm. The detailed procedures of the adaptive Metropolis-Hastings MCMC algorithm are thoroughly presented by authors of [6] and [9], where the authors provide a comprehensive framework for its implementation and application. Based on guidelines of [6], we further manage a scale matrix to attain ideal acceptance rates of $20\%$ to $60\%$.

In a Bayesian framework, we need to set up the initial values for each parameter group. For autoregressive coefficient parameters, ${\mathit{\varphi }}_i^{\left(J_t\right)}=(0,0,0)$; for degrees of freedom ${\nu }_i^{\left(J_t\right)}=50$; ${\rho }_i^{\left(J_t\right)}=0.5$; $d_0=3$; ${\omega }_{i0}^{\left(J_t\right)}={\alpha }_{i1}^{\left(J_t\right)}={\beta }_{i1}^{\left(J_t\right)}={\delta }_{i1}^{\left(J_t\right)}=0.1$; and $({\theta }_1^{\left(J_t\right)},{\theta }_2^{\left(J_t\right)})=(0.05,0.2)fori,J_t=1,2$; thresholds $r_L$ and $r_U$ are established at the 33rd and 67th percentiles, respectively; and we set $p=20$ to make certain of enough observations in each regime for a valid inference. For the remainder of the analysis, we specify $S=3$.

4. Forecasting Marginal Expected Shortfall and Value at Risk

Value-at-Risk (VaR) and Marginal Expected Shortfall (MES) are now considered systemic risk assessments for financial risk management. The authors of [1] define MES as a financial firm’s marginal contribution to the financial system’s expected shortfall. The authors of [4] define MES at the alpha level for a financial institution at time t given ${\mathcal{F}}_{t-1}$ as follows:

$$\begin{array}{ccc}\hfill & & {MES}_t=E\left[y_{2,t}|y_{1,t}<{VaR}_{1,t}\left(\alpha \right);{\mathcal{F}}_{t-1}\right],\hfill \end{array}$$

(10)

where ${VaR}_{1,t}\left(\alpha \right)$ is the VaR of $y_{1,t}$ at the $\alpha$-level such that $P\left(y_{1,t}<{VaR}_{1,t}\left(\alpha \right)\mid {\mathcal{F}}_{t-1}\right)=\alpha$. Here, $y_{2,t}$ stands for the stock return of a financial institution, whereas $y_{1,t}$ stands for the market return.

To produce ${MES}_t$, we estimate one-step-ahead quantiles and volatilities for $y_{1,n+1}$ from the investigated model described in (4) with forecast origin $t=n$. We get quantiles from the posterior predictive distribution, which is:

$$p\left({\mathit{y}}_{n+1}\right|{\mathcal{F}}_n)={\int }_{\mathit{\theta }\in \Theta }p\left({\mathit{y}}_{n+1}\right|{\mathcal{F}}_n,\mathit{\theta })p\left(\mathit{\theta }\right|{\mathcal{F}}_n)d\mathit{\theta }.$$

(11)

Suppose that $\{{\mathit{\theta }}^{\left[r\right]},r=n_0+1,\dots ,N\}$ are rth MCMC draw from the posterior density $p\left(\mathit{\theta }\right|{\mathcal{F}}_n)$ after the $n_0$ burn-in sample. Thus, we can sample $\{{\mathit{y}}_{n+1}^{\left[r\right]},r=n_0+1,\dots ,N\}$ from the marginal predictive distribution, $p\left({\mathit{y}}_{n+1}\right|{\mathcal{F}}_n)$, by sampling the following conditional distribution:

$$\begin{array}{ccc}\hfill {\mathit{y}}_{n+1}^{\left[r\right]}|{\mathcal{F}}_n,{\mathit{\theta }}^{\left[r\right]}& \sim & {\mathcal{T}}_2^{*}({\mathbf{\mu }}_{n+1}^{\left[r\right]},{\Sigma }_{n+1}^{*\left[r\right]},{\mathbf{\nu }}^{\left[r\right]})\hfill \\ \hfill {\mathcal{T}}_2^{*}\left({\mathit{y}}_{n+1}^{\left[r\right]}\right|{\mathbf{\mu }}_{n+1}^{\left[r\right]},{\Sigma }_{n+1}^{*\left[r\right]},{\mathbf{\nu }}^{\left[r\right]})& =& {\int }_0^{\infty }{\int }_0^{\infty }{N}_2\left({\mathit{y}}_{n+1}^{\left[r\right]}\left|{\mathbf{\mu }}_{n+1}^{\left[r\right]},{\mathbf{\lambda }}_{n+1}^{-1/2}{\Sigma }_{n+1}^{*\left[r\right]}{\mathbf{\lambda }}_{n+1}^{-1/2}\right.\right)\hfill \\ & \times & \prod _{i=1}^{2}Ga\left({\lambda }_{i,n+1}\left|{\displaystyle \frac{{\nu }_i^{\left[r\right]}}{2}},{\displaystyle \frac{{\nu }_i^{\left[r\right]}}{2}}\right.\right)d{\lambda }_{1,n+1}d{\lambda }_{2,n+1},\hfill \end{array}$$

(12)

where ${\mathbf{\mu }}_{n+1}^{\left[r\right]}$ and ${\Sigma }_{n+1}^{*\left[r\right]}=\mathrm{diag}(h_{1,n+1}^{\left[r\right]},h_{2,n+1}^{\left[r\right]}){\Sigma }_{n+1}^{\left[r\right]}$ are a conditional mean and covariance of $p\left({\mathit{y}}_{n+1}\right|{\mathcal{F}}_n,\mathit{\theta })$ at the rth iteration. To assess the correctness of a VaR performance, we calculate the violation rate (VRate). The accuracy of a VaR performance is verified by recording the failure rate; that is, the violation rate:

$$\begin{array}{c}\hfill \mathrm{Violation}\mathrm{rate}={\displaystyle \frac{1}{h_0}}\sum _{t=n+1}^{n+h_0}I(r_t<-{\mathrm{VaR}}_t)\end{array}$$

(13)

where $h_0$ is the out-of-sample period size and $r_t$ is the return at time t. We use two tests to assess the validity of the VaR forecasts: the conditional coverage (CC) test created by [11] and the unconditional coverage (UC) test created by [18]. The CC test is conducted to investigate the null hypothesis that the violations are independently distributed, whereas the UC test is suggested to determine whether the percentage of violations is equivalent to the VaR significance level.

5. Simulation Study

In order to access the effectiveness of the Bayesian approach, we run two simulations of the suggested models in this section. Model 1 is the bivariate HAR(1)- GJR - GARCH(1,1) model and Model 2 is the bivariate HAR(1) - QGARCH(1,1) model. Model 1 is given by::

$$\begin{array}{ccc}\hfill {\mathit{y}}_t& =& {\Phi }_0^{\left(J_t\right)}+{\Phi }_1^{\left(J_t\right)}{\mathit{y}}_{t-1}+{\mathit{a}}_t,\hfill \\ \hfill {\mathit{a}}_t& =& diag\left(\sqrt{h_{1,t}},\sqrt{h_{2,t}}\right){\mathit{ϵ}}_t,{\mathit{ϵ}}_t\sim {\mathcal{T}}_2^{*}\left(\mathbf{0},{\Sigma }_t,\mathbf{\nu }\right)\hfill \\ \hfill {\mathit{ϵ}}_t|{\Lambda }_t& \sim & N_2(\mathbf{0},{\Lambda }_t^{-1/2}{\Sigma }_t{\Lambda }_t^{-1/2}),\hfill \end{array}$$

(14)

$$\begin{array}{c}{\Phi }_0^{\left(1\right)}=\left[\begin{array}{c}{\varphi }_{10}^{\left(1\right)}\\ {\varphi }_{20}^{\left(1\right)}\end{array}\right]=\left[\begin{array}{c}-0.10\\ -0.10\end{array}\right],{\Phi }_1^{\left(1\right)}=\left[\begin{array}{cc}{\varphi }_{11}^{\left(1\right)}& {\varphi }_{12}^{\left(1\right)}\\ {\varphi }_{21}^{\left(1\right)}& {\varphi }_{22}^{\left(1\right)}\end{array}\right]=\left[\begin{array}{cc}0.20& 0.25\\ 0.25& 0.30\end{array}\right]\hfill \\ \hfill \\ {\Phi }_0^{\left(2\right)}=\left[\begin{array}{c}{\varphi }_{10}^{\left(2\right)}\\ {\varphi }_{20}^{\left(2\right)}\end{array}\right]=\left[\begin{array}{c}-0.08\\ -0.15\end{array}\right],{\Phi }_1^{\left(2\right)}=\left[\begin{array}{cc}{\varphi }_{11}^{\left(2\right)}& {\varphi }_{12}^{\left(2\right)}\\ {\varphi }_{21}^{\left(2\right)}& {\varphi }_{22}^{\left(2\right)}\end{array}\right]=\left[\begin{array}{cc}0.30& 0.35\\ 0.35& 0.30\end{array}\right]\hfill \end{array}$$

$$\begin{array}{ccc}\hfill with{h}_{1,t}& =& \left\{\begin{array}{cc}0.07+0.20a_{1,t-1}^2+0.20h_{1,t-1}+0.40I_{1,t-1}\{{\mathit{ϵ}}_{t-1}<0\}{a}_{1,t-1}^2\hfill & \mathrm{if}{J}_t=1,\hfill \\ 0.03+0.20a_{1,t-1}^2+0.25h_{1,t-1}+0.55I_{1,t-1}\{{\mathit{ϵ}}_{t-1}<0\}{a}_{1,t-1}^2\hfill & \mathrm{if}{J}_t=2,\hfill \end{array}\right.\hfill \end{array}$$

(15)

$$\begin{array}{ccc}\hfill h_{2,t}& =& \left\{\begin{array}{cc}0.04+0.25a_{2,t-1}^2+0.10h_{2,t-1}+0.40I_{2,t-1}\{{\mathit{ϵ}}_{t-1}<0\}{a}_{2,t-1}^2\hfill & \mathrm{if}{J}_t=1,\hfill \\ 0.02+0.30a_{2,t-1}^2+0.15h_{2,t-1}+0.40I_{2,t-1}\{{\mathit{ϵ}}_{t-1}<0\}{a}_{2,t-1}^2\hfill & \mathrm{if}{J}_t=2,\hfill \end{array}\right.\hfill \end{array}$$

(16)

$$\begin{array}{ccc}\hfill {\Sigma }_t& =& \left\{\begin{array}{cc}(1-0.40-0.10){\Sigma }^{\left(1\right)}+0.40{\Sigma }_{t-1}+0.10{\Psi }_{t-1}\hfill & \mathrm{if}{J}_t=1,\hfill \\ (1-0.50-0.20){\Sigma }^{\left(2\right)}+0.50{\Sigma }_{t-1}+0.20{\Psi }_{t-1}\hfill & \mathrm{if}{J}_t=2,\hfill \end{array}\right.\hfill \end{array}$$

(17)

$$\begin{array}{ccc}\hfill where{J}_t& =& \left\{\begin{array}{ccc}1\hfill & \mathrm{if}\hfill & y_{1,t-d}<-0.5,\hfill \\ J_{t-1}\hfill & \mathrm{if}\hfill & -0.45\le y_{1,t-d}\le 0.1,\hfill \\ 2\hfill & \mathrm{if}\hfill & y_{1,t-d}>0.1,\hfill \end{array}\right.\hfill \\ \hfill and{\Sigma }^{\left(1\right)}& =& \left[\begin{array}{cc}1& 0.65\\ 0.65& 1\end{array}\right],{\Sigma }^{\left(2\right)}=\left[\begin{array}{cc}1& 0.8\\ 0.8& 1\end{array}\right],\mathbf{\nu }=\left[\begin{array}{c}{\nu }_1\\ {\nu }_2\end{array}\right]=\left[\begin{array}{c}8.0\\ 10.0\end{array}\right],\mathrm{and}d=1.\hfill \end{array}$$

(18)

Model 2 is describled as follows:

$$\begin{array}{ccc}\hfill {\mathit{y}}_t& =& {\Phi }_0^{\left(J_t\right)}+{\Phi }_1^{\left(J_t\right)}{\mathit{y}}_{t-1}+{\mathit{a}}_t,\hfill \\ \hfill {\mathit{a}}_t& =& diag\left(\sqrt{h_{1,t}},\sqrt{h_{2,t}}\right){\mathit{ϵ}}_t,{\mathit{ϵ}}_t\sim {\mathcal{T}}_2^{*}\left(\mathbf{0},{\Sigma }_t,\mathbf{\nu }\right)\hfill \end{array}$$

(19)

$$\begin{array}{c}{\Phi }_0^{\left(1\right)}=\left[\begin{array}{c}{\varphi }_{10}^{\left(1\right)}\\ {\varphi }_{20}^{\left(1\right)}\end{array}\right]=\left[\begin{array}{c}-0.10\\ -0.08\end{array}\right],{\Phi }_1^{\left(1\right)}=\left[\begin{array}{cc}{\varphi }_{11}^{\left(1\right)}& {\varphi }_{12}^{\left(1\right)}\\ {\varphi }_{21}^{\left(1\right)}& {\varphi }_{22}^{\left(1\right)}\end{array}\right]=\left[\begin{array}{cc}0.32& 0.30\\ 0.37& 0.35\end{array}\right]\hfill \\ \hfill \\ {\Phi }_0^{\left(2\right)}=\left[\begin{array}{c}{\varphi }_{10}^{\left(2\right)}\\ {\varphi }_{20}^{\left(2\right)}\end{array}\right]=\left[\begin{array}{c}-0.08\\ -0.08\end{array}\right],{\Phi }_1^{\left(2\right)}=\left[\begin{array}{cc}{\varphi }_{11}^{\left(2\right)}& {\varphi }_{12}^{\left(2\right)}\\ {\varphi }_{21}^{\left(2\right)}& {\varphi }_{22}^{\left(2\right)}\end{array}\right]=\left[\begin{array}{cc}0.35& 0.30\\ 0.33& 0.37\end{array}\right]\hfill \end{array}$$

$$\begin{array}{ccc}\hfill \mathrm{with}{h}_{1t}& =& \left\{\begin{array}{cc}0.07+0.20a_{1,t-1}^2+0.10h_{1,t-1}+0.40a_{1,t-1}\hfill & \mathrm{if}{J}_t=1,\hfill \\ 0.03+0.30a_{1,t-1}^2+0.10h_{1,t-1}+0.35a_{1,t-1}\hfill & \mathrm{if}{J}_t=2,\hfill \end{array}\right.\hfill \end{array}$$

(20)

$$\begin{array}{ccc}\hfill h_{2t}& =& \left\{\begin{array}{cc}0.04+0.40a_{2,t-1}^2+0.05h_{2,t-1}+0.30a_{2,t-1}\hfill & \mathrm{if}{J}_t=1,\hfill \\ 0.02+0.30a_{2,t-1}^2+0.10h_{2,t-1}+0.20a_{2,t-1}\hfill & \mathrm{if}{J}_t=2,\hfill \end{array}\right.\hfill \end{array}$$

(21)

$$\begin{array}{ccc}\hfill {\Sigma }_t& =& \left\{\begin{array}{cc}(1-0.40-0.35){\Sigma }^{\left(1\right)}+0.40{\Sigma }_{t-1}+0.35{\Psi }_{t-1}\hfill & \mathrm{if}{J}_t=1,\hfill \\ (1-0.55-0.15){\Sigma }^{\left(2\right)}+0.55{\Sigma }_{t-1}+0.15{\Psi }_{t-1}\hfill & \mathrm{if}{J}_t=2,\hfill \end{array}\right.\hfill \end{array}$$

(22)

$$\begin{array}{ccc}\hfill \mathrm{where}{J}_t& =& \left\{\begin{array}{ccc}1\hfill & \mathrm{if}\hfill & y_{1,t-d}<-0.45,\hfill \\ J_{t-1}\hfill & \mathrm{if}\hfill & -0.45\le y_{1,t-d}\le 0.1,\hfill \\ 2\hfill & \mathrm{if}\hfill & y_{1,t-d}>0.1,\hfill \end{array}\right.\hfill \\ \hfill \mathrm{and}{\Sigma }^{\left(1\right)}& =& \left[\begin{array}{cc}1& 0.5\\ 0.5& 1\end{array}\right],{\Sigma }^{\left(2\right)}=\left[\begin{array}{cc}1& 0.85\\ 0.85& 1\end{array}\right],\mathbf{\nu }=\left[\begin{array}{c}{\nu }_1\\ {\nu }_2\end{array}\right]=\left[\begin{array}{c}8.0\\ 10.0\end{array}\right],\mathrm{and}d=1.\hfill \end{array}$$

(23)

Models 1 and 2 are created utilizing the actual values shown in Table 1 and Table 2. For each time series, we set up the sample size $n=2,000$. We carry out $N=30,000$ MCMC iterations and discard the first $M=10,000$ burn-in iterates. For the hyper - parameters, we choose the initial values for all parameters of the investigated model to be ${\mathbf{\mu }}_{0i}=0$, $\mathrm{diag}\left({\Sigma }_{0i}\right)=0.1$, $l_1={\wp }_{20}$, $u_2={\wp }_{80}$, $c_L=2$, $c_U=60$, and $d_0=3$.

Results for the parameter estimates of the simulation study are shown in Table 1 and Table 2. The tables present the posterior means, medians, standard deviations, and $95\%$ credible ranges for Models 1-2 over the 200 replications. We observe that the $95\%$ credible interval contains the corresponding true value for each parameter. The posterior means and medians in each case are fairly close to the true parameter values. The posterior modes of lag d are demonstrated, and it can be explained that the posterior mode of d provides a reliable estimate of the delay parameter because the posterior probability for $d=1$ is nearly equal to one. To check the convergence of MCMC, we examine the ACF plots of all coefficients. For compactness, we present only the autocorrelation function (ACF) plots based on Model 2, omitting ACF plots of Model 1 to conserve space. Figure 1 and Figure 2 provide additional evidence supporting the convergence of the MCMC algorithm. Based on these diagnostic checks, we conclude that the proposed models are well-suited for implementation within the Bayesian framework.

6. Emperical Study

The empirical analysis in this study is based on daily closing prices of four major financial indices: the S&P500, Bank of America (BAC), Intercontinental Exchange (ICE), and Goldman Sachs (GS). The data span from January 4, 2006, to December 30, 2021, encompassing a total of 4026 trading days. These data were retrieved from Yahoo Finance, a widely recognized source for historical market data. To construct the return series, we compute the continuously compounded returns (log-returns) using the formula $y_t=log\left(p_t\right)-log\left(p_{t-1}\right)$, where $p_t$ denotes the asset’s closing price at time t.

Table 3 defines three target datasets: DS1 {S&P500, GS}, DS2 {S&P500, ICE}, and DS3 {S&P500, BAC}. It also presents the descriptive statistics of the corresponding return series, along with the results of two multivariate normality tests: Mardia’s test and the Henze-Zirkler test (see [12,13]). The return distributions are clearly skewed and have high kurtosis, especially showing strong positive skewness. Due to the noticeable asymmetry and the presence of heavy tails in the return data, we recommend using asymmetric models with fat-tailed multivariate error distributions instead of models that assume multivariate normal errors. Figure 3 presents the time series plot of daily returns for the selected financial assets. As shown, the sample period spans several significant market events, notably the Global Financial Crisis, which officially began on September 15, 2008, following the bankruptcy of Lehman Brothers. For the purpose of estimation and out-of-sample evaluation, the dataset is divided into two distinct sub-periods. The first segment, consisting of 3726 daily observations, is used to estimate the model parameters. The remaining 300 observations are reserved for out-of-sample forecasting and performance assessment.

This section’s hyper-parameters correspond to those in the simulation study. Table 4, Table 5, Table 6 and Table 7 present a summary of Bayesian estimates for three datasets for the BHAR(1) - GJR - GARCH(1,1) and the BHAR(1) - QGARCH(1,1) models. The significant value of ${\varphi }_{12}^{\left(1\right)}$ in Table 4 and Table 6 indicate that the performance of the previous day’s return of Goldman Sachs stock has a considerable negative impact on the S&P 500’s returns in the lower regimes. We can see that the parameter estimates for ${\varphi }_{12}^{\left(1\right)}$, ${\varphi }_{12}^{\left(2\right)}$, ${\varphi }_{21}^{\left(1\right)}$, and ${\varphi }_{21}^{\left(2\right)}$ are identical in both fitted models when we look at datasets DS2 and DS3 in Table 5 and Table 7. To assess the validity of the proposed models, we further employ the Geweke convergence diagnostic (see [16]). The p-values reported in Table 8 and Table 9 suggest that the MCMC chains generated from the models have converged. As there is no statistical evidence of non-convergence, we conclude that the proposed models are appropriately specified and reliable for inference.

To evaluate the accuracy of the models using Value-at-Risk (VaR), we present VaR forecasts along with the results of VaR backtesting at the 1% and 5% significance levels. Specifically, Table 10 and Table 11 report the p-values of the Unconditional Coverage (UC) and Conditional Coverage (CC) tests for the two proposed models: bivariate HAR(1)-GJR-GARCH(1,1) and HAR(1)-QGARCH(1,1) as well as the benchmark bivariate HAR(1)-GARCH(1,1) model. When evaluating DS1, DS2, and DS3 across the three models, the violation rates (VRate) for the S&P 500 tend to be significantly higher than the nominal 1% level, suggesting a slight underestimation of tail risk. In contrast, the VRates for Bank of America (BAC), Intercontinental Exchange (ICE), and Goldman Sachs (GS) indicate a tendency toward risk overestimation. Nevertheless, the backtesting results show that all three models perform adequately as risk models. At the 5% significance level, both the proposed models and the benchmark BHAR(1)-GARCH(1,1) model yield UC and CC test p-values above 5%, indicating no statistical evidence of model misspecification. These findings confirm that the proposed models deliver reliable and independent risk forecasts. Figure 4 and Figure 5 display the VaR forecasts based on the bivariate HAR(1)-GJR - GARCH(1,1) and HAR(1)-QGARCH(1,1) models, which clearly show that the models are capable of identifying volatility spikes in returns, despite infrequent violations of the forecast bounds.

To understand how well the proposed models can capture the expected shortfall movement, we present the backtesting measures of the MES forecasts proposed by authors of [14] in Table 10. The model with the smallest values in the boxes is the best. These findings indicate that the proposed models are the best.

7. Conclusions

This paper investigates the MHAR–GJR–GARCH and MHAR – QGARCH models by incorporating asymmetric volatility dynamics, conditional correlations, and a hysteresis variable to control regime switching and dynamic delays. Bayesian inference is employed for efficient estimation of the model parameters. A comparative analysis of backtesting results for the VaR and MES forecasts is conducted. We also include the benchmark model MHAR - GARCH with adapted multivariate Student-t errors and compare backtesting measures of Value-at-Risk (VaR) and Marginal Expected Shortfall (MES) forecasts. The backtesting measures indicate that, in general, the proposed models demonstrate reliable capabilities in capturing tail risk behavior and delivering accurate risk predictions. Notably, the proposed models deliver significantly improved performance over the original MHAR–GARCH errors model, particularly in capturing asymmetric tail risks and providing more accurate risk forecasts.