Continuous-time Markov-switching GARCH Process with Robust and Efficient State Path and Volatility Estimation

TL;DR

This paper introduces a robust Bayesian method for estimating a continuous-time Markov-switching GARCH model, effectively handling irregular data and multiple volatility states, with applications to finance and biomedical signals.

Contribution

It develops a novel COMS-GARCH model with a multi-path sampling scheme and Bernoulli noise injection for improved robustness and efficiency in state and volatility estimation.

Findings

Accurately predicts volatility regimes in simulated data.

Demonstrates effective application to real currency exchange rate data.

Shows robustness and efficiency improvements through Bernoulli noise injection.

Abstract

We propose a continuous-time Markov-switching generalized autoregressive conditional heteroskedasticity (COMS-GARCH) process for handling irregularly spaced time series (TS) with multiple volatilities states. We employ a Gibbs sampler in the Bayesian framework to estimate the COMS-GARCH model parameters, the latent state path and volatilities. To improve the inferential robustness and computational efficiency for obtaining the maximum a posteriori estimates for the state path and volatilities, we suggest a multi-path sampling scheme and incorporate the Bernoulli noise injection in the computational algorithm. We provide theoretical justifications for the improved stability and robustness with the Bernoulli noise injection through the concept of ensemble learning and the low sensitivity of the objective function to external perturbation in the TS. We apply the proposed COMS-GARCH process and the computational procedure to simulated TS, a real currency exchange rate TS, and a real blood volume amplitude TS. The empirical results demonstrate that the COMS-GARCH process and the computational procedure are able to predict volatility regimes and volatilities in a TS with satisfactory accuracy.

Equations

$$d{G}_t=Y_t={\sigma }_{t}d{L}_t(s_t)$$

$$d{\sigma }_t^2=\alpha (s_t)dt-\beta (s_t){\sigma }_{t-}^{2}dt+\lambda (s_t){\sigma }_{t-}^{2}d[L,L{]}_{t-}$$

$$\mathrm{Pr}(s_t=j|s_{t-}=k)={\eta }_{jk}dt+o(dt)\text{mbox}forj\ne k$$

$$\mathrm{Pr}(s_t=k|s_{t-}=k)=1-{\sum }_{j\ne k}{\eta }_{jk}dt+o(dt).$$

$$G_i-G_{i-1}=Y_i={\sigma }_{i-1}\sqrt{\Delta t_i}{ϵ}_i,$$

$${\sigma }_i^2=\alpha (s_i)\Delta t_i+\left({\sigma }_{i-1}^2+\lambda (s_i)Y_i^2\right)\exp \left(-\beta (s_i)\Delta t_i\right),$$

$$\mathrm{Pr}(s_i=j|s_{i-1}=k)=1-\exp (-{\eta }_{jk}\Delta t_i)\text{mbox}forj\ne k,$$

$$\mathrm{Pr}(s_i=k|s_{i-1}=k)={\sum }_{j\ne k}\exp (-{\eta }_{jk}\Delta t_i).$$

$$f(Y_i|Y_1,\dots ,Y_{i-1},s_1,\dots ,s_i)=N(0,{\rho }_i^2),\text{mbox}where$$

$${\rho }_i^2$$

$$\approx {\sigma }_{i-1}^2\Delta t_i=\alpha (s_{i-1})\Delta t_{i-1}\Delta t_i+\Delta t_i\left({\sigma }_{i-2}^2+\lambda (s_{i-1})Y_{i-1}^2\right)\exp \left(-\beta (s_{i-1})\Delta t_{i-1})\right).\text{vspace}-12pt$$

$$\text{mbox}E(Y_i^2|Y_1,\dots ,Y_{i-1},s_1,\dots ,s_i)=\text{mbox}V(Y_i|Y_1,\dots ,Y_{i-1},s_1,\dots ,s_i)={\rho }_i^2\approx {\sigma }_{i-1}^2\Delta t_i,\text{vspace}-12pt$$

$$f(\Theta |\mathbf{\eta },\mathbf{Y},\Delta \mathbf{t},S)\propto \pi (\Theta )L(\Theta ,\mathbf{\eta }|\mathbf{Y},S)=\pi (\Theta ){\prod }_{i=1}^n{\rho }_i^{-1}\exp \left(-Y_i^2/(2{\rho }_i^2)\right),$$

$$f({\eta }_{1k},\dots ,{\eta }_{\nu k}|\Theta ,\mathbf{Y},\Delta \mathbf{t},S)=f({\eta }_{1k},\dots ,{\eta }_{\nu k}|S,\Delta \mathbf{t})\text{mbox}fork=1,\dots ,\nu ,$$

$$\propto \pi ({\eta }_{1k},\dots ,{\eta }_{\nu k})\prod _{s_{i+1}=k,s_i=k}^{n-1}\left(2-\nu +\sum _{j\ne \nu }\exp (-{\eta }_{jk}\Delta t_{i+1})\right)\prod _{j\ne k}\prod _{s_{i+1}=j,s_i=k}^{n-1}(1-\exp (-{\eta }_{jk}\Delta t_{i+1})),$$

$$f(s_i|S_{-i},\Theta ,\mathbf{\eta },\mathbf{Y},\mathbf{\Delta t})\propto {\xi }_{s_i,s_{i-1}}{\xi }_{s_{i+1},s_i}{\prod }_{t=i}^n{\rho }_t^{-1}\exp \left(-Y_t^2/(2{\rho }_t^2)\right),$$

$$\xi_{s_{i},s_{i-1}}\!=\!\begin{cases}\begin{aligned} \textstyle 2-\nu+\sum_{k\neq s_{i-1}}\exp(-\eta_{k,s_{i-1}}\Delta t_{i})&\mbox{ when }s_{i}=s_{i-1}\\ 1-\exp(-\eta_{s_{i},s_{i-1}}\Delta t_{i})&\mbox{ when }s_{i}\neq s_{i-1}\\ \end{aligned}\end{cases}\!\!\!\!\!;\mbox{ similarly for $\xi_{s_{i+1},s_{i}}$. When there }$$

$$f({\eta }_{21}|S,\mathbf{Y},\Delta \mathbf{t})\propto$$

$$f({\eta }_{12}|S,\mathbf{Y},\Delta \mathbf{t})\propto$$

$$f(s_i|S_{-i},\Theta ,\mathbf{\eta },,\Delta \mathbf{t})$$

$$\left|\frac{{\prod }_{j=i}^{n}f(Y_j|{\mathbf{Y}}_{j-1},s_1,\dots ,s_i=k_1,\dots ,s_j)}{{\prod }_{j=i}^{n}f(Y_j|{\mathbf{Y}}_{j-1},s_1,\dots ,s_i=k_2,\dots ,s_j)}-\frac{{\prod }_{j=i}^{i+b}f(Y_j|{\mathbf{Y}}_{j-1},s_1,\dots ,s_i=k_1,\dots ,s_j)}{{\prod }_{j=i}^{i+b}f(Y_j|{\mathbf{Y}}_{j-1},s_1,\dots ,s_i=k_2,\dots ,s_j)}\right|<ϵ$$

$${\sum }_{k=b}^{n-i}\left(p^{k-b}(1-p{)}^{b-1}{\sum }_{\tilde{\mathbf{Y}}\in \tilde{\mathcal{Y}}}f({\tilde{s}}_i|{\tilde{S}}_{-i},\Theta ,\mathbf{\eta },\tilde{\mathbf{Y}})\right).\text{vspace}-5pt$$

$$\frac{{\prod }_{j=i}^{i+b}f(Y_j|{\mathbf{Y}}_{j-1},s_1,\dots ,s_i=k_1,\dots ,s_j)}{{\prod }_{j=i}^{i+b}f(Y_j|{\mathbf{Y}}_{j-1},s_1,\dots ,s_i=k_2,\dots ,s_j)}\left|\frac{{\prod }_{j=i+b+1}^{n}f(Y_j|{\mathbf{Y}}_{j-1},s_1,\dots ,s_i=k_1,\dots ,s_j)}{{\prod }_{j=i+b+1}^{n}f(Y_j|{\mathbf{Y}}_{j-1},s_1,\dots ,s_i=k_2,\dots ,s_j)}-1\right|<ϵ$$

$$\log (L(\tilde{S},\Theta ,\mathbf{\eta },{\tilde{Y}}_j=Y_i+Y_{i+1}))=\frac{1}{2}\log ({\tilde{\rho }}_j^2)+\frac{(Y_i+Y_{i+1}{)}^2}{2{\tilde{\rho }}_j^2}+\text{mbox}const.$$

$$\log (L(\tilde{S},\Theta ,\mathbf{\eta },{\tilde{Y}}_j^{\prime }=Y_i^{\prime }+Y_{i+1}^{\prime }))=\frac{1}{2}\log ({\tilde{\rho }}_j^2)+\frac{(Y_i^{\prime }+Y_{i+1}^{\prime }{)}^2}{2{\tilde{\rho }}_j^2}+\text{mbox}const.,$$

