A time-varying GARCH mixed-effects model for isolating high- and low- frequency volatility and co-volatility

Abstract

This article studies long-term, short-term volatility and co-volatility in stock markets by introducing modelling strategies to the multivariate data analysis that deal with serially correlated innovations and cross-section dependence. In particular, it presents an innovative mixed-effects model through a GARCH process, allowing for heterogeneity effects and time-series dynamics. We propose a non-parametric regression model of the penalized low-rank smoothing spline to present time trends into the variance and covariance equations. The strategy provides flexible modelling of the low-frequency volatility and co-volatility in equity markets. The decomposed low-frequency matrix was modelled using the modified Cholesky factorization. The Hamiltonian Monte Carlo technique is implemented as a Bayesian computing process for estimating parameters and latent factors. The advantage of our modelling strategy in empirical studies is highlighted by examining the effect of latent financial factors on a panel across 10 equities over 110 weekly series. The model can differentiate non-parametrically dynamic patterns of high and low frequencies of variance–covariance structural equations and incorporate economic features to predict variabilities in stock markets regarding time-series evidence.

1 Introduction

Most effort in the literature of financial and economic multivariate data with serially correlated time-series has been devoted to the slow-moving variation and covariation topics. Related studies have frequently been centred on the long-run variation. For example, Engle and Rangel (2008) introduced a notable spline-GARCH model for low-frequency volatility and global macroeconomic events. Audrino and Bühlmann (2009) offered a spline regression model in the non-parametric setting to predict the volatility of financial time-series. Some useful strategies rendering high-performance model fittings are introduced to mixed modelling with smoothing splines (Wahba (1990); Eilers and Marx (1996); Wand (2003)). Speckman and Sun (2003) applied full Bayesian smoothing splines, thin-plate splines, L-splines and several models using intrinsic autoregressive priors. Vogt and Linton (2017) investigated the longitudinal data analysis with non-parametric regression functions in the case of varying individual-specific effects.

The treatment of volatilities is conventionally based on variances, while the modelling of covariance functions can provide better representation. It is reasonable since the volatility in a market can cause volatility in other markets. Thus, most research recommends fitting the multivariate generalized autoregressive conditional heteroscedastic (MGARCH) model. It can describe the volatility in financial markets and efficiently estimate large time-varying covariance matrices. Rangel and Engle (2012) characterize the correlation matrix's high- and low-frequency components by combining a factor model with other specifications that capture the dynamic behaviour of volatilities and covariances between a common factor and idiosyncratic returns. Engle and Sokalska (2012) introduced an intra day volatility forecasting model to some assets by decomposing the volatility of high-frequency returns into several components.

Extensions of multivariate GARCH models to cross-sectional time-series data appeal to further investigations since they can conveniently consider volatilities and co-volatilities in correlated cross-section financial data. So far, the existing models’ methodology has been concentrating on volatility using the quadratic-spline method. There is a gap in the literature for flexible low-frequency models to deliver a more precise prediction. These issues motivate us to introduce the spline-GARCH mixed-effects model by extending the univariate case of Engle and Rangel (2008). We study low-frequency changes using semi-parametric mixed-effects models having a time-varying conditional covariance matrix. We let parameters across assets heterogeneous to control possibly correlated, time-invariant and heterogeneity effects.

Frequently, the volatility of an asset affects other assets through conditional covariances over time. We expect that a shock in a market increases the volatility of other associated markets. However, the influence of negative and positive shocks does not lead to the same fluctuation. Our study aims to present the cross-asset dependence between asset returns, which changes over time intervals, and the amount of volatility during a specific period. Our findings reveal that cross-asset dependence may not increase in the long term. It may occur because of the globalization of financial markets. The dependence effect can be investigated directly using cross-section and time-series modelling specified by the dynamic behaviour of covariances or correlations.

In addition to preserving appropriate features of the GARCH model, its time-varying extension can sufficiently isolate high- and low-frequency volatility and co-volatility. This extension leads to coherent predictions since short-term volatility immediately affects the log-return, while long-term volatility possesses a significant role in future transactions. It also helps financial traders with proper forecasting decisions.

We organize the article as follows. In Section 2, we introduce the mixed-effects spline MGARCH model. Section 3 implements an advanced numerical estimation method, so-called Hamiltonian or Hybrid Monte Carlo (HMC), to provide inference. In Section 4, we conducted simulation studies to investigate the main properties of our model and presented an empirical study in predicting low- and high-frequency volatility of log-returns data. Section 5 includes concluding remarks.

2 Specification of the proposed model

We first specify the heterogeneous MGARCH model and then extend the mixed-effects methodology to time-varying MGARCH models aiming to describe low-frequency volatilities and co-volatilities.

2.1 Heterogeneous MGARCH model

For a cross-section of $N$ assets and $T$ periods, let the $r_{\mathit{it}}$ , for $i=1,2,\dots ,N$ and $t=1,2,\dots ,T$ be the return at time $t$ on the $i$ th asset. Assume that the unexpected returns follow the multivariate GARCH model

$${ϵ}_t={\mathrm{r}}_t-E({\mathrm{r}}_t|I_{t-1})={\displaystyle {\sum }_t^{1/2}{e}_t},$$

(2.1)

where the innovation term ${\mathbf{e}}_t$ follows the multivariate normal distribution with zero-mean vector and identity covariance matrix. Then, the unexpected returns series ${\mathbf{\epsilon }}_t$ are normally distributed with zero means and heterogeneous variances and covariances ${\mathrm{\Sigma }}_t=\left\{{\sigma }_{\mathit{ij},t}\right\}$ for all $i,j=1,\dots ,N$ . Assume that there is no autocorrelation and no contemporaneous cross-sectional correlations in innovations ${\mathit{e}}_t$ . Consider the following equations for variances and covariances

$${\sigma }_{ii,t}={\omega }_{ii}+\gamma {ϵ}_{i,t-1}^2+\hspace{0.17em}\delta {\sigma }_{ii,t-1},$$

(2.2)

$${\sigma }_{ij,t}={\omega }_{ij}+\lambda {ϵ}_{i,t-1}{ϵ}_{j,t-1}+\theta {\sigma }_{ij,t-1},$$

(2.3)

for $t>1$ and $i\ne j$ , where ${\omega }_{\mathit{ii}}>0$ and ${\omega }_{\mathit{ij}}\in \mathbb{ℝ}$ . A convenient covariance matrix for the start-up process is ${\mathrm{\Sigma }}_1=\left\{{\sigma }^2{\omega }_{\mathit{ij}}\right\}$ , where ${\sigma }^2$ is a constant positive parameter. This setting accounts for the correlation between the elements of ${\mathrm{\Sigma }}_1$ and subsequent components (2.2) and (2.3). Also, it provides valid inference since the estimation process uses information of all periods. The sufficient conditions for the conditional variance and covariance processes to converge a positive definite (PD) matrix are that $\gamma +\delta <1$ and $\lambda +\theta <1$ . Also, under conditions $\gamma ⩾\lambda$ , $\delta ⩾\theta$ , and positive definiteness of random-effects matrix $\left\{{\omega }_{\mathit{ij}}\right\}$ , the matrix ${\mathrm{\Sigma }}_t$ will be PD for all $t$ .

2.2 Modelling techniques to isolate volatility and co-volatility

We now extend mixed-effects modelling of the heterogeneous MGARCH using the spline-GARCH technique while describing low-frequency volatilities and co-volatilities in the equities market. Commonly, stock returns possess time-varying variances and covariances due to daily news intensity that (a) depend on the macroeconomic events or time-varying conditions that some observed covariates can explain, (b) respond to several latent factors that cannot often be recognized, at least without having additional information. Flexible modelling of the low-frequency volatility and co-volatility, say ${\tau }_{\mathit{ij},t}$ , as variance and covariance components, can deliver the information of all observed and latent factors. It motivated us to develop the non-parametric model fittings as follows. First, we reparameterize (2.2) and (2.3) to the decomposition form