$$\text{mbox}E\left(\log (L(S,\Theta ,\mathbf{\eta },{\tilde{Y}}_j=Y_i+Y_{i+1}))\right)=\frac{1}{2}\log ({\tilde{\rho }}_j^2)+\text{mbox}E\left(\frac{(Y_i+Y_{i+1}{)}^2}{2{\tilde{\rho }}_j^2}\right)=\frac{1}{2}\log ({\tilde{\rho }}_j^2)+\frac{1}{2}$$

$$\text{mbox}E\left(\log (L(S,\Theta ,\mathbf{\eta },{\tilde{Y}}_j^{\prime }=Y_i^{\prime }+Y_{i+1}^{\prime })\right)=\frac{1}{2}\log ({\tilde{\rho }}_j^2)+\text{mbox}E\left(\frac{(Y_i^{\prime }+Y_{i+1}^{\prime }{)}^2}{2{\tilde{\rho }}_j^2}\right)=\frac{1}{2}\log ({\tilde{\rho }}_j^2)+\frac{1}{2}+\frac{{\epsilon }^2}{{\tilde{\rho }}_j^2}.$$

$${\epsilon }^2/({\tilde{\rho }}_j^2)={\epsilon }^2/({\rho }_i^2+{\rho }_{i+1}^2)$$

$$\log \left(L(S,\Theta ,\mathbf{\eta },Y_{i+1})\right)=\frac{1}{2}\log ({\rho }_{i+1}^2)+\frac{Y_{i+1}^2}{2{\rho }_{i+1}^2}\text{mbox}and\log \left(L(S,\Theta ,\mathbf{\eta },Y_i)\right)=\frac{1}{2}\log ({\rho }_i^2)+\frac{Y_i^2}{2{\rho }_i^2}.$$

$$\text{mbox}E\left(\log \left(L(S,\Theta ,\mathbf{\eta },Y_{i+1})\right)\right)=$$

$$\text{mbox}E\left(\log \left(L(S,\Theta ,\mathbf{\eta },Y_i)\right)\right)=$$

$$\text{mbox}E\left(\log \left(L(S,\Theta ,\mathbf{\eta },Y_{i+1}^{\prime })\right)\right)=\frac{1}{2}\log ({\rho }_{i+1}^2)+\text{mbox}E\left(\frac{Y_{i+1}^{{}^{\prime }2}}{2{\rho }_{i+1}^2}\right)=\frac{1}{2}\log ({\rho }_{i+1}^2)+\frac{1}{2}+\frac{{\epsilon }^2}{2{\rho }_{i+1}^2}$$

$$\text{mbox}E\left(\log \left(L(S,\Theta ,\mathbf{\eta },Y_i^{\prime })\right)\right)=\frac{1}{2}\log ({\rho }_i^2)+\text{mbox}E\left(\frac{Y_i^{{}^{\prime }2}}{2{\rho }_i^2}\right)=\frac{1}{2}\log ({\rho }_i^2)+\frac{1}{2}+\frac{{\epsilon }^2}{2{\rho }_i^2}.$$

$${\epsilon }^2/(2{\rho }_i^2)+{\epsilon }^2/(2{\rho }_{i+1}^2)$$

$${\sum }_{k=0}^r{\epsilon }^2/(2{\rho }_k^2)\text{mbox}and{\epsilon }^2/(2{\sum }_{k=0}^r{\rho }_{i+k}^2),$$

$$\displaystyle\sum_{k=0}^{r}\frac{\varepsilon^{2}}{2{\rho}_{i+k}^{2}}\bigg{/}\frac{r\varepsilon^{2}}{2\sum_{k=0}^{r}{\rho}_{i+k}^{2}}=r^{-1}\sum_{k=0}^{r}{\rho}_{i+k}^{-2}\sum_{k=0}^{r}{\rho}_{i+k}^{2}=r^{-1}\sum_{k=0}^{r}\frac{\sum_{k=0}^{r}{\rho}_{i+k}^{2}}{{\rho}_{i+k}^{2}}>1$$

$$\text{mbox}{E}_S\left(L(\Theta ,S|\mathbf{\Delta }\mathbf{t},\mathbf{Y},{\Theta }^{(l-1)},{\mathbf{\eta }}^{(l-1))}\right)$$

$$=$$

$$\propto$$

$$E_S\left(L(\mathbf{\eta },S|\mathbf{\Delta }\mathbf{t},\mathbf{Y},{\Theta }^{(l-1)},{\mathbf{\eta }}^{(l-1)}\right)$$

Full text

Keywords: Bernoulli noise injection; continuous-time; ensemble learning; irregularly (unevenly) spaced Markov-switching; Maximum a posterior; stability and robustness

1 Introduction

Heteroskedasticity is a common issue in time series (TS) data. The generalized autoregressive conditional heteroskedasticity (GARCH) model is a popular discrete-time TS model that accommodates heteroskeasticity and estimates the underlying stochastic volatility. The GARCH process has been extensively studied from both the theoretical and practical perspectives (Nelson, 1990b, ; Nelson,, 1991; Glosten et al.,, 1993; Zakoian,, 1994; Sentana,, 1995). Depending on the settings and problems, practical TS data can be recorded at irregularly spaced time points, creating a demand for continuous-time models. Nelson, 1990a derived the conditions under which the discretized-time GARCH model converges in distribution to a bivariate non-degenerate diffusion process as the length of the discrete time intervals goes to zero. The fact that the limiting process consisting of two independent Brownian motions (that drives the underlying volatility process and the accumulated TS, respectively) contradicts the GARCH model’s intuition that large volatilities are feedback of large innovations. Corradi, (2000) applied different parameterizations as a function of the discrete time interval to GARCH$(1,1)$ and obtained both degenerate and non-degenerate diffusion limits. Wang, (2002) further showed the asymptotic non-equivalence between the GARCH model and the continuous-time bivariate diffusion limit except for the degenerate case in Corradi, (2000). Kluppelberg et al., (2004) proposed a COntinuous-time GARCH (CO-GARCH) model that replaces the Brownian motions by a single Lévy process and incorporates the feedback mechanism by modeling the squared innovation as the quadratic variation of the Lévy process. Regarding the inferences for the CO-GARCH process, there exist several approaches, such as the quasi-likelihood (Bollerslev et al.,, 1994), method of moments (MoM) (Kluppelberg et al.,, 2004), pseudo-likelihood (MALLER et al.,, 2008; Marín et al.,, 2015), and Markov chain Monte Carlo (MCMC) procedures (Müller,, 2010).

Though the GARCH process and its variants, including the CO-GARCH process, account for the conditional heteroskedasticity, studies by Lamoureux and Lastrapes, (1990); Mikosch and Starica, (2004); Bauwens et al., 2014a and others have shown that volatility predictions by GARCH-type models may fail to capture the true variation in the presence of regime changes in the volatility dynamics. To solve this issue, Gray, (1996) developed Markov-switching GARCH (MS-GARCH) model that employs a hidden discrete Markov chain to assign a state to each time point. This model generates a temporal state path to allow variations in volatility. Regarding the inferences for the MS-GARCH models, since the state path is unobservable whereas the GARCH model parameters are state dependent, regular likelihood-based approaches require summing over exponentially many possible paths and can be computationally unfeasible. Several alternatives exist that can deal with the problem more efficiently. For example, the collapsing procedures (Kim,, 1994; Gray,, 1996; Dueker,, 1997; Klaassen,, 2002) based on simplified versions of the MS-GARCH model incorporates recombination mechanisms of the state space are popular approaches. Haas et al., (2004) developed a new MS-GARCH model that is analytically tractable and allows the derivation of stationarity condition and the process properties. Augustyniak, 2014a employed a Markov Chain Expectation Maximization (MC-EM) approach. Bauwens et al., (2010) proposed a Bayesian MCMC method but it can be slow in convergence. Recent methods focus on efficient sampling of state paths. Elliott et al., (2012) introduced a Viterbi-based technique to sample state paths; Bauwens et al., 2014b proposed a particle MCMC algorithm; Billioa et al., (2016) used a multi-point sampler in combination with the forward filtering backward sampling technique. Both the likelihood-based and the MCMC estimations have been implemented in software, such as the R package MSGARCH (Ardia et al.,, 2019), to facilitate the applications of MS-GARCH models in practice.

To the best of our knowledge, there does not exist a MS-GARCH model for irregularly spaced TS nor a CO-GARCH model to handle multiple states. On the other hand, there is a practical need for continuous-time Markov-switching GARCH processes to analyze TS data that exhibit heteroskedasticity and multiple regimes, and are collected in irregularly spaced time points or on a near-continuous time scale. For example, a heart rate variability (HRV) TS can have multiple volatilities regimes due to different types of activities or stress levels, and are typically recorded on a millisecond scale. Seismic waves are another example of TS data that consist of wave types of different magnitudes and are collected with high-frequency (milliseconds). Seismic waves inform our understanding of earth’s interior structure and are also used to predict earthquakes. Financial data are often irregularly spaced in time due to weekend and holiday effects and are known to exhibit different volatility states – changing behaviour drastically from steadily trending to extremely volatile after a major event or news.

To fill the methodological gap and respond to the practical needs, we propose a COntinuous-time Markov-switching GARCH process, COMS-GARCH for short. The COMS-GARCH model employs the Lévy process to model volatility within each different state and the continuous-time hidden Markov chain to model switching between states. For inferences about the COMS-GARCH parameters, we propose a Gibbs sampler in the Bayesian framework. For maximum a posterior (MAP) estimation on the state path and volatilities, we develop a computational procedure with a multi-path sampling scheme and the Bernoulli noise injection (NI) to accelerate the optimization and improve the robustness of the predicted state path and volatilities. We provide theoretical justifications to the employment of the Bernoulli NI, which includes ensemble learning of the state path and lowered sensitivity of the objective function to small random external perturbation in the TS.

The proposed COMS-GARCH process can be used to analyze multi-state irregularly spaced TS, estimate the historical volatilities, identify an optimal state path, and forecast volatilities and states. Though there exist simple and model-free approaches such as realized volatility to estimate the historical volatility, these methods cannot systematically identify different volatility states. In addition, for some TS, such as the HRV TS and seismic wave TS, it may not be meaningful to calculate aggregate measures across time points to calculate the realized volatilities. Though the existing MS-GARCH model can also be used for identify different volatility regimens, it cannot analyze irregularly spaced TS. Finally, we expect the estimated state path and historical volatilities via our proposed COMS-GARCH model to be more robust, largely due to the computational procedure we design specifically for obtaining the MAP estimation of the state path and volatilities.

In what follows, we introduce the COMS-GARCH process in Sec 2 and the inferential approach and computational algorithm in Sec 3, including the Bernoulli noise injection technique and its the theoretical properties, and forecasting through the COMS-GARCH model. We demonstrate the applications of the COMS-GARCH process and the proposed computational procedure via the simulation studies in Sec 4, a real currency exchange rate TS and a real blood volume amplitude TS in Sec 5. The paper concludes in Sec 6 with some final remarks.

2 COMS-GARCH Process

We propose the COMS-GARCH process based upon the CO-GARCH(1,1) process (Kluppelberg et al.,, 2004). To the best of our knowledge, CO-GARCH$(1,1)$ is the only Lévy-process driven CO-GARCH model that has analytical solutions for the model parameters from the stochastic differential equations and is inference-capable in the context of pseudo-likelihood. Brockwell et al., (2006) theoretically analyzed the CO-GARCH$(p,q)$ model driven by the Lévy process for general $p$ and $q$ values, but unable to obtain inferences for the model parameters.

Let $G_{t}$ for $t\in(0,T)$ denote the observed TS, and $L$ refer to the innovation, modeled by a Lévy process. Our proposed COMS-GARCH process $(G,\sigma^{2},S)=(\{G_{t}\},\{\sigma^{2}_{t}\},\{s_{t}\})$ for $t>0$ is the solution to the following set of stochastic differential equations

[TABLE]

The increment of the Lévy process $dL_{t}(s_{t})$ in Eqn (1) is assumed standardized with mean 0 and variance 1, and $[L,L]_{t-}$ in Eqn (2) is its quadratic variation process. $\sigma^{2}_{t}$ is the underlying volatility process governing the state of $G_{t}$ at time $t$. Eqns (3) and (4) represent the hidden continuous-time Markov chain with $\nu$ discrete states and transition parameters ${\boldsymbol{\eta}}=\{\eta_{jk}\}$ that model the regime switching in the TS for $j,k\in\{1,\cdots,\nu\}$ (note that $\eta_{jk}$ in Eqns (7) and (8) is not a probability, and the parameter space for $\eta_{jk}$ is $(0,\infty)$ instead of $\in(0,1)$).

Next, we define a family of discrete-time processes that approximates the above continuous-time process $(G,\sigma^{2},S)$, following the methodological framework in MALLER et al., (2008). There are a couple of reasons for obtaining a discretized process. First, real-life observed TS data are often recorded in discrete time, whether irregularly spaced or regardless of how fine the time scale is. Second, the discretization allows us to take advantage of the well developed inferential approaches for discrete-time GARCH processes. We will show the discretized process converge to the COMS-GARCH process.

The discretization is defined over a finite time interval $[0,T]$ for $T>0$. Let $0=t_{0}<t_{1}<\cdots<t_{i}<\cdots<t_{n}=T$ be a deterministic sequence that divides $[0,T]$ into $n$ sub-intervals of lengths $\Delta t_{i}=t_{i}-t_{i-1}$ for integers $i=1,\ldots,n$. Let $G_{0}\!=\!0$ and $\epsilon_{i}$ be a first-jump approximation of the Lévy process (Kluppelberg et al.,, 2004). A discretized COMS-GARCH process $(G_{n},\sigma_{n}^{2},s_{n})=(\{G_{i}\},\{\sigma^{2}_{i}\},\{s_{i}\})$ satisfies

[TABLE]

Since $Y_{i}$ is obtained by differencing the observed $G_{i}$, it is also observed. To ensure the positivity of Eqn (6), we require $\alpha(k)$ and $\lambda(k)$ to be non-negative for all states $k\!=\!1,\ldots,\nu$. To reflect the general belief that dependence between two quantities at two time points diminishes as the time gap increases, we also impose positivity on $\beta(k)$ for all states. As $n\rightarrow\infty$, $\Delta t_{i}\rightarrow 0$ and the discretized COMS-GARCH process in Eqns (5) to (8) converges in probability to the COMS-GARCH process defined in Eqns (1) to (4), as stated in Lemma 1.

Lemma 1** (Convergence of discretized COMS-GARCH process).**

Let $(G,\sigma^{2},s)$ be the COMS-GARCH process on time interval $[0,T]$, and $(G_{n},\sigma^{2}_{n},s_{n})$ be its discretized process. As $n\rightarrow\infty$, $\Delta t_{i}\rightarrow 0$ for $i=1,\ldots,n$ and $(G_{n},\sigma^{2}_{n},s_{n})$ converges in probability to $(G,\sigma^{2},s)$ in that the Skorokhod distance $D_{S}((G_{n},\sigma^{2}_{n},S_{n}),(G,\sigma^{2},s))\overset{p}{\rightarrow}0\mbox{ as }n\rightarrow\infty$.

The converges in probability in Lemma 1 also implies the convergence of $(G_{n},\sigma^{2}_{n},s_{n})$ in distribution to $(G,\sigma^{2},s)$. Lemma 1 is an extension of the theorem on the convergence of a discretized CO-GARCH process in Maller et al., (2008). The additional complexity in the COMS-GARCH process case is that it accommodates multiple states and has state-dependent GARCH parameters, both of which have no material impact on the discretization of the process and the underlying conditions for the convergence of discretized GARCH process. Therefore, the theoretical result of the convergence of the discretized COGARCH process n Maller et al., (2008) can be directly extended to the discretized COMS-GARCH process.

Remark 1**.**

The CO-GARCH model can be regarded as a special case of the COMS-GARCH process with the number of states $\nu=1$. The inferential approaches proposed for COMS-GARCH in Sec 3.1 and 3.2 and some theoretical results in Sec 3.4 also apply to CO-GARCH.

3 Inferences for COMS-GARCH Process

The parameters in the COMS-GARCH process include $\Theta=\{\alpha(k),\beta(k),\lambda(k)\}\;\forall\;k=1,\ldots,\nu$ and transition parameters $\boldsymbol{\eta}$. In addition to $\Theta$ and $\boldsymbol{\eta}$, we are also interested in learning the latent state $s_{i}$ and volatility $\sigma^{2}_{i}$ for $i=1,\ldots,n$ so to better understand an observed TS and to aid prediction of future states and volatilities.

As presented in Sec 1, there exist several approaches for obtaining inferences about the CO-GARCH and MS-GARCH models, respectively, such as the MoM, the maximum likelihood estimation (MLE) based on quasi-likelihood and the pseudo-likelihood, and the Bayesian MCMC framework for the CO-GARCH model, and the collapsing procedures, MLE via the EM or MC-EM algorithms, and Bayesian MCMC algorithms for the MS-GARCH model. For inferences about the proposed COMS-GARCH process, which has both the “MS” and the “CO” components, we propose a Bayesian Gibbs sampler coupled with the pseudo-likelihood as the inferential approach for the COMS-GARCH process. The rationales are as follows. The MoM for the CO-GARCH model does not work as it cannot deal with irregularly spaced TS. The collapsing procedure introduces an approximation into the MS-GARCH model and a bias in the likelihood (Augustyniak, 2014b, ). The MLE based on the pseudo-likelihood, pseudo-likelihood, and from the EM algorithms do not provide straightforward solutions on latent states $s_{i}$ and volatilities $\sigma^{2}_{i}$ which are also of inferential interests in addition to $\Theta$ and $\boldsymbol{\eta}$, whereas the Bayesian framework is rather convenient to achieve that goal.

In what follows, we first present the pseudo-likelihood for the COMS-GARCH process in Sec 3.1, then develop a Gibbs sampler for obtaining Bayesian inferences for the model parameters and estimating the latent states and volatilities from the COMS-GARCH process in Sec 3.2. We propose a robust and efficient computational procedure for obtaining the MAP estimates for the state path and volatilities in Sec 3.3 with multi-path sampling and the Bernoulli NI. A theoretical analysis of the effects of the Bernoulli NI on learning COMS-GARCH from a TS in Sec 3.4. A foresting procedure on future states and volatilities are provided in Sec 3.5.

3.1 Pseudo-likelihood for COMS-GARCH Process

Extending the pseudo-likelihood framework in Maller et al., (2008) for CO-GARCH to COMS-GARCH with multiple states, we have for $i=1,\ldots,n$, that

[TABLE]

Eqn (11) is obtained by taking the first-order Taylor expansion of Eqn (10) around $\Delta t_{i}=0$ and substituting $\sigma_{i-1}^{2}$ in Eqn (6). Eqns (9) to (11) suggest that

[TABLE]

which will be useful for foresting future volatilities as demonstrated in Sec. 3.2.

3.2 Gibbs Sampler for Bayesian Inferences about Model Parameters

If the goal of fitting the COMS-GARCH process to a TS is to obtain inferences for $\Theta$ and $\boldsymbol{\eta}$ given the pseudo-likelihood, the EM and MC-EM algorithms would make good choice to obtain Pseudo-MLE for $\Theta$ and $\boldsymbol{\eta}$, integrating out the latent states (we provide the EM and MC-EM algorithms in the Supplementary Materials for interested readers). As mentioned, our interest is also in estimating and forecasting states and volatilities, methods that rely on integrating out the latent states do not work well, whereas the Bayesian framework provides a convenient and straightforward approach. We propose a Gibbs sampler to obtain Bayesian inferences for the COMS-GARCH process.

Define $\boldsymbol{\Delta t}\!=\!(\Delta t_{1},\ldots,\Delta t_{n}),\mathbf{Y}\!=\!(Y_{1},\ldots,Y_{n}),S\!=\!(s_{1},\ldots,S_{n})$, and $\mathcal{S}$ is the set of all possible state paths. Denote the priors for $\Theta$ and $\boldsymbol{\eta}$ by $\pi(\Theta,\boldsymbol{\eta})$ and assume $\pi(\Theta,\boldsymbol{\eta})\!=\!\pi(\Theta)\pi(\boldsymbol{\eta})$. The conditional posterior distributions of $\Theta,\boldsymbol{\eta}$, and the states are respectively

[TABLE]

where $\rho_{i}^{2}$ in Eqn (13) is a function of $\Theta$ (see Eqn (11)), $\sum_{j=1}^{\nu}\eta_{jk}=1$ in Eqn (14), and

[TABLE]

are two states ($\nu=2$), Eqns (14) and (15) can be simplified to

[TABLE]

where $\xi_{1,s_{i-1}}\!=\!\exp(-\eta_{21}\Delta t_{i})$ if $s_{i-1}\!=\!1$, and $1-\exp(-\eta_{21}\Delta t_{i})$ if $s_{i-1}\!=\!2$; $\xi_{1,s_{i}}\!=\!\exp(-\eta_{21}\Delta t_{i+1})$ if $s_{i}=1$, and $1-\exp(-\eta_{21}\Delta t_{i+1})$ if $s_{i}=2$.

The Gibbs sampler draws samples on $\Theta,\boldsymbol{\eta}$ and $s_{i}$ for $i=1,\ldots,n$ alternatively from Eqns (13), (14), and (15). Upon convergence, after burning and thinning, we will have multiple, say $M$, sets of posterior samples of $\Theta,\boldsymbol{\eta}$, based on which their posterior inferences can be obtained. We will also have $M$ sets of samples on state $s_{i}$ and can calculate the posterior volatility $\sigma_{i}^{2}$ at each time point via Eqn (6). Connecting the states across the $n$ times points from each set of the state posterior samples leads to a state path. Due to the large sample space (totally $\nu^{n}$ possible paths), it is difficult to identity the MAP estimate for the state path out of the $M$ paths with acceptable accuracy unless $M\!>\!\!>\!\mu^{n}$ and a significant portion of paths have close-to-0 posterior probabilities with a few paths having significantly higher probabilities compared to the rest. To solve this issue, we design a new computational algorithm (reSAVE) as detailed next.

3.3 MAP Estimation for State Path and Volatility

To deal with the computational challenge in obtaining the MAP estimates for state path via the Gibbs sampler in Sec 3.2, we propose an inferentially Robust and computationally Efficient iterative procedure for State path And Volatility Estimation (reSAVE). The inferential robustness for the MAP estimates of the reSAVE procedure is brought by the Bernoulli NI implemented in each iteration of procedure, leading to both ensemble learning and improved stability of the object function from which the MAP estimates are obtained (see Sec 3.4 for details). The computational efficiency of the reSAVE can be attributed to a couple of factors: the Bernoulli NI that generates a sub-TS (smaller data size) in each iteration, of sampling of a small set of state $m$ to calculate MAP estimates for the state path, and the employment of a maximization-maximization scheme to obtain the MAP estimates of the parameters $\Theta$ and $\boldsymbol{\eta}$ and for the state path in each iteration. The steps of the reSAVE procedure are listed in Algorithm 1. Remarks 2 and 3 provide further remarks and explanations on some steps of the procedure.

Remark 2**.**

Regarding the number of iterations $N$, either $N$ can be prespecified or a convergence criterion can be used, such as the $l_{1}$ distances on the MAP estimates (e.g., $\{|\hat{\boldsymbol{\eta}}^{(l+1)}-\hat{\boldsymbol{\eta}}^{(l)}|,|\hat{\Theta}^{(l+1)}-\hat{\Theta}^{(l)},|S^{(l+1)}-S^{(l)}|\}$) or the objective functions between two consecutive iterations. If the distances are below a certain threshold, then the algorithm converges and can stop. Regarding the number of sampled state paths $m>1$, $m$ too small will not lead to stable MAP estimates; $m$ too large would increase the computational costs. In the simulation and case studies in Secs 4 and 5, we used $m=6$, which is deemed good enough. In general, we suspect the choice of $m$ relates to $n$ – the number of time points in the observed TS; the larger $n$ is, the larger $m$ is. The ensemble size $b$ refers to the number of observations following a given time point $i$ that are used to update the conditional posterior distribution of $s_{i}$. The specification of $b$ is mainly for computational efficiency consideration and is optional (see Proposition 1).

Remark 3**.**

The MAP estimates of $\Theta$ and $\boldsymbol{\eta}$ can be determined either through direct optimization of their respective conditional posterior distributions, or via MC approaches based on samples of $\Theta$ and $\boldsymbol{\eta}$ from their conditional posterior distributions. The MAP estimate of the state path in each iteration is defined as the path, out of the sampled $m$ paths, that maximizes the conditional posterior distribution of $S$, which is proportional to $\prod_{i=1}^{n}\rho_{i}^{-1}\exp(-y_{i}^{2}/(2\rho_{i}^{2}))\eta_{i,i-1}$ (Bauwens et al.,, 2010), given the latest MAP estimates $\Theta$ and $\boldsymbol{\eta}$ .

In Algorithm 1, if we set $m$ at 1 and remove line 10 for identifying the MAP estimate for $S$, and replace the calculation of conditional MAP estimates $\hat{\Theta}_{\mbox{\tiny{MAP}}}$ and $\hat{\boldsymbol{\eta}}_{\mbox{\tiny{MAP}}}$ in line 3 by random sampling $\Theta$ and $\boldsymbol{\eta}$ from their respective conditional posterior distributions, then it basically is the Gibbs sampler presented in Sec 3.2 for obtaining the Bayesian inferences about $\Theta$ and $\boldsymbol{\eta}$, with one important difference – the usage of Bernoulli NI to generate a sub-TS in every iteration. Though the Bernoulli NI is designed more for achieving ensemble learning and improving the stability of the objective functions for the state path optimization (see Sec 3.4), we expect its usage also helps in making inferences on $\Theta$ and $\boldsymbol{\eta}$ more robust, especially if $\nu$ is relatively large.

A key step in Algorithm 1 is the generation of sub-TS via the Bernoulli NI. The rationale for creating sub-TS when estimating the state path and volatilities is that the estimation can be sensitive to the TS data for the COMS-GARCH process. We employ the Bernoulli NI to create an ensemble of sub-TS’ of considerable diversity among the ensemble members across iterations so to reduce the sensitivity (see Sec 3.4 for more details). Algorithm 2 lists the steps of the Bernoulli NI. Specifically, in each iteration of Algorithm 1, we apply the Bernoulli NI to the observed TS $\mathbf{G}=\{G_{i}\}_{i=0}^{n}$ to obtain a sub-TS $\tilde{\mathbf{G}}$. Denote the sequence of time gaps in $\tilde{\mathbf{G}}$ by $\Delta\tilde{\mathbf{t}}$. With sub-TS $\{\tilde{\mathbf{G}},\Delta\tilde{\mathbf{t}}\}$, we only need to sample and update the states $\tilde{S}$ at the retained time points in each iteration, and the states of the dropped time points are kept at the values from the previous iteration, saving costs computationally.

Remark 4**.**

The differenced $\tilde{\mathbf{Y}}$ in a sub-TS after the Bernoulli NI is a summation of a sequence of differenced $\mathbf{Y}$ formed with the dropped observations in the original TS.

Remark 4 is a simple but interesting fact. For example, if $G_{i+1}$ gets dropped from the sequence of $\ldots,G_{i},G_{i+1},G_{i+2},\ldots$, then $\tilde{Y}_{i^{\prime}}=G_{i+2}-G_{i}=(G_{i+2}-G_{i+1})+(G_{i+1}-G_{i})=Y_{i+2}+Y_{i+1}$; say $r$ observations are dropped between $G_{i}$ and $G_{i+r+1}$, then $\tilde{Y}_{i^{\prime}}=G_{i+r+1}-G_{i}=(G_{i+r+1}-G_{i+r})+(G_{i+r}-G_{i+r-1})+\cdots+(G_{i+1}-G_{i})=Y_{i+r+1}+\cdots+Y_{i+1}$. This fact is used in the proof of Proposition 2 in Sec 3.4. Since the NI rate $p$ is usually small and the times points are dropped from the original TS randomly, with the fine time scale on which the TS is collected, the COMS-GARCH process can “digests” these “missing” time points effortlessly, without needing an ad-hoc approach to handle these dropped time points. The full conditional distributions of $\Theta,\boldsymbol{\eta}$, and states $\{s_{i}\}$ given the sub-TS in each iteration are given in Eqns (13) and (15) by replacing the original TS $(\mathbf{Y},\Delta{\mathbf{t}})$ with the sub-TS $(\tilde{\mathbf{Y}},\Delta{\tilde{\mathbf{t}}})$.

The Bernoulli NI for COMS-GARCH is inspired by the dropout technique for regularizing neural networks (NNs) (Srivastava et al.,, 2014), which injects Bernoulli noises to the NN structure (hidden and input nodes), leading to the $l_{2}$ regularization on model parameters. The Bernoulli NI we propose here is different procedurally in that it is applied to the observed data and drops randomly selected time points in the original TS in each iteration, rather than generating sub-models; in other word, the COMS-GARCH model remains the same throughout. The benefits of applying the Bernoulli NI in Algorithm 1 include reduced computational cost and its connection with ensemble learning and inferential stability and robustness. The Bernoulli NI also bears some similarity to bagging (Breiman,, 1996), a well-known ensemble learning algorithm, but also differs from the latter in two aspects. First, the Bernoulli NI leads to a random sub-TS (without replacement) of the original TS in each iteration of Algorithm 1 whereas bagging often generates a bootstrapped sample set with replacement that is of the same size as the training data. Second, bagging often generates multiple sets of samples, trains a model on each set in parallel, and then ensembles them into a meta-model, whereas the ensemble learning brought by the Bernoulli NI to the MAP estimation for COMS-GARCH is implicit, iterative, and realized sequentially. Finally, the Bernoulli NI we propose is not the same as the down-sampling technique in TS signal processing (Oppenheim et al.,, 1999). Signal processing aims at extracting useful features from collected signals, where down-sampling is used there for data reduction, compression, memory conservation, among others.

To choose Bernoulli NI $p$, we can apply a $k$-fold cross-validation (CV) procedure as listed in Algorithm 3.

Remark 5**.**

We predict $Y_{i}^{2}\approx\hat{\sigma}_{i-1}^{2}\Delta t_{i}$ via Eqn (12), where $\hat{\sigma}_{i-1}^{2}$ is the volatility at the immediately preceding time point that belongs to the training set and can be calculated from Eqn (6) given $s_{i-1}$ and $Y^{2}_{i-1}$. Once $\hat{Y}^{2}_{i}$ is obtained, $\hat{\sigma}^{2}_{i}$ can be back-calculated from Eqn (6) given $\hat{Y}^{2}_{i}$ and the interpolated state $s_{i}$.

Remark 6**.**

We apply the one-standard-error rule when choosing the Bernoulli NI rate $p$ instead of using the one minimizes the CV error due to two reasons. First, the one-standard-error rule is commonly used in selecting tuning parameters through CV in statistical machine learning to further mitigate over-fitting and improve generalization of the trained model. Second, there are always concerns on dependency between training and validating sets no matter what procedure is used for partitioning a TS into training and validating sets when developing CV procedures for TS data. While some partitioning methods might lead to less dependency than others, this is often achieved by throw away some data points (Hjorth,, 1982; Marron,, 1987; Chu and Marron,, 1991). We conjecture that the application of the one-standard-error rule (or something even harsher) helps alleviate the dependency concerns and leads to more generalizable parameter estimation and more robust prediction.

3.4 Theoretical Analysis on Inferential Benefits of Bernoulli NI

We state briefly in Sec 3.3 that the Bernoulli NI aims at improving the efficiency and robustness of the MAP estimates of the state path and thus the volatilities. In this section, we investigate theoretically the inferential benefits of the Bernoulli NI in two aspects. First, we show that the Bernoulli NI assists in the MAP estimation in an ensemble learning fashion; second, we establish that Bernoulli NI stabilizes the objective function in the presence of random external perturbation in the original TS.

3.4.1 Ensemble Learning of State Path

Ensemble learning aims for better prediction by assembling or combining a diverse group of learned models given a training data set (Opitz and Maclin,, 1999; Sollich and Krogh,, 1996; Kuncheva and Whitaker,, 2003; Polikar,, 2006). Many machine learning techniques are built on ensemble learning such as boosting and bagging. We show that, through the iterative Algorithm 1, the Bernoulli NI leads to sequential and implicit ensemble learning of the parameters and states for the COMS-GARCH model. The formal results are presented in Proposition 1.

Proposition 1** (ensemble learning of state path with Bernoulli NI).**

Assume that for $\forall\;\epsilon>0$, $\exists b\in N+$ such that

[TABLE]

$\forall k_{1}\neq k_{2}\in\{1,\ldots,\nu\},i\leq n-b,\mathbf{Y}_{j-1}=(Y_{1},\ldots,Y_{j-1})$. There exist $C_{k-1}^{b-1}$ ways to yield a set of b observations from a sequence of $k\in[b,n-i]$ consecutive observations. Denote the ensemble of the resultant $C_{k-1}^{b-1}$ sub-TS’ by $\tilde{\mathcal{Y}}$. Given a Bernoulli NI rate $p$, the conditional posterior distribution of $s_{i}$ given the ensemble $\tilde{\mathcal{Y}}$ is

[TABLE]

The proof of Proposition 1 is straightforward. Eqn (6) implies that the conditional distribution of $Y_{j}$ depends only on its variance since its mean is fixed at 0. Eqn (9) suggests that the impact of state $s_{i}$ on $\sigma^{2}_{j-1}$ (and thus $\rho_{j}^{2}$) decreases as $i$ departs from $j$ given the recursive formula on $\sigma^{2}$. Taken together, it implies that the state at time $t_{i}$ has minimal effect on the distribution of $Y_{j}$ at a future time point $t_{j}$ once the distance $t_{j}-t_{i}$ surpasses a certain threshold, which we use $b$ to denote. Mathematically, it means the ratio between $\prod_{j=i+b+1}^{n}f(Y_{j}|\mathbf{Y}_{j-1},s_{1},\ldots,s_{i}\!=\!k_{1},\ldots,s_{j})$ and $\prod_{j=i+b+1}^{n}f(Y_{j}|\mathbf{Y}_{j-1},s_{1},\ldots,s_{i}\!=\!k_{2},\ldots,s_{j})$ is arbitrarily close to 1, or

[TABLE]

for any $\epsilon>0$, leading to Eqn (19). $k$ given $b$ and $p$ follows a negative binomial distribution, leading directly to Eqn (20).

Taken together with the posterior distribution of $s_{i}$ in Eqn (15), Eqn (19) implies the posterior distribution of $s_{i}$ can be almost surely determined by the $b$ observations in TS $\mathbf{Y}$ that immediately follow $t_{i}$; that is, $f(s_{i}|S_{-i},\Theta,\boldsymbol{\eta},\mathbf{Y})=f(s_{i}|S_{-i},\Theta,\boldsymbol{\eta},Y_{i},Y_{i+1},\ldots,Y_{i+b})$. This narrow focus on just $b$ observations is undesirable from an inferential perspective especially when $b$ is small because the inference about the state path can be unstable and highly sensitive to even insignificant fluctuation in the TS. The Bernoulli NI helps mitigate this concern by diversifying the set of the $b$ observations.Specifically, after the Bernoulli NI, the posterior probability of $s_{i}$ in Algorithm 1 is a weighted average of the posterior distributions over multiple sets of $b$ observations with different compositions across iterations, as suggested by Eqn (20), leading to more robust state estimation.

Figure 1 provides a visual illustration on the ensemble effect achieved through the Bernoulli NI. The dashed line in each plot is associated with the right $y$-axis, which is the size (how many numbers) of an ensemble. The solid lines are associated with the left $y$-axis, which represent the weights assigned to ensembles of different sizes.

When there is no Bernoulli NI, $p=0$ and $k=b$, the ensemble is of size 1 (the first point on the dashed line in each plot). For $p>0$, we have more than one way of generating the set of $b$ observations; and the actual ensemble set size depends on $p$ and $k$. In brief, for a fixed $b$, as $k$ increases, the size of the ensemble set $\tilde{\mathcal{Y}}$, $C_{k-1}^{b-1}$, increases dramatically (the dashed line within each plot), implying more sub-TS’ are involved to obtain the posterior distribution of $s_{i}$. In addition, the ensemble set also increases dramatically with $b$ for a fixed $k-b$ value (the trend of the dashed lines across the 3 plots). The separated lines for different $p$ suggest that ensembles of different sizes are not weighted equally toward the conditional posterior distribution of $s_{i}$: the larger an ensemble, a smaller the weight it carries, especially for small $p$. As stated in Sec 3, the goal of employing the Bernoulli NI is to create an ensemble of sub-TS’ of considerable diversity among the ensemble members so to achieve more robust estimation. Figure 1 implies that $p$ as small as $O(0.01)$ can create an ensemble of sub-TS’ of enough diversity. A large $p$ leads to a larger ensemble, but it also implies higher computational cost which could overshadow the improved diversity. In addition, a large $p$ could also drop too many time points and lead to too much fluctuation in the sub-TS’ from iteration to iteration, leading to possibly large bias or large variance in the estimation. Our empirical studies in Sec 4 explore further the effect $p$ on the estimation of state path.

3.4.2 Stability Improvement in Objective Function

The Bernoulli NI also improves the stability of the objective functions from which the MAP estimation is obtained in the presence of random perturbation in the TS. The formal result is presented in Proposition 2. The proof is given in the Appendix

Proposition 2** (improved stability of objective function).**

Let $Y_{i}^{\prime}=Y_{i}+z_{i}$, with $z_{i}\overset{\text{ind}}{\sim}N(0,\varepsilon^{2})$ for $i=1,\ldots,n$, be an externally perturbed observation to the original observation $Y_{i}$ from TS $\mathbf{Y}$; and $Y_{i}^{\prime}$ comprises the perturbed TS $\mathbf{Y}^{\prime}$. Let $\tilde{\mathbf{Y}}$ and $\tilde{\mathbf{Y}}^{\prime}$ denote a sub-TS of $\mathbf{Y}$ and $\mathbf{Y}^{\prime}$, respectively, after implementing the Bernoulli NI in an iteration of Algorithm 1. The difference in the objection function (negative log-likelihood function or negative log-posterior distributions of $(\Theta,\boldsymbol{\eta}),S$ and $\boldsymbol{\sigma}^{2}$) given $\tilde{\mathbf{Y}}^{\prime}$ vs. that given $\tilde{\mathbf{Y}}$ after the Bernoulli NI is on average smaller than the difference obtained without the Bernoulli NI.

Since the objective function is more stable with the Bernoulli NI than without in the presence of random external perturbations, the MLEs of $(\Theta,\boldsymbol{\eta})$ and the MAP estimates of $(\Theta,\boldsymbol{\eta},S,\boldsymbol{\sigma}^{2})$ based on the former are expected to be more stable as well.

3.5 Foresting of Future State and Volatility

Foresting is often of major interest in TS analysis as they provide insights into future trends and are useful for decision making (e.g., developing option trading strategies in financial markets, predicting earthquakes). Algorithm 4 lists the steps for the $h^{th}$-step-ahead prediction of future volatilities and states through a trained COMS-GARCH process.

4 Simulation Studies

In this section, we demonstrate through simulation studies the applications of the COMS-GARCH process and illustrate the reSAVE algorithm in obtaining the MAP estimates of the model parameters, state path, and volatilities.

4.1 Simulation Study 2

We use this simulation study to show the improvement in the robustness of MAP estimates of volatilities brought by the Bernoulli NI in the reSAVE procedure in the COMS-GARCH process with one state. The COMS-GARCH process was simulated using the Lévy process as innovations (Kluppelberg et al.,, 2004), realized through the first-jump approximation that is normalized to have mean [math] and variance $1$. We characterized the time gap $\Delta t$ between two consecutive observations by a Poisson process with rate $\zeta$, that is, $\mbox{E}(\Delta t)=\zeta^{-1}$; and examined $\zeta=2.5,5,10,20$ (i.e., $\mbox{E}(\Delta t)=0.4,0.2,0.1,0.05$, respectively). The COMS-GARCH parameters are $\alpha=0.1\zeta,\beta=\log(10)\zeta,\lambda=\zeta$ The generated COMS-GARCH process is a discretized realization of the continuous-time process per Lemma 1 as $\Delta t\rightarrow 0$. We generated 50 TS repetitions and set the number of time points in each simulated TS at $n=500$; we also examined several Bernoulli NI rates $p=0.03,0.06,0.09,0.12$. To examine whether and how the Bernoulli NI help to improve the robustness of the volatility estimation, we perturbed each simulated TS by adding independently sampled noise from $N(0,0.1\!\cdot\!\mbox{SD}_{\mathbf{Y}})$ to each observation $Y_{i}$ for $i=1,\ldots,n$ in TS $\mathbf{Y}$, where SDY is the sample standard deviation of $\mathbf{Y}$.

The COMS-GARCH model was fitted to the the unperturbed and perturbed TS data in each repetition under each simulation scenario, following the steps of a simplified version of Algorithm 1 without the steps for state imputation or estimation of the transition parameters $\boldsymbol{\eta}$ as there is only one state. We imposed flat priors on $\Theta\!=\!(\alpha,\beta,\lambda)$. Per proposition 1, we used $b=20$ observations $\tilde{Y}_{i},\ldots,\tilde{Y}_{\min\{i+b,n\}}$ instead of the whole TS $\tilde{Y}$ when sampling $s_{i}$ from its full conditional posterior probability in each sub-TS to further save computational time. The number of iterations $N$ was set at 300. The convergence of the reSAVE procedure was examined by visual inspection of the trace plots on the MAP estimates of $\Theta$.

Figure 2 presents an example on the estimated volatility superimposed on top of the true volatility from a randomly chosen repetition (left plot) and how much the Bernoulli NI can help reduce the relative %|bias| of the estimated volatility ($l_{1}$-distance scaled by the true volatility and averaged over the 500 time points per TS and the 50 TS repetitions) (right plot).

In the left plot, the estimated volatilities almost completely overlap with the true volatility across all the time points, implying the reSAVE procedure performs well for this TS repetition. There are several observations from the right plot. First, the bias in the volatility estimation depends on $\zeta$ and thus $\mbox{E}(\Delta t)$. Second, the estimated volatilities can be sensitive to even mild fluctuation in the TS data as the bias with externally perturbed TS (crosses at $p=0$) increases compared to that at no external perturbation (circles at $p=0$). Third, the Bernoulli NI helps bring the bias down for all the examined $\zeta$ values, and in the case of $\zeta=2.5\;(\mbox{E}(\Delta t)=0.4)$ and $\zeta=5\;(\mbox{E}(\Delta t)=0.2)$, down to the level comparable to the unperturbed TS. There are not as much differences in the bias across $p$.

4.2 Simulation Study 2

In this simulation study, we demonstrate the inferential robustness and computational efficiency of the proposed reSAVE procedure in a two-state COMS-GARCH process, compared the estimation vs that from the MSGARCH(1,1) process. We simulated 50 TS repetitions. In simulating each TS ($n=1,000$), we used the same the CO-GARCH process as the first simulation study, but separately in each of the two state. Specifically, we characterized time gap $\Delta t$ by a Poisson process with rate $\zeta=10,40$, respectively. We set the GARCH model parameters at $\alpha(k)\!=\!\zeta c_{k},\beta(k)\!=\!-\zeta\log(c_{k}),\lambda(k)\!=\!\zeta$ for $k=1,2$. We let $c_{2}\!=\!2.5c_{1}$ and examined two values of $c_{1}$ at 0.1 and 0.025, respectively. The transitions of the states between two adjacent time points were modeled by a hidden Markov process and the transition probabilities were calculated per Eqns (7) and (8) given the transition parameters $\boldsymbol{\eta}$ and simulated time gaps from the Poisson process. We examined two sets of $\boldsymbol{\eta}=(\eta_{12},\eta_{21})$ at (0.1, 0.1) and (0.25, 0.25), respectively. For the Bernoulli NI rate, we examined $p=0.01,0.02,0.03$. The externally perturbed TS data were obtained by adding independently sampled noise from $N(0,0.1\!\cdot\!\mbox{SD}_{i})$ to each $Y_{i}$ for $i\!=\!1,\ldots,n$, where SDi is the sample SD of the subset of $\mathbf{Y}$ that are of the same state as $s_{i}$.

Algorithm 1 was applied to fit the COMS-GARCH model to each simulated TS in each simulation setting. The number of iterations was set at 2,000 and the number of sampled paths $m$ was set at 6. Per proposition 1, we used $b=20$ observations $\tilde{Y}_{i},\ldots,\tilde{Y}_{\min\{i+b,n\}}$ instead of the whole TS $\tilde{Y}$ when sampling $s_{i}$ from its full conditional posterior probability in each sub-TS to further save computational time. Non-informative flat priors were imposed on $\boldsymbol{\theta}$ and $\eta$. The MAP estimates of $(\Theta,\boldsymbol{\eta})$ in each iteration were obtained directly via an optimizer. The convergence of the algorithm was assessed by visual inspection of the trace plots on the MAP estimates of $(\Theta,\boldsymbol{\eta})$ and the logarithm of the joint posterior distribution of $(\Theta,\boldsymbol{\eta},\boldsymbol{\sigma}^{2},S)$. The benchmark method MS-GARCH(1, 1) was fitted to the TS data using R package MSGARCH. Since the MS-GARCH model assumes evenly-space TS data, we first applied linear interpolation to each simulated TS to obtain the equally spaced TS data before the fitting the model.

First, to illustrate the effects of the Bernoulli NI (with $b=20$) and the multiple path sampling scheme ($m=6$) in accelerating the convergence of the reSAVE procedure for estimating the MAP of the state path and volatilities, we compared the following three settings on convergence: 1) $m=6$ with Bernoulli NI at $p=0.02$ ($m=6,p=0.02$); 2) $m=6$ without Bernoulli NI ($m=6,p=0$); 3) $m=1$ (single path sampling) without Bernoulli NI ($m=1,p=0$). The results from 4 simulation scenarios are presented in Figure 3.

($m=6,p=0.02$) reached convergence within the least iterations ($\sim 800$); ($m=6,p=0$) reached convergence around iteration $1,000$ in the state path estimation and around iteration $1,600$ in the volatility estimation. For ($m=1,p=0$), convergence is reached around $1,800$ for $c_{1}=0.1$, but does not seem to converge at iteration $2,000$. In terms of the actual computational time for the reSAVE procedure, it is $O(Nd)$, where $N$ is the number of iteration and $d$ is the computational time per iteration. Though $b$ might be larger for $m>1$ with multi-path sampling and identifying the MAP path out of $m$, the longer per-iteration computational time in general does not overshadow the amount of time saved with the reduction in $N$, especially that $m$ is often a small number and parallel computing can used for sampling. In this simulation, the total computational time was less for $m=6$ on average than for $m=1$. We provide in the Supplementary Materials some examples on the MAP estimates for the state path and volatilities from a randomly selected TS analyzed in these three settings. It is obvious from the plot the setting of ($m=6,p=0.02$) has the most accurate estimation for the state path and volatilities across the time points.

We provide in Figure 4 an example on the MAP estimates via the reSAVE procedure on the state path and volatilities from a randomly selected TS simulated at $\zeta=10,c_{1}=0.1,\eta_{12}=\eta_{21}=0.1$ with Bernoulli NI ($p=0.02$) vs. without, along with the benchmark results from the MS-GARCH model. When there is no external perturbation in the TS, the Bernoulli NI does not negatively affect the volatility and state estimation; when there is, the Bernoulli NI improves the accuracy in the state and volatility estimation. This observation suggests that the Bernoulli NI is a “intelligent” technique, only doing its tricks when needed and is silent otherwise. As for the MS-GARCH model, the volatilities tend to be over-estimated and there are more mis-predicted states compared to using the COMS-GARCH process.