$${\sigma }_{ii,t}=h_{ii,t}{\tau }_{ii,t},$$

(2.4)

where the ${\tau }_{\mathit{ii},t}$ are the low-frequency volatility components. The parameter $h_{\mathit{ii},t}$ , for each $i$ , is a non-negative time-series, for example, GARCH(1,1) with unit unconditional mean, expressed as

$$\begin{array}{cc}{h}_{\mathit{ii},t}={\omega }_{\mathit{ii}}+\gamma \frac{{\epsilon }_{i,t-1}^2}{{\tau }_{\mathit{ii},t-1}}+\delta h_{\mathit{ii},t-1}.\hfill & \hfill \end{array}$$

(2.5)

Similarly,

$$\begin{array}{cc}{\sigma }_{\mathit{ij},t}=h_{\mathit{ij},t}{\tau }_{\mathit{ij},t},\hfill & \hfill \end{array}$$

(2.6)

for $i\ne j$ , where the ${\tau }_{\mathit{ijt}}$ are the low-frequency co-volatility components and

$$\begin{array}{cc}{h}_{\mathit{ij},t}={\omega }_{\mathit{ij}}+\lambda \frac{{\epsilon }_{i,t-1}{\epsilon }_{j,t-1}}{\sqrt{{\tau }_{\mathit{ii},t-1}{\tau }_{\mathit{jj},t-1}}}+\theta h_{\mathit{ij},t-1}.\hfill & \hfill \end{array}$$

(2.7)

Equation (2.7) provides a well-defined description of co-volatility. We let (2.5) and (2.7) be heterogeneous by assuming random elements ${\omega }_{\mathit{ii}}$ and ${\omega }_{\mathit{ij}}$ as discussed later. This strategy extends the previous work on volatility modelling from the univariate to multivariate settings, in which a slowly moving unconditional volatility is incorporated into a GARCH-type model. The unconditional volatility directly depends on the macroeconomic information in a GARCH model. This article further studies the low-frequency co-volatility in the correlated cross-sectional individuals. Adopting random-effects models, the correlation between the $i$ th and $j$ th log-return is

$$Corr({ϵ}_{it},{ϵ}_{jt})=\frac{{\omega }_{ij}}{\sqrt{{\omega }_{ii}{\omega }_{jj}}}\frac{{\tau }_{ij,t}}{\sqrt{{\tau }_{ii,t}{\tau }_{jj,t}}}\frac{1-\gamma -\delta }{1-\lambda \frac{{\tau }_{ij,t}}{\sqrt{{\tau }_{ii,t}{\tau }_{jj,t}}}-\theta },$$

(2.8)

for $i,j=1,\dots ,N$ , $i\ne j$ . Equation (2.8) illustrates that the correlation coefficient is proportional to associated random effects and low-frequency factors. A new presentation of the originally proposed model is necessary to provide inference. Denote the $N\times 1$ vector ${\text{r}}_t=(r_{1t,}{r}_{2t},\dots ,r_{\mathit{Nt}}{)}^{\prime }$ for $t=1,\dots ,T$ . The N-dimensional unexpected returns vector ${\mathbf{\epsilon }}_t$ is conditionally normal with zero-mean and the $N\times N$ covariance matrix ${\mathrm{\Sigma }}_t$ having diagonal and off-diagonal elements (2.4) and (2.6), respectively.

The representation of unexpected returns in the vector notation is similar to the familiar MGARCH model, while random effects are present in the model specification. To estimate model parameters in frequentist or Bayesian frameworks, we may face two serious challenges: (a) there is no guarantee that the conditional covariance matrix of innovations, ${\mathrm{\Sigma }}_t$ for each $t=1,\dots ,T$ , is PD and (b) there is a large number of parameters which needs to be estimated. The following prerequisites are necessary before we present conditions to guarantee that ${\mathrm{\Sigma }}_t$ is PD, the following prerequisites are necessary. The covariance matrix ${\mathrm{\Sigma }}_t$ can easily be rewritten as

$$\begin{array}{cc}{\mathrm{\Sigma }}_t={\mathbf{\tau }}_t\odot {\text{H}}_t,\hfill & \hfill \end{array}$$

(2.9)

where $\odot$ denotes the Hadamard product of two matrices, ${\mathbf{\tau }}_t$ is the low-frequency volatility matrix with elements ${\tau }_{\mathit{ij},t}$ , and ${\text{H}}_t=\left\{h_{\mathit{ij},t}\right\}$ . If ${\mathbf{\tau }}_t$ and ${\text{H}}_t$ become positive semi-definite matrices, then so is ${\mathrm{\Sigma }}_t$ . The model can be represented in the Vech Matrix-Diagonal form to guarantee the positive definiteness of the conditional covariance matrix. Rewrite ${\text{H}}_t$ in the form

$$\begin{array}{cc}{\text{H}}_t={\text{CC}}^{\prime }+{\text{AA}}^{\prime }\odot ({\mathbf{ϵ}}_{t-1}{\mathbf{ϵ}}_{t-1}^{\prime }ø{\mathbf{\tau }}_{t-1}{\mathbf{\tau }}_{t-1}^{\prime })+{\text{BB}}^{\prime }\odot {\text{H}}_{t-1},\hfill & \hfill \end{array}$$

(2.10)

where $ø$ denotes Hadamard division, and $\text{C}$ , $\text{A}$ and $\text{B}$ denote $N\times N$ parameter matrices with either symmetric or the lower triangular form. Let ${\text{CC}}^{\prime }=\text{W}=\left\{{\omega }_{\mathit{ij}}\right\}$ ,

$${\text{AA}}^{\prime }=(\gamma -\lambda ){\text{I}}_N+\lambda {\text{J}}_N,\hspace{0.17em}\hspace{0.17em}{\text{BB}}^{\prime }=(\delta -\theta ){\text{I}}_N+\theta {\text{J}}_N,$$

(2.11)

where ${\text{I}}_N$ denotes the $N\times N$ identity matrix and ${\text{J}}_N=1_N{1}_N^{\prime }$ is a square unit matrix with $1_N=(1,\hspace{0.17em}1,\dots ,1{)}^{\prime }$ . The general structure for matrices ${\text{AA}}^{\prime }$ , ${\text{BB}}^{\prime }$ , and ${\text{CC}}^{\prime }$ imposes $3N(N+1)/2$ non-redundant parameters to the model. A simple Diagonal-Vech representation, proposed in the literature, is to restrict these coefficient matrices to diagonal. Also, the number of parameters is dramatically reduced because each conditional covariance depends only on its past values and innovations, and no interaction is included between conditional variances and covariances (Bollerslev et al.(1988)). While several structures are available, we examine compound-symmetry representation for autoregressive coefficients instead. It is a reasonable structure in the variance-covariance analysis, which incorporates homogeneous effects via $\gamma$ and $\delta$ for the variance equation and $\lambda$ and $\theta$ for the covariance equation. Also, we let the covariance matrix be heterogeneous by setting the elements of $\text{W}$ to follow a random-intercept model. This strategy sensibly reduces the number of coefficients, and heterogeneity in $\text{W}$ justifies variability. As seen later, it helps us to handle computational difficulties in empirical studies. Operating preceding conditions for all parameters, the matrix ${\mathrm{\Sigma }}_t$ will be PD for all $t$ if the initial covariance matrix ${\mathrm{\Sigma }}_1$ is PD.

Now, we present a practical trick to confirm the positive definiteness of the low-frequency matrix. Assume that

$$\begin{array}{cc}{\mathbf{\tau }}_t={\text{Q}}_t{\text{F}}_t^2{\text{Q}}_t^{\prime },\hfill & \hfill \end{array}$$

(2.12)

for $t=1,\dots ,T$ , where ${\text{F}}_t$ denotes a diagonal matrix of positive structural volatility elements and ${\text{Q}}_t$ is the lower triangular matrix of co-volatility components with unit diagonal elements and other elements determined by P-spline functions. Equation (2.12) leads to the modified Cholesky decomposition (Pourahmadi (2007)) for modelling low-frequency volatility.