To examine the effectiveness of the reSAVE procedure in improving the state path and volatility estimation, we summarized in Figure 5 the state mis-prediction rate and the relative %|bias| of the estimated volatilities over the 1,000 time points per TS and the 50 TS repetitions . In all the examined simulation scenarios, the accuracy of the path identification and volatility estimation is significantly improved with a proper Bernoulli NI rate $p$ than without NI. The smallest mis-prediction rate (8% $\sim$ 17%) and volatility estimation bias (12% $\sim$ 25%) are achieved around $p$ at $0.01\sim 0.02$ in most scenarios; and further increasing $p$ does not seem to improve the prediction accuracy, implying that selecting a proper choice of $p$ – such as via the CV procedure in Algorithm 3 – is important. In addition, how much the Bernoulli NI helps in reducing the prediction bias relates to $\zeta$ (gaps between two adjacent time points), $\boldsymbol{\eta}$ (transition parameters), and the values of $\Theta$.

The biases and root mean squared errors (MSE) for the MAP estimates of COMS-GARCH parameters $\Theta$ and $\boldsymbol{\eta}$ are presented in Figure 6. The estimates for $\alpha(k)$ and $\beta(k)$ are generally accurate and are noticeable for some simulation scenarios for $\lambda(k)$ and $\boldsymbol{\eta}$. The relatively large bias for $\boldsymbol{\eta}$ can be at least partially attributed to the low transition probabilities between different states, leading to data sparsity in estimating $\boldsymbol{\eta}$ (we provide in the Supplementary Materials the histograms of $\Pr(s_{1}|s_{2})=\Pr(s_{2}|s_{1})$ for $\eta_{12}=\eta_{21}=0.1$ and 0.25, respectively. In both cases, the mode transition probabilities is close to 0). Regarding the estimation of $\Theta$, the existence of bias in the GARCH, CO-GARCH and MS-GARCH parameters is rather a common problem than something unique to the Gibbs sampler or the reSAVE procedure we propose for the COMS-GARCH process. Bollerslev et al., (1994) suggested that the MLE for the parameters from the GARCH model are biased; the estimation bias for the GARCH parameters for both the quasi-MLE and the constrained M-estimators (more robust) can be as large as 20% in the empirical studies in Mendes and De, (2000); that for the parameters for the CO-GARCH process in MALLER et al., (2008) can be as large as 30%. How to improve the statistical accuracy for the model parameters is on our future research agenda (more details are provided in Sec 6).

5 Applications

We apply the proposed COMS-GARCH and the reSAVE algorithm in two real data sets. The first is the exchange rate between the US dollar and Canadian dollar collected from July to December 2008 and can be downloaded from https://www.histdata.com. The second is blood volume amplitude (BVA) TS data and can be downloaded from https://archive.ics.uci.edu/ml/datasets/PPG-DaLiA#. In the exchange rate application, the volatility process of interest is the daily volatility of the log-return. To keep the data at a manageable size, we took every $90$-th observation and performed a log transformation on the exchange rate. The final TS $\mathbf{G}$ we worked with contains $n=1501$ times points from July to December. Figure 7(a) depicts the log-return TS $\mathbf{Y}$ (obtained by differencing $\mathbf{G}$), time gap $\Delta\mathbf{t}$ in days, and the histogram $\Delta\mathbf{t}$. The TS $\mathbf{Y}$ exhibits a change in volatility somewhere between September and October, an indication of multiple volatility states. The histogram suggests that $\Delta\mathbf{t}$ has a median around 0.07 days, but can be as short as just a couple of hours (the majority) and as long as $\geq 2$ days due to the weekend and holiday effect. In the BVA application, we extracted the measurements from the photoplethysmograph of the blood volume pulse (64 Hz; i.e., 64 times per second) by taking the valley and peak pulse values in each cycle and scaling them by $0.01$. We then took the first $1,000$ BVA measurements from a random patient as the input TS $\mathbf{Y}$. Figure 7(b) depicts the TS $\mathbf{Y}$, time gap $\Delta\mathbf{t}$ and its histogram in the BVA TS. The plot of TS $\mathbf{Y}$ suggests different volatilities, which can occur when a subject is in different physical conditions or experiences different emotional episodes. In summary, Figure 7 suggests that the time points are irregularly spaced, and the TS’ have multiple regimens in both cases, justifying the employment of the COMS-GARCH process to analyze such data.