Denote ${\text{F}}_t=\mathit{diag}(f_1(t),\dots ,f_N(t))$ and ${\text{g}}_t=(g_1(t),\dots ,g_P(t){)}^{\prime }$ as a vector of order $P=N\times (N-1)/2$ , which contains the lower-triangular elements of ${\text{Q}}_t=\left\{q_{\mathit{ij},t}\right\}$ stacked by its rows, where $q_{\mathit{ij},t}=g_p(t)$ for $p=i(i-1)/2-(i-j-1)$ . This novel formulation reduces the multivariate dynamics to a univariate volatility process and guarantees the PD of ${\mathbf{\tau }}_t$ .

We use a penalized splines technique, so-called P-splines (Ruppert et al.(2003); Crainiceanu et al.(2005)). Let two functions $f_i$ and $g_j$ be deterministic low-rank regression spline functions. These functions can be derived by dividing the domain of each time point $t$ into continuous intervals and representing the unknown functions $f_i$ and $g_j$ by a separate function in each interval. Also, each function should join smoothly at the fixed knots ${\kappa }_k$ for $k=1,\dots ,K$ and ${\varkappa }_s$ for $s=1,\dots ,S$ . Thus, the low-rank thin-plate splines must be continuous on the whole interval $[1,T]$ . For fixed knots, we have

$$f_i\left(t\mid {\left\{{\kappa }_k\right\}}_{k=1}^K\right)={\varphi }_{i0}+{\varphi }_{i1}t+\sum _{k=1}^K{v}_{ik}{\left|t-{\kappa }_k\right|}^3,$$

(2.13)

$$g_j\left(t\mid {\left\{{\varkappa }_s\right\}}_{s=1}^S\right)={\psi }_{j_0}+{\psi }_{j1}t+\sum _{s=1}^S{u}_{js}{\left|t-{\varkappa }_s\right|}^3,$$

(2.14)

where the regression coefficients ${\text{V}}_k=(v_{1k},\dots ,v_{\mathit{Nk}}{)}^{\prime }$ and ${\text{U}}_s=(u_{1s},\dots ,u_{\mathit{Ps}}{)}^{\prime }$ represent the change in gradient between consecutive line segments. Smooth estimates of $f_i$ and $g_j$ may be obtained by allowing $v_{\mathit{ik}}$ and $u_{\mathit{js}}$ be random coefficients, independently distributed as $v_{ik}\stackrel{iid}{\sim }N\left(0,{\sigma }_v^2\right)$ and $u_{js}\stackrel{iid}{\sim }N\left(0,{\sigma }_u^2\right)$ (Gurrin et al.(2005)). Furthermore, $\{{\kappa }_1=0<{\kappa }_2<\dots <{\kappa }_K=T\}$ and $\{{\varkappa }_1=0<{\varkappa }_2<\dots <{\varkappa }_s=T\}$ , respectively, represent partitions of the time horizon $T$ into $K$ and $S$ equally spaced intervals. Values $K$ and $S$ denote the number of knots in the spline mixed-effects model of associated variances and covariances. The analyst can employ some information criteria to determine $K$ and $S$ optimal choices, which govern the cyclical pattern in the low-frequency volatility trend. Large values of $K$ and $S$ provide more frequent cycles. Coefficient vectors govern the sharpness of each cycle ${\mathbf{\varphi }}_i=({\varphi }_{i0},{\varphi }_{i1}{)}^{\prime }$ and ${\mathbf{\psi }}_j=({\psi }_{j0},{\psi }_{j1}{)}^{\prime }$ . Following Ruppert (2002), one can consider the number of knots large enough to ensure the desired flexibility, where $K$ and $S$ are the sample quantiles of time $t$ corresponding to probabilities $k/(K+1)$ and $s/(S+1)$ for all $k$ and $s$ . However, other choices of knots are also available. For (2.13), the basis function set is $\{1,t,|t-{\kappa }_1{|}^3,\dots ,|t-{\kappa }_K{|}^3\}$ , since any low-rank thin-plate splines with the given knots is a linear combination of this set, while for (2.14) the basis set includes $\left\{1,t,{\left|t-{\varkappa }_1\right|}^3,\dots ,{\left|t-{\varkappa }_S\right|}^3\right\}$ .

Smoothing splines require a certain background to use suitable spline functions. In practical applications, regression splines yield similar estimates with fewer knots than smoothing splines, particularly when estimated by discrete penalized techniques. Various splines can be formed by altering the choice of knots and changing how roughness in the estimated regression function is penalized. An approach allows a knot at each value of the variable, which, given an appropriate choice of roughness penalty, leads to a natural cubic smoothing spline. It is necessary for fitting regression splines to carefully choose knot locations and basis functions to have subjectivity in the model fitting process. A greater smoothness can be achieved by shrinking the estimated coefficients towards zero.

Let ${\text{X}}_t=\{1,t\}$ , ${\mathbf{Z}}_{1t}=\left\{\right|t-{\kappa }_1{|}^3,\dots ,|t-{\kappa }_K{|}^3\}$ , and $Z_{2t}=\left\{\right|t-{\varkappa }_1{|}^3,\dots ,|t-{\varkappa }_S{|}^3\}$ . Non-parametric regression equations $f_i$ and $g_j$ can be represented by the mixed-effects models

$$f_i(t)={\text{X}}_t{\varphi }_i+{\text{Z}}_{1t}{\text{V}}_k,$$

(2.15)

$$g_j(t)={\text{X}}_t{\psi }_j+{\text{Z}}_{2t}{\text{U}}_s,$$

(2.16)

for $i=1,\dots ,N$ and $j=1,\dots ,P$ . The first terms in (2.15) and (2.16) specify fixed effects, and the second terms smooth unknown functions $f_i$ and $g_j$ , respectively. This modelling strategy imposes estimating a large number of unknown parameters, especially when the number of assets increases. It may render low-frequency volatility models to be unidentifiable. To handle these issues in empirical studies, we consider assumptions for ${\mathbf{\varphi }}_i$ and ${\mathbf{\psi }}_j$ . In modelling $f_i$ and $g_j$ , we let the following forms of parameter heterogeneity:

Homogeneous intercepts and slopes: Coefficient vectors ${\mathbf{\varphi }}_i$ and ${\mathbf{\psi }}_j$ reduce to $\mathbf{\varphi }$ and $\mathbf{\psi }$ , restricting common intercepts and slopes for assets.

Heterogeneous intercepts and homogeneous slopes: Coefficient vectors ${\mathbf{\varphi }}_i$ and ${\mathbf{\psi }}_j$ reduce to $({\varphi }_{i0},{\varphi }_1{)}^{\prime }$ and $({\psi }_{j0},{\psi }_1{)}^{\prime }$ , such that ${\varphi }_{i0}$ and ${\psi }_{j0}$ for $i=1,\dots ,N$ and $j=1,\dots ,P$ are individual random effects that follow known distributions $G_{\mathrm{\varphi }}(·)$ and $G_{\mathrm{\psi }}(·)$ , respectively. The low-frequency volatility model for variance and covariance components includes common slopes.

Both intercepts and slopes heterogeneous: Let coefficient vectors ${\mathbf{\varphi }}_i$ and ${\mathbf{\psi }}_j$ for $i=1,\dots ,N$ and $j=1,\dots ,P$ be random effects that follow distributions $G_{\mathrm{\varphi }}(·)$ and $G_{\mathrm{\psi }}(·)$ .

Since the proposed model contains many latent variables, it complicates the estimation process and enforces advanced numerical approaches. The Hamiltonian Monte Carlo is a useful Bayesian computing tool to deal with such issues.

3 Bayesian computational aspects

Bayesian computational tools are convenient to fit complex financial models that appear in hierarchical representations. Bayesian inference requires the joint posterior distribution of all unknown model parameters. The conditional density of ${\mathbf{ϵ}}_t$ is given by

$$\begin{array}{cc}\pi \left({\mathbf{\epsilon }}_t\right|{\mathrm{\Sigma }}_t)=(2\pi {)}^{-\frac{N}{2}}{\left|{\mathrm{\Sigma }}_t\right|}^{-\frac{1}{2}}\text{exp}(-\frac{1}{2}{\mathbf{\epsilon }}_t{}^{\prime }{\mathrm{\Sigma }}_t^{-1}{\mathbf{\epsilon }}_t).\hfill & \hfill \end{array}$$

(3.1)

For the complete dataset, the conditional likelihood function is $\pi \left(\mathbf{ϵ}\right|{\mathrm{\Sigma }}_1,\dots ,{\mathrm{\Sigma }}_T)={\prod }_{t=1}^T\pi \left({\mathbf{ϵ}}_t\right|{\mathrm{\Sigma }}_t),$ where matrix $\mathbf{ϵ}=({\mathbf{ϵ}}_1,\dots ,{\mathbf{ϵ}}_T)$ . As (2.5) and (2.7) involve some restrictions, the Bayesian implementation requires assigning priors of bounded supports.

3.1 Prior specification for the high-frequency volatility model

The Bayesian analysis necessities to determine prior distributions for parameters under conditions ${\omega }_{\mathit{ii}}>0$ , $\gamma +\delta <1$ , and $\lambda +\theta <1$ . These are sufficient for the conditional variance and covariance processes to converge to some fixed numbers or a certain fixed PD matrix. Furthermore, the PD of ${\text{AA}}^{\prime }$ , ${\text{BB}}^{\prime }$ , and ${\text{CC}}^{\prime }$ in (2.11) is sufficient for ${\text{H}}_t$ ’s positive semi-definiteness. We notice $\gamma -\lambda \ge 0$ and $\delta -\theta \ge 0$ guarantee the positive semi-definiteness of ${\text{AA}}^{\prime }$ and ${\text{BB}}^{\prime }$ , respectively. Thus, algebraic operations directly assign joint priors with bounded support $0<\lambda \le \gamma <1-\delta \le 1-\theta$ or $0<\theta \le \delta <1-\gamma \le 1-\lambda$ . However, such priors are inapplicable unless we transform constraints into simple inequalities. As detailed later in the data analysis section, we recommend that each parameter be uniformly distributed with restriction to the associated stationary regions.

To study the impact of heterogeneity on the variance and covariance equations, suppose that random effects matrix, $\text{W}=\left\{{\omega }_{\mathit{ij}}\right\}$ , follows the Wishart distribution with $N\times N$ scale matrix ${\text{S}}_{\omega }$ and a known degrees-of-freedom $N$ . This assumption confirms the PD of ${\text{H}}_t$ in (2.10). Moreover, we adopt a compound symmetry structure for ${\text{S}}_{\omega }$ , with diagonals ${\sigma }_a^2+{\sigma }_b^2$ and off-diagonal elements ${\sigma }_b^2$ , such that ${\sigma }_a$ and ${\sigma }_b$ follow half-Cauchy $(0,5)$ , see Gelman (2006).

3.2 Prior specification for the low-frequency volatility model

We adopt $v_{ik}\stackrel{iid}{~}N\left(0,{\sigma }_v^2\right)$ , $k=1,\dots ,K$ , $i=1,\dots ,N$ , and ${\sigma }_v\sim$ half-Cauchy $(0,5)$ , similarly $u_{js}\stackrel{iid}{~}N\left(0,{\sigma }_u^2\right)$ for $s=1,\dots ,S$ , $j=1,\dots ,P$ , and ${\sigma }_u\sim$ half-Cauchy $(0,5)$ , as proposed by Gelman (2006) and Polson and Scott (2012). To allow intercepts heterogeneity and slopes homogeneity, let each coefficient of the low-rank spline functions follow independently normal distribution, that is, ${\varphi }_{i0}\stackrel{iid}{~}N\left(0,{\sigma }_{\varphi }^2\right)$ , ${\sigma }_{\varphi }\sim$ half-Cauchy $(0,5)$ , ${\varphi }_1\sim N(0,1)$ , ${\psi }_{j0}\stackrel{iid}{~}N\left(0,{\sigma }_{\psi }^2\right)$ , ${\sigma }_{\psi }\sim$ half-Cauchy $(0,2.5)$ and ${\psi }_1\sim N(0,1)$ (Ruppert et al.(2003); Durban and Currie (2003); Zhang et al.(1998)).

3.3 Parameter estimation via Hamiltonian Monte Carlo

Hamiltonian Monte Carlo is an effective computational approach (Duane et al.(1987); Neal (1994)) when dealing with complex modelling problems. It uses an approximate Hamiltonian dynamics simulation based on derivatives of density functions to generate an efficient sample from posterior density functions. Neal (2011) and Betancourt and Girolami (2015) detailed some theoretical and practical topics.

Let $\mathbf{\vartheta }$ be the vector of all parameters in (2.9) and $\mathbf{\varphi }$ be a vector of auxiliary momentum variables, independent of $\mathbf{\vartheta }$ , which jointly drawn from $\pi (\mathbf{\vartheta },\mathbf{\varphi })=\pi \left(\mathbf{\varphi }\right|\mathbf{\vartheta })$ $\pi (\mathbf{\vartheta })$ . The $\mathbf{\varphi }$ is commonly assumed to follow $N(0,{\mathrm{\Sigma }}_{\varphi })$ . The covariance matrix ${\mathrm{\Sigma }}_{\varphi }$ plays the role of Euclidean metric to rotate and scale the target density, which may be set to the identity matrix, estimated from warm-up draws, or restricted to a diagonal matrix ${\sigma }_{\varphi }^2\text{I}$ (Betancourt and Stein (2011)). The total energy or Hamiltonian function is

$$\begin{array}{ccc}\mathcal{H}(\mathbf{\vartheta },\mathbf{\varphi })\hfill & =\hfill & -\text{log}\pi (\mathbf{\vartheta },\mathbf{\varphi })\hfill \\ \hfill & =\hfill & -\text{log}\pi \left(\mathbf{\varphi }\right|\mathbf{\vartheta })-\text{log}\pi (\mathbf{\vartheta })\mathrm{notag}\hfill \\ \hfill & =\hfill & \mathcal{F}\left(\mathbf{\varphi }\right|\mathbf{\vartheta })+\mathcal{L}(\mathbf{\vartheta }),\mathrm{notag}\hfill \end{array}$$

(3.2)

where $\mathcal{F}\left(\mathbf{\varphi }\right|\mathbf{\vartheta })=-\text{log}\pi \left(\mathbf{\varphi }\right|\mathbf{\vartheta })$ denotes the kinetic energy and $\mathcal{L}(\mathbf{\vartheta })=-\text{log}\pi (\mathbf{\vartheta })$ is called potential energy. For the independent case $\pi \left(\mathbf{\varphi }\right|\mathbf{\vartheta })=\pi (\mathbf{\varphi })$ , the kinetic energy equals $\mathbf{\varphi }{\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\Sigma }}_{\varphi }^{-1}$ ${\mathbf{\varphi }}^{\prime }$ and the joint density reduces to

$$\pi (\vartheta ,\varphi )\propto \pi (\vartheta )\times \text{exp}(-\phi {\text{Σ}}_{\phi }^{-1}{\phi }^{\prime })\propto \text{exp}(-\mathcal{H}(\vartheta ,\phi )).$$

For the continuous time $t$ , the new momentum is evolved via the following Hamilton equations

$$\frac{\partial \vartheta }{\partial t}=\frac{\partial \mathcal{H}(\vartheta ,\varphi )}{\partial \varphi }=\frac{\partial \mathcal{F}\text{(}\varphi \text{|}\vartheta )}{\partial \varphi }={\Sigma }_{\varphi }^{-1}\varphi ,$$

(3.3)

$$\frac{\partial \varphi }{\partial t}=-\frac{\partial \mathcal{H}(\vartheta ,\varphi )}{\partial \vartheta }=-\frac{\partial \mathcal{L}(\vartheta )}{\partial \vartheta }=-{\nabla }_{\vartheta }\mathcal{L}(\vartheta ),$$