The $\mathbf{Y}$ plot in the exchange rate application suggests that there are two volatility states. We thus adopted a 2-state COMS-GARCH process to estimate the state path and volatilities via the reSAVE algorithm with $p=0.02,b=20,m=6$ and $N=1,000$ iterations. The $\mathbf{Y}$ plot in the BVA application suggests that there are two volatility states. We adopted a 2-state COMS-GARCH process to estimate the state path and volatilities via the reSAVE algorithm with $p=0.02,b=20,m=6$ and $N=1,500$ iterations. The results from both applications are provided in Figure 8.

The estimated states and volatilities in general reflect well the two expected volatility states in both applications. For the exchange rate application, the estimated volatilities can be used to measure how far the exchange rate moves away from its mean value and help define risks and assess the volatility of the market when making investment decisions. For the BVA application, the estimated volatility and the state path can be used for evaluating how quickly a patient returns to the baseline (low) volatility state after spending a certain amount of time in the high-volatility state. A delay in returning hints potential cardiac risks in the patient. Compared to the simple methods such as the moving average (Hamburg and Benjamin,, 2009), the state and volatility predictions via the COMS-GARCH model are more interpretable as it assigns a definite state to each time point and does not have time-lagging issues when the TS switches regime. In both applications, the estimated state paths and volatilities can also be used as inputs for further analysis.