(3.4)

where ${\nabla }_{\mathbf{\vartheta }}\mathcal{L}(\mathbf{\vartheta })$ is the gradient of $\mathcal{L}(\mathbf{\vartheta })$ for $\mathbf{\vartheta }$ . Performing the simulation process demands a vector of auxiliary momentum variables. Then, the first derivative computation of the density function proceeds for unknown quantities to generate an efficient exploration of the posterior distribution. Advanced numerical algorithms are needed to solve (3.3) and (3.3). The HMC benefits from a numerical integration algorithm, so-called the leapfrog or Stormer-Verlet integrator, to produce reliable solutions of Hamiltonian equations (Leimkuhler and Reich (2004)). The leapfrog algorithm discretizes a short time interval of $\epsilon$ and starts drawing a new momentum sample of the parameter or the previous momentum value. The Hamiltonian dynamics alternates half-step updates of the momentum and full-step updates of the position as the following steps

$$\begin{array}{ccc}\hfill \mathbf{\phi }& \hfill \leftarrow \hfill & \mathbf{\phi }-\frac{\epsilon }{2}{\nabla }_{\mathbf{\vartheta }}\mathcal{L}(\mathbf{\vartheta }),\hfill \\ \hfill \mathbf{\vartheta }& \hfill \leftarrow \hfill & \mathbf{\vartheta }+\epsilon {\mathbf{\Sigma }}_{\phi }^{-1}\mathbf{\phi },\hfill \\ \hfill \mathbf{\phi }& \hfill \leftarrow \hfill & \mathbf{\phi }-\frac{\epsilon }{2}{\nabla }_{\mathbf{\vartheta }}\mathcal{L}(\mathbf{\vartheta }),\hfill \end{array}$$

for some small step-size $\epsilon >0$ . These steps repeat $m$ times, and simulation results $({\mathbf{\vartheta }}^{*},{\mathbf{\varphi }}^{*})$ are obtained from a total of $m\times \epsilon$ simulations. In practice, the Metropolis-acceptance step keeps the proposed $({\mathbf{\vartheta }}^{*},{\mathbf{\varphi }}^{*})$ with probability