6 Discussion

We propose the COMS-GARCH process for handling irregularly spaced TS data with multiple volatility states. We also introduce the reSAVE procedure with the Bernoulli NI for obtaining the MAP estimates for model parameters, state path, and volatilities. The computational efficiency and inferential robustness of the reSAVE procedure are established and illustrated theoretically or empirically.

As discussed briefly when presenting the simulations results, there is a lack of in-depth theoretical investigation on the asymptotic properties of the MLE and MAP estimation, such as consistency, based on the pseudo-likelihood or quasi-likelihood as $n\rightarrow\infty$ and $T\rightarrow\infty$ for the recently developed MS-GARCH and CO-GARCH processes; and the existence of bias in the estimation of *-GARCH parameters is well acknowledged. This also creates opportunities for us to develop more accurate inferential procedures for *-GARCH models (COMS-GARCH included), coupled with investigation on their theoretical properties.

We conjecture that the reSAVE procedure is applicable not only to the COMS-GARCH and CO-GARCH processes but also to other solvable CO-*-GARCH processes. For example, it will make an interesting future topic to develop the COMS-Exponential-GARCH and COMS-Integrated-GARCH processes and examine the performance of the reSAVE procedure in these settings. We also expect that the reSAVE procedure can be used in the COMS-ARMA process for trend estimation, yielding some types of weighted $l_{2}$ regularization on the ARMA parameters. More work will be needed to prove the theoretical conjectures formally.

Appendix

Proof of Proposition 2

The conditional distribution of $Y_{i}$ given $\mathbf{Y}_{i-1},S,\Theta,\eta$ in the original TS is $N(0,\rho^{2}_{i})$ (Eqn 9), and that of $Y^{\prime}_{i}$ from given $\mathbf{Y}^{\prime}_{i-1},S,\Theta,\eta$ in the perturbed TS is $N(0,\varrho^{2}_{i})$, where $\mathbf{Y}^{\prime}_{i-1}=(Y^{\prime}_{1},\ldots,Y^{\prime}_{i-1})$, and $\varrho^{2}_{i}=\rho_{i}^{2}+\varepsilon^{2}$ for $i=1,\ldots,n$. Similarly, $f(Y_{i+1}|\mathbf{Y}_{i},S,\Theta,\eta)=N(0,\rho^{2}_{i+1})$ and $f(Y^{\prime}_{i+1}|\mathbf{Y}^{\prime}_{i},S,\Theta,\eta)=N(0,\varrho^{2}_{i+1})$. WLOG, assume $\tilde{Y}_{j}=Y_{i}+Y_{i+1}$ after Bernoulli NI in the original TS and the corresponding perturbed version is $\tilde{Y}^{\prime}_{j}=Y^{\prime}_{i}+Y^{\prime}_{i+1}$.

When analyzing the perturbed TS and original perturbed TS data after the NI via the same COMS-GARCH process, the contribution of observation $\tilde{Y}^{\prime}_{j}$ and $\tilde{Y}_{j}$ to the overall negative log-pseudo-likelihood value given $S,\Theta,\boldsymbol{\eta}$ is

[TABLE]

respectively, and their expectations are

[TABLE]

The difference between Eqns (.2) and (.1)

[TABLE]

Without Bernoulli NI, the contribution of $Y_{i+1}$ and $Y_{i}$ to the overall negative log-pseudo-likelihood on $S,\Theta,\boldsymbol{\eta}$, given the original perturbed data is, respectively

[TABLE]

And their expectations are

[TABLE]

respectively. Similarly, in the perturbed TS $\mathbf{Y}^{\prime}$,

[TABLE]

The difference in the sum of the negative conditional log-pseudo-likelihood given $Y^{\prime}_{i}$ and $Y^{\prime}_{i+1}$ vs. that given $Y_{i}$ and $Y_{i+1}$, that is, Eqns [(.6)+(.7)]-[(.4)+(.5)], is

[TABLE]

Finally, the ratio between Eqns (.3) and (.8) is $\frac{(2\rho_{i})^{-2}+(2\rho_{i+1})^{-2}}{\left(\rho_{i}^{2}+\rho_{i+1}^{2}\right)^{-1}}=\frac{2\rho_{i}^{2}\rho_{i+1}^{2}}{\left(\rho_{i}^{2}+\rho_{i+1}^{2}\right)^{2}}<1.$; in other words, the difference in the expected loss function between the original TS vs. that subject to random perturbation with Bernoulli NI is smaller that without Bernoulli NI.

In general, when applying the Bernoulli NI during the reSAVE procedure, a time point is dropped from the TS randomly with a probability $p$. Say $r$ observations are dropped between $G_{i-1}$ and $G_{i+r}$, then $\tilde{Y}_{j}=Y_{i+r}+\cdots+Y_{i}$. Eqns (.8) and (.3) now become

[TABLE]

respectively and their ratio is

[TABLE]

for any general $r\geq 1$ in each sub-TS that involves a dropped observation. Putting all the sub-TS’ together, the overall difference in the expected negative log-likelihood function with vs without perturbation, after Bernoulli NI, is smaller than that without Bernoulli NI.

In the framework of Bayesian modelling, the objective function becomes the log-posterior distribution, that is the sum of the log-likelihood function and the log-prior. Since the prior is the same with vs. without the external perturbation, it is cancelled out when the difference between two log-posterior distributions is taken, and we arrive at the same equations as Eqns (.3), (.8), and (.9); and the conclusion also holds in the Bayesian framework. ∎

**Supplementary Materials

**

\justify

Appendix S1 EM and MC-EM Algorithms for Parameter Estimation in COMS-GARCH process

The expectation step in iteration $l$ of the EM algorithm comprises the calculation of

[TABLE]

where $\boldsymbol{\Delta t}=(\Delta t_{1},\ldots,\Delta t_{n}),\mathbf{Y}=(Y_{1},\ldots,Y_{n})$, and $\mathcal{S}$ is the set of all possible paths. The M-step maximizes the expectations in Eqns (S1.1) and (S1.2) to obtain MLEs $\Theta^{(l)}$ and $\boldsymbol{\eta}^{(l)}$.

The classical EM estimation procedure can be computationally intensive or even unfeasible when $n$ is large. In the E-step, one would sum over $\nu^{n}$ possible state paths in Eqns (S1.1) and (S1.2). For example, if $\nu=2$ and $n=100$, the number of possible state paths is $3.27\times 10^{150}$, an astronomical figure to deal with in practice. On top of that, $\sigma^{2}_{i}$ and $\rho^{2}_{i}$ for a given path are calculated recursively for $i=1,\ldots,n$ (Eqns (6) and (11)). To circumvent the computational issue, the MC-EM algorithm can be used to obtain numerical approximations to the expected likelihood. Specifically, the expectations over $S\in\mathcal{S}$ in Eqns (S1.1) and (S1.2) are replaced by their respective averages over a set of path samples $S_{j}\!=\!(s_{1j},\ldots,s_{nj})$ for $j\!=\!1,\ldots,m$ in each iteration. Specifically, in the $l$-th iteration, Eqns (S1.1) and (S1.2) will be replaced by

[TABLE]

The path samples $S_{j}$ for $j=1,\ldots,m$ are drawn from the conditional distribution of state $s_{i}$ for $i=1,\ldots,n$ given $\mathbf{Y},\mathbf{\Delta t},\Theta^{(l)},\boldsymbol{\eta}^{(l-1)}$, and the states at other time points $S^{(l-1)}_{-i}=\left\{s^{(l-1)}_{1},\ldots,s^{(l-1)}_{i-1},s^{(l-1)}_{i+1},\ldots,s^{(l-1)}_{n}\right\}$. Specifically,

[TABLE]

Appendix S2 Additional Results From the Simulation Studies