$$\text{min}[1,\text{exp}(\mathcal{H}(\mathbf{\vartheta },\mathbf{\phi })-\mathcal{H}({\mathbf{\vartheta }}^{*},{\mathbf{\phi }}^{*})],$$

otherwise, the previous parameter value is returned for the next draw and used as an initial value for the next iteration. A summary of the HMC algorithm is as follows

initialize parameters, set initial values ${\mathbf{\vartheta }}^{(0)}$ and let $l=1$ ,

generate ${\mathbf{\phi }}^{*}\sim N(0,\text{I})$ and $u\sim U(0,1),$

set $(\mathbf{\vartheta },\mathbf{\phi })=({\mathbf{\vartheta }}^{(l-1)},{\mathbf{\phi }}^{*})$ and $\mathcal{H}=\mathcal{H}(\mathbf{\vartheta },\mathbf{\phi }),$

repeat the leapfrog solution $m$ times

${\mathbf{\phi }}^{*}={\mathbf{\phi }}^{*}-\frac{\epsilon }{2}{\nabla }_{\mathbf{\vartheta }}$ $\mathcal{L}({\mathbf{\vartheta }}^{(l-1)}),$

${\mathbf{\vartheta }}^{(l-1)}={\mathbf{\vartheta }}^{(l-1)}+\epsilon$ ${\mathbf{\phi }}^{*},$

${\mathbf{\phi }}^{*}={\mathbf{\phi }}^{*}-\frac{\epsilon }{2}{\nabla }_{\mathbf{\vartheta }}$ $\mathcal{L}({\mathbf{\vartheta }}^{(l-1)}),$

set $({\mathbf{\vartheta }}^{*},{\mathbf{\phi }}^{*})=({\mathbf{\vartheta }}^{(l-1)},{\mathbf{\phi }}^{*})$ and ${\mathcal{H}}^{*}=\mathcal{H}({\mathbf{\vartheta }}^{*},{\mathbf{\phi }}^{*}),$

compute $\rho =\text{min}[1,\text{exp}(\mathcal{H}-{\mathcal{H}}^{*})],$

set ${\mathbf{\vartheta }}^{(l)}={\mathbf{\vartheta }}^{*}$ if $\rho >u$ , otherwise set ${\mathbf{\vartheta }}^{(l)}=\mathbf{\vartheta },$

let $l=l+1$ and return to step 2 until convergence.

To perform data analysis of the empirical studies, we employ Stan’s probabilistic programming language, which uses the HMC algorithm. A Stan code is initially compiled to a C++ programme by the facility of Stan compiler, the so-called stanc. Next, a self-contained platform-specific executable recompiles the code. We use the RStan library, the R interface to Stan, see Carpenter et al.(2017), and ‘GitHub, https://github.com/stan-dev/rstan/wiki/RStan-Getting-Started’. RStan shares initialization, sampling, tuning controls and approximation to provide the posterior analysis.

3.4 Evaluating the predictive accuracy

The posterior predictive distribution (PPD) has often been used to check model performance. For $N$ assets and total time $T$ , let $\mathbf{\vartheta }$ denote the parameter vector and $\mathbf{\epsilon }$ the observed log-return. The posterior density of $\mathbf{\vartheta }$ , given $\mathbf{ϵ}$ , measures the uncertainty of estimates. Let ${\mathbf{ϵ}}_{\mathit{new}}$ be a vector of unobserved log-returns at a forthcoming time. The PPD is computed by

$$\begin{array}{c}\pi \left({ϵ}_{new}\mid ϵ\right)={{\displaystyle \int }}_{\vartheta }\pi \left({ϵ}_{new}\mid \vartheta ,ϵ\right)\pi (\vartheta \mid ϵ)d\vartheta \\ \approx \sum _{l=1}^L{L}^{-1}\pi \left({ϵ}_{new}\mid {\vartheta }^{(l)},ϵ\right)\end{array}$$

where $L$ is the number of essential Monte Carlo samples to reach convergence (Gelman et al.(2013)).

3.5 Optimization-based inference

Denoting $\pi \left(\mathbf{\vartheta }\right|\mathbf{\epsilon })$ the posterior density, the maximum a posteriori (MAP) estimate is derived by

$$\begin{array}{ccc}\hfill {\widehat{\mathbf{\vartheta }}}_{\mathit{MAP}}& \hfill =\hfill & \text{arg max}{}_{\mathbf{\vartheta }}(\pi \left(\mathbf{\vartheta }\right|\mathbf{ϵ}))\hfill \\ \hfill & \hfill =\hfill & \text{arg max}{}_{\mathbf{\vartheta }}(\pi \left(\mathbf{ϵ}\right|{\mathrm{\Sigma }}_1,\dots ,{\mathrm{\Sigma }}_T)\times \pi (\mathbf{\vartheta })),\hfill \end{array}$$

where $\text{arg}{\text{max}}_{\vartheta }(·)$ determines the value of $\mathbf{\vartheta }$ , maximizing the desired function. The maximum likelihood estimation (MLE) and the MAP estimation are closely related. By setting a flat prior over entire values of $\mathbf{\vartheta }$ , the posterior mode corresponds to the MLE of parameters. In some situations, such as predicting volatility in stock markets, the MAP estimation may produce more accurate results than the MLE. The MAP estimation procedure uses previous information about parameter values of a particular form. Technical details are given by Bassett and Deride (2019).

4 Simulation and empirical application

We first conduct simulation studies to investigate the statistical properties of the time-varying GARCH mixed-effects model and then analyse the stock market dataset to model volatility and co-volatility.

4.1 Simulation experiments

We simulate the unexpected returns, conditioned on latent variables and innovation terms, and generate some proposed time-varying multivariate model paths. Financial time-series data exhibit specific features that should be included in the simulation study. Specifically, the simulated series variability is taken into account by time-varying coefficients. We expect that non-stationary processes with high variation intend extreme peaks. Also, variance components control volatility and co-volatility in the sample paths.

We first conduct a preliminary simulation study to investigate the Bayesian sensitivity analysis. Particularly, to evaluate the priors effect and hyperparameters choice on the posterior distribution, we examine a series of prior distributions, especially for variance components introduced by Gelman (2006), and study convergence by generating parallel chains. We also defined flat priors for all parameters respecting the already mentioned bounded supports. Simulation findings revealed that priors in Section 3 and flat priors produced almost the same estimation results, implying the model fitting was less sensitive to the choice of priors. For fixed hyperparameters in low-frequency equations, results of the simulated experiments were obtained by computing a series of means and mean-squared errors (MSEs) for all parameters of high-frequency equations. We organize several plans for simulations to investigate various properties of our proposed model:

Simulation 1: When the number of individuals varies, the model behaviour may be interesting. We designed some simulation scenarios, comprising dimensions $N=10$ , $15$ , $20$ and $40$ over length $T=180$ . After the successful convergence, Table 1 reports our findings, including the parameters estimate and the latent effects with their MSEs. The bias monotonously decreases as $N$ gets large, that is, it decreases or is nearly equal to the true values. Hence, results verify the consistency of estimates for various $N$ .

Table 1

Simulation results for the proposed model with different $N$ and $T=180$

		$\gamma$	$\lambda$	$\delta$	$\theta$	${\sigma }_a$	${\sigma }_b$
	True	0.2500	0.1000	0.4000	0.2000	0.1000	0.1000
$N=10$	Mean	0.2552	0.1410	0.3692	0.1823	0.2238	0.0552
	MSE	0.0012	0.0044	0.0046	0.0109	0.0159	0.0033
$N=15$	Mean	0.2734	0.1138	0.4008	0.1968	0.1012	0.0901
	MSE	0.0013	0.0012	0.0023	0.0095	0.0001	0.0005
$N=20$	Mean	0.2239	0.1017	0.3670	0.1249	0.1123	0.1094
	MSE	0.0012	0.0005	0.0036	0.0122	0.0002	0.0005
$N=40$	Mean	0.2524	0.0955	0.3859	0.1665	0.0891	0.1403
	MSE	0.0002	0.0002	0.0012	0.0120	0.0002	0.0019

Simulation 2: The model assessment was performed for four simulation experiments for $N=10$ and $T=180$ . For all series in Table 2, we observed that the estimates were close to their true values with MSEs near to zero. Thus, results confirm that our proposed model performs satisfactorily in terms of unbiasedness and precision.

Table 2

Simulation results for the proposed model with $N=10$ and $T=180$

		$\gamma$	$\lambda$	$\delta$	$\theta$	${\sigma }_a$	${\sigma }_b$
Series I	True	0.5000	0.2000	0.4000	0.1000	0.2000	0.2000
	Mean	0.4893	0.2252	0.4152	0.1268	0.2134	0.2511
	MSE	0.0019	0.0043	0.0013	0.0074	0.0012	0.0086
Series II	True	0.2500	0.1000	0.2500	0.1000	0.1000	0.1000
	Mean	0.2351	0.1382	0.2655	0.1308	0.1038	0.0913
	MSE	0.0015	0.0041	0.0078	0.0089	0.0002	0.0007
Series III	True	0.2500	0.1000	0.4000	0.2000	0.1000	0.1000
	Mean	0.2552	0.1410	0.3692	0.1823	0.2238	0.0552
	MSE	0.0012	0.0044	0.0046	0.0109	0.0159	0.0033
Series IV	True	0.5000	0.1000	0.2500	0.2000	0.2000	0.3000
	Mean	0.5751	0.0932	0.2549	0.1755	0.1607	0.2887
	MSE	0.0079	0.0022	0.0017	0.0051	0.0026	0.0056

Simulation 3: To compare the proposed model with some existing multivariate GARCH models, we conduct an additional simulation study. We first generate $N=10$ vectors of unexpected returns from the time-varying GARCH mixed-effects model ${\mathbf{ϵ}}_t={\mathrm{\Sigma }}_t^{1/2}{\mathit{e}}_t$ with $T=180$ , and the covariance matrix ${\mathrm{\Sigma }}_t$ is given by (2.9) and the innovation term ${\mathit{e}}_t$ follows the contaminated-normal distribution ${\mathit{CN}}_N(0,\text{I},c_1,c_2)$ , where $c_1$ and $c_2$ respectively determine the contamination probability and the degree of contamination (Arellano-Valle et al.(2018)). Set $\gamma =0.25$ , $\lambda =0.1$ , $\delta =0.25$ , $\theta =0.1$ , ${\sigma }_{\varphi }=0.1$ , ${\sigma }_{\psi }=0.25$ , ${\sigma }_u=0.25$ , ${\sigma }_v=0.1$ , ${\sigma }_a=0.1$ , ${\sigma }_b=0.1$ , ${\varphi }_1=0.1$ and ${\psi }_1=0.1$ . We also assume $c_1\sim \mathit{Uniform}(0,1)$ and $c_2\sim \mathit{Uniform}(1,5)$ . Next, we fitted the following competing multivariate models to the generated samples:

M1: The proposed model in Section 2.

M2: The proposed model with eliminating co-volatilities by assuming ${\mathrm{\Sigma }}_t=\mathit{diag}({\tau }_{\mathit{ii},t}{h}_{\mathit{ii},t})$ , where $h_{\mathit{ii},t}$ is defined in (2.5) such that ${\omega }_{\mathit{ii}}\sim \mathit{Gamma}(s/3,s/3)$ , $s\sim \mathit{Inv}-\mathit{Gamma}(3,3)$ , and ${\tau }_{\mathit{ii},t}=\mathit{exp}\left(f_i(t)\right)$ , with $f_i(t)$ in (2.14).

M3: An extension of the univariate spline GARCH model proposed by Engle and Rangel (2008) to the multivariate case; for which, ${\mathrm{\Sigma }}_t=\mathit{diag}({\tau }_{\mathit{ii},t}{h}_{\mathit{ii},t})$ , where $h_{\mathit{ii},t}=(1-\gamma -\delta )+\gamma \left({ϵ}_{i,t-1}^2/{\tau }_{\mathit{ii},t-1}\right)+\delta h_{\mathit{ii},t-1}$ , and ${\tau }_{\mathit{ii},t}=\text{exp}({\varphi }_{i0}+{\varphi }_{i1}t+{\sum }_{k=1}^K{v}_{\mathit{ik}}(t-{\kappa }_k{)}_{+}^2$ , the sharpness of each cycle is govern by fixing coefficients ${\varphi }_{i0}$ and ${\varphi }_{i1}$ with defining the basis function $(t-{\kappa }_k{)}_{+}=(t-{\kappa }_k)$ if $t>{\kappa }_k$ and $0$ otherwise. ${\kappa }_k$ is the sample quantiles of time $t$ corresponded to probabilities $k/(K+1)$ for $k=1,\dots ,K$ . We set $K=7$ . The difference between models M2 and M3 is that for M3, variance and covariance equations are homogeneous, and quadratic spline basis functions are used.

M4: The Matrix-Diagonal model proposed by Ding and Engle (2001), with ${\mathrm{\Sigma }}_t={\text{CC}}^{\prime }+{\text{AA}}^{\prime }\odot {\mathbf{ϵ}}_{t-1}{\mathbf{ϵ}}_{t-1}^{\prime }+{\text{BB}}^{\prime }\odot {\mathrm{\Sigma }}_{t-1}$ , where ${\text{CC}}^{\prime }$ , ${\text{AA}}^{\prime }$ and ${\text{BB}}^{\prime }$ are all positive semi-definite, implying the PD of ${\mathrm{\Sigma }}_t$ , for all $t>1$ provided that the initial covariance matrix ${\mathrm{\Sigma }}_1$ being PD. We adopt ${\mathrm{\Sigma }}_t$ to be compound symmetry and set ${\text{CC}}^{\prime }=({\omega }_1-{\omega }_2){\text{I}}_N+{\omega }_2{\text{J}}_N$ , ${\text{AA}}^{\prime }=(\gamma -\lambda ){\text{I}}_N+\lambda {\text{J}}_N$ and ${\text{BB}}^{\prime }=(\delta -\theta ){\text{I}}_N+\theta {\text{J}}_N$ . Also, fixing ${\omega }_1$ , $\gamma$ and $\delta$ for the variance and ${\omega }_2$ , $\lambda$ and $\theta$ makes the covariance matrix homogeneous.

M5: The heterogeneous multivariate GARCH model (Cermeño and Grier (2006)). Heterogeneity in $\text{W}$ justifies variability and practicality. The model differs from M4 since the covariance matrix lets to be heterogeneous by setting ${\text{CC}}^{\prime }$ compound symmetry. In contrast, we assume ${\text{CC}}^{\prime }=\text{W}$ follows the Wishart distribution, which guarantees the ${\mathrm{\Sigma }}_t$ to be PD for all $t$ as the ${\mathrm{\Sigma }}_1$ being PD.

Existing MGARCH models have specific features, while some, such as M4 and M5, can not capture high- and low-frequency volatilities. The routine practice uses information criteria for models choice. First, let $\mathbf{\vartheta }$ be a $k\times 1$ vector of all parameters estimated in a fitted model and the deviance $D(\mathbf{\vartheta })=-2\text{log}(\pi (\mathbf{\vartheta }|\mathbf{ϵ}))$ . The most used criteria are (a) Akaike information criterion, $\mathit{AIC}(\mathbf{\vartheta })=D(\mathbf{\vartheta })+2k$ , (b) Bayesian information criterion, $\mathit{BIC}(\mathbf{\vartheta })=D(\mathbf{\vartheta })+\text{log}(T)k$ , (c) Deviance information criterion, $\mathit{DIC}(\mathbf{\vartheta })=-2\mathit{lpd}(\mathbf{\vartheta })+2k_{\mathit{DIC}}$ , where $\mathit{lpd}(\mathbf{\vartheta })$ denotes the log-predictive density of the data given a point estimate of the fitted model (Spiegelhalter et al.(2002); Gelman et al.(2014)), $\mathit{lpd}(\mathbf{\vartheta })=\text{log}\pi (ϵ|\mathbf{\vartheta })$ , and $k_{\mathit{DIC}}=2\{\text{log}\pi (ϵ|\mathbf{\vartheta })-E_{\mathit{post}}(\text{log}\pi (ϵ|\vartheta ))\}$ is the penalty function, (d) Watanabe-Akaike information criterion, $\mathit{WAIC}(\mathbf{\vartheta })=-2\mathit{lppd}(\mathbf{\vartheta })+2k_{\mathit{WAIC}}$ , where the log point-wise posterior predictive density $\mathit{lppd}(\mathbf{\vartheta })={\sum }_{t=1}^T\text{log}{\int }_{\mathbf{\vartheta }}\pi ({\mathbf{\varkappa }}_t|\mathbf{\vartheta },{\mathbf{ϵ}}_{t-1}){\pi }_{\mathit{post}}(\mathbf{\vartheta })d\mathbf{\vartheta }$ , and $k_{\mathit{WAIC}}=2{\sum }_{t=1}^T\{\text{log}{E}_{\mathit{post}}(\pi ({\mathbf{ϵ}}_t|\mathbf{\vartheta },{\mathbf{ϵ}}_{t-1}))-E_{\mathit{post}}(\text{log}\pi ({\mathbf{ϵ}}_t|\mathbf{\vartheta },{\mathbf{ϵ}}_{t-1})\}$ . WAIC estimate the out-of-sample prediction accuracy from a fitted Bayesian model using the log-likelihood evaluated at the posterior simulations of the parameter values (Watanabe and Opper (2010); Vehtari et al.(2017)).

Table 3

Information criteria for fitted models

	M1	M2	M3	M4	M5	True Model
AIC	$-323.80$	$-304.06$	$-256.45$	$-160.38$	$-233.01$	$-354.40$
BIC	$-320.49$	$-302.27$	$-250.32$	$-158.59$	$-231.22$	$-346.57$
DIC	$-1115.58$	$-1071.20$	$-1034.44$	$-628.68$	$-825.23$	$-1264.52$
WAIC	$-1169.45$	$-1063.92$	$-1023.60$	$-616.92$	$-812.24$	$-1293.70$

For each model, we ran Stan codes until successful convergence. Table 3 shows that the proposed model M1 performs satisfactorily compared to competing models since it has the least difference of criteria values with the true model. Notice that the data generating process was upon the contaminated normal distribution that involves outlying points. Thus, our results generally reveal that the fitted models involving cross-covariance, especially those accommodating more flexibility and variability through random effects, could deal with contaminated data, implying robustness to possible outlying points. The worst is M4 since it does not account for heterogeneity and volatility or co-volatility isolation.

4.2 A study on volatility and co-volatility in 10 top US stock indexes

Stock market data analysis has often been desirable in economic and financial studies. As emphasized, the long-term and short-term changes in a stock price cause significant effects on other stocks and macroeconomic factors. Modern regression volatility models, such as widely used GARCH, are available to investigate the impact of such effects. This section analyses the weekly log-returns of the 10 top stock indices from the Consumer Staples sector. Our study contains Walmart (WMT), Procter and Gamble (PG), Coca-Cola (KO), Philip Morris International (PM), Altria Group (MO), Estee Lauder (EL) and Colgate-Palmolive (CL) from the New York Stock Exchange (NYSE), and PepsiCo (PEP), Costco (COST) and Mondelez International (MDLZ) from the National Association of Securities Dealers Automated Quotations exchange (NASDAQ). The data cover 110 weeks from 6 November, 2017 to 2 December, 2019, and is freely available on https://finance.yahoo.com. We aim to model high- and low-frequency volatility for each stock index appropriately and study the effect of changing one stock price on the other stocks. Figure 1 presents the time-series plots of the log-returns of 10 stocks. None of the stock datasets follows the stationary process.

Figure 1

Time-series plots for the log-returns of 10 stocks

Now, we assess the proposed model performance in empirical applications compared to current models for cross-correlated financial data. It is instructive to signify model's ability to isolate high- and low-frequency volatilities. Estimation results, shown in Table 4, reveal that our model is the best fitting. Although the heterogeneous multivariate GARCH model fits fairly satisfactorily, it cannot isolate frequencies, unlike our model.

Table 4

Model comparison for fitted to stock market data

	M1	M2	M3	M4	M5
AIC	$-3035.67$	$-2822.87$	$-2777.78$	$-2935.41$	$-3023.49$
BIC	$-3035.19$	$-2822.61$	$-2776.88$	$-2935.14$	$-3023.23$
DIC	$-10873.65$	$-10161.52$	$-10141.59$	$-10633.13$	$-10830.17$
WAIC	$-10771.08$	$-10153.55$	$-10145.61$	$-9091.47$	$-10839.71$

As mentioned in Section 3.3, the RStan library in R utilizes the HMC simulation approach. Stan provides the mode of posterior or the MAP estimate. The Newton-Raphson algorithm is the default optimization of Stan (Carpenter et al.(2017)). Also, by adopting non-informative priors, the MAP estimates are equivalent to the MLEs. Table 5 presents the estimation results for the given period.

Table 5

Bayesian results, standard deviations, MLE and MAP estimates, and rhat values

Parameters	Mean	sd	MLE	MAP	rhat
$\gamma$	0.0995	0.0344	0.1473	0.1090	1.0018
$\lambda$	0.0220	0.0200	0.0269	0.0111	1.0017
$\delta$	0.2520	0.1311	0.2067	0.1777	1.0004
$\theta$	0.1569	0.1436	0.2067	0.0571	1.0004
${\sigma }_{\varphi }$	0.0728	0.0411	0.1871	0.0428	1.0003
${\sigma }_{\psi }$	0.4171	0.0944	0.4211	0.4001	0.9999
${\sigma }_u$	0.2576	0.1217	0.2347	0.1785	1.0002
${\sigma }_v$	0.0953	0.0681	0.3273	0.0666	1.0220
${\sigma }_a$	0.0163	0.0057	0.0100	0.0125	1.0002
${\sigma }_b$	0.0416	0.0209	0.0152	0.0277	1.0041
${\varphi }_1$	0.0013	0.0700	0.5149	0.0010	1.0000
${\psi }_1$	0.3549	0.1364	0.2138	0.3571	1.0551
$\mathit{lp}$ –	5016.7250	148.9654	–	–	0.9999

The random-effects ${\omega }_{\mathit{ij}}$ , $i,j=1,\dots ,N$ , deliver the heterogeneity in the variance and covariance equations. We have a special concern about their features in our analysis. The matrix of random effects is commonly supposed to follow the Wishart distribution. This traditional case plays an important role in the multivariate analysis, connected to inference about covariance matrices. Our modelling methodology uses a practical algorithm introduced by Odell and Feiveson (1966) to derive an $N$ -dimensional Wishart distribution with the degrees-of-freedom $N$ and the scale matrix ${\text{S}}_w$ . In the Bayesian contexts, ${\text{S}}_w$ interprets uncertainty about $\text{W}$ , ${\text{S}}_w$ evaluates the deviation or confusion, and $N$ interprets a degree of belief in ${\text{S}}_w$ . To operationalize the computation process, let ${\text{S}}_w$ be compound symmetry with parameters ${\sigma }_a^2$ and ${\sigma }_b^2$ , which allows for setting up a specific dispersion of modelling heterogeneity in the variance and the same dispersion for the heterogeneity in modelling the covariance structure. We have $E(\text{W}|{\text{S}}_w,N)=N{\text{S}}_w$ , implying that $E(W_{\mathit{ii}})=N({\sigma }_a^2+{\sigma }_b^2)$ and $E(W_{\mathit{ij}})=N{\sigma }_b^2$ , for all $i$ and $j$ . From Table 5, both parameters ${\sigma }_a^2$ and ${\sigma }_b^2$ are significant, confirming the stock-index heterogeneity. It means that the analyst should consider the effect of volatility in stock on other indices.

A convergence diagnostic in Stan is the rhat. It is close to one in Table 5, showing that all chains have converged to the same stationary distribution. Stan uses a more conservative version of rhat than its usual form in other packages, such as Coda (Plummer et al.(2006)). By default, it divides each chain in half to identify non-stationary chains (Gelman et al.(2013); Carpenter et al.(2017)). Trajectories of all parameters estimated over simulations and changes within a reasonable tolerance level (not shown here) showed a well-mixed chain. Figure 2 shows the posterior predictive histogram.

Figure 2

Predictive posterior histograms

The individual high-frequency volatility $h_{\mathit{ii},t}$ is modelled via GARCH, while the low-frequency volatility, ${\tau }_{\mathit{ii},t}$ , is modelled via the low-rank thin-plate spline. Asset returns are affected by politics and economic conditions, while unobserved factors can cause volatility. Figures 3–5, respectively, illustrate the variances ${\sigma }_{\mathit{ii},t}$ , high-frequencies $h_{\mathit{ii},t}$ and low-frequencies ${\tau }_{\mathit{ii},t}$ for the log-return of 10 stocks. Over the given period, our results show that the smoothing has correctly depicted the ${\sigma }_{\mathit{ii},t}$ function and ${\tau }_{\mathit{ii},t}$ properly makes smoothness the covariance function ${\sigma }_{\mathit{ii},t}$ .

Figure 3

Time-varying variance components for the log-return of 10 US stocks

Figure 4

High-frequency volatility in the log-return of 10 US stocks

Figure 5

Low-frequency volatility in the log-return of 10 US stocks ( ${10}^3$ scaled up)

Figure 6

(A) Time-varying covariance and (B) Low-frequency co-volatility of the WMT and 9 other log-returns ( ${10}^3$ scaled up)

To study the co-volatility in the covariance equation, consider the WMT index as an example. Figure 6 illustrates the low-frequency co-volatility and the time-varying covariance for the log-return of Walmart's stock and the nine other stocks. Low-frequency plots display the smoothness of covariances patterns. Figures describe the long-term changes of the WMT when other stocks have presented. Short-term volatility immediately affects the log return, while long-term volatility plays a role in the decision for future trades. The covariance structure is more influenced by long-term than short-term changes for this dataset. As seen here, coefficients in the low-frequency model are significant, unlike in the high-frequency model. It means that the volatility in one stock has a long-term impact on other stocks. When unexpected changes in one stock immediately impact other stocks, there is higher economic instability in the financial markets. Figure 7 presents the temporal heatmap for the time-varying covariance matrix for some selected days over this period.

Figure 7

Temporal heatmap for the time-varying covariance matrix

Forecasting volatility is particularly important for companies concerned with portfolio management, and brokers will trade stocks more confidently. There are various methods to estimate or predict volatility. The accurate estimate of a future market trend confirms the proposed model is appropriately fitted. To evaluate prediction accuracy, we first remove the vector observed at the time $T+1$ from the dataset, fit the model for the current period, and then estimate the observation vector, ${\mathbf{ϵ}}_{T+1}$ . The difference between predicted and true values can evaluate the appropriate fitting (Gelman et al.(2013)). Table 6 shows the observed values fall between $25\%$ and $75\%$ of the credible bounds. The $p$ -value for the Chi-squared test equals $0.075$ , implying the difference between observed and predicted values is significant.

Table 6

The observed and predicted bounds of log-return

Name of Stock	${\mathbf{ϵ}}_{T+1}$	$25\%$	$75\%$
WMT	0.0060	$-0.0323$	0.0460
PG	$-0.0015$	$-0.0180$	0.0450
KO	0.0018	$-0.0124$	0.0179
PEP	0.0062	$-0.0047$	0.0739
PM	0.0045	0.0158	0.2640
COST	0.0027	$-0.0202$	0.0400
MO	0.0020	0.0064	0.1047
MDIZ	0.0002	$-0.0120$	0.0691
EL	0.0035	$-0.0185$	0.0892
CL	0.0026	$-0.0041$	0.0678

Figure 8

Heatmap of the predicted covariance matrix

Moreover, the multivariate analysis can predict the cross-covariance of stocks involving their latent short- or long-run effects. Figure 8 displays a one-step-ahead forecast of the covariance matrix. Closer pairs of stocks indicate a higher covariance and darker colours over the heatmap. It graphically represents cross-dependence at time $T+1$ ; the covariance magnitude becomes larger as the colours turn stronger.

5 Concluding remarks

The article introduced a flexible model for describing high- and low-frequency variations with novel data analysis methods. Specifically, we offered modelling strategies for cross-sectional time-series analysis with serially correlated innovation terms and cross-individual dependence. The related conditional variance and covariance components followed the spline-GARCH process. In particular, we proposed a non-parametric regression model of the penalized low-rank smoothing spline form to account for time trends in the variance and covariance equations. Our proposed strategy can provide flexible modelling on the low-frequency volatility in equity markets. The model can differentiate non-parametrically dynamic patterns of high- and low-frequencies of variance-covariance structural equations and incorporate economic features to predict stock market volatility based on the time-series evidence. Many financial and macroeconomic time-series data are often conditionally heteroscedasticity, meaning that traditional estimation methods produce inefficient estimates. To analyse such data, taking advantage of advanced numerical methods is required. The HMC sampler can estimate related variances and covariances components in the normal case using computational Bayes. Future studies aim to present flexible multivariate distributions to analyse such correlated data with skewed structures.

Supplementary materials

Supplementary materials for this article containing software codes are available from http://www.statmod.org/smij/archive.html, and the datasets are freely available from https://finance.yahoo.com.