Fisher’s z Distribution-Based Mixture Autoregressive Model

Abstract

We generalize the Gaussian Mixture Autoregressive (GMAR) model to the Fisher’s z Mixture Autoregressive (ZMAR) model for modeling nonlinear time series. The model consists of a mixture of K-component Fisher’s z autoregressive models with the mixing proportions changing over time. This model can capture time series with both heteroskedasticity and multimodal conditional distribution, using Fisher’s z distribution as an innovation in the MAR model. The ZMAR model is classified as nonlinearity in the level (or mode) model because the mode of the Fisher’s z distribution is stable in its location parameter, whether symmetric or asymmetric. Using the Markov Chain Monte Carlo (MCMC) algorithm, e.g., the No-U-Turn Sampler (NUTS), we conducted a simulation study to investigate the model performance compared to the GMAR model and Student t Mixture Autoregressive (TMAR) model. The models are applied to the daily IBM stock prices and the monthly Brent crude oil prices. The results show that the proposed model outperforms the existing ones, as indicated by the Pareto-Smoothed Important Sampling Leave-One-Out cross-validation (PSIS-LOO) minimum criterion.

1. Introduction

Many time series indicate non-Gaussian characteristics, such as outliers, flat stretches, bursts of activity, and change points (Le et al. 1996). Several methods have been proposed to deal with the presence of bursts and outliers such as applying robust or resistant estimation procedures (Martin and Yohai 1986) or omitting the outliers based on the use of diagnostics (Bruce and Martin 1989). Le et al. (1996) introduced a Mixture Transition Distribution (MTD) model to capture non-Gaussian and nonlinear patterns, using the Expectation–Maximization (EM) algorithm as its estimation method. The model was applied to two real datasets, i.e., the daily International Business Machines (IBM) common stock closing price from 17 May 1961 to 2 November 1962 and the series of consecutive hourly viscosity readings from a chemical process. The MTD model appears to capture the features of the data better than the Autoregressive Integrated Moving Average (ARIMA) models.

The Gaussian Mixture Transition Distribution (GMTD), which is a special form of MTD, was generalized to a Gaussian Mixture Autoregressive (GMAR) model by Wong and Li (2000). The model consists of a mixture of K Gaussian autoregressive components and is able to model time series with both heteroscedasticity and multimodal conditional distribution. It was applied to both the daily IBM common stock closing price from 17 May 1961 to 2 November 1962 and the Canadian lynx data for the period 1821–1934. The results indicated that the GMAR model was better than the GMTD, ARIMA, and Self-Exciting Threshold Autoregressive (SETAR) models.

The use of the Gaussian distribution in the GMAR model still leaves problems, because it is able to capture only short-tailed data patterns. Some methods developed to overcome this problem include the use of distributions other than Gaussian, e.g., the Logistic Mixture Autoregressive with Exogenous Variables (LMARX) model (Wong and Li 2001), Student t-Mixture Autoregressive (TMAR) model (Wong et al. 2009), Laplace MAR model (Nguyen et al. 2016), and a mixture of autoregressive models based on the scale mixture of skew-normal distributions (SMSN-MAR) model (Maleki et al. 2020). Maleki et al. (2020) proposed the finite mixtures of autoregressive processes assuming that the distribution of innovations belongs to the class of Scale Mixture of Skew-Normal (SMSN) distributions. This distribution innovation can be employed in data modeling that has outliers, asymmetry, and fat tails in the distribution simultaneously. However, the SMSN distribution’s mode was not stable in its location parameters (Azzalini 2014). In this paper, we propose a new MAR model called the Fisher’s z Mixture Autoregressive (ZMAR) model which assumes that the distribution of innovations belongs to the Fisher’s z distributions (Solikhah et al. 2021). The ZMAR model consists of a mixture of K-component Fisher’s z autoregressive models, where the numbers of components are based on the number of modes in the marginal density. The Fisher’s z distribution’s mode is stable in its location parameters, whether it is symmetrical or skewed. Therefore, Fisher’s z uses the errors in each component of the MAR model to capture the ‘most likely’ mode value—(not the mean, median, or quantile) of the conditional distribution Y_t given the past information. The conditional mode may be a more useful summary than the conditional mean when the conditional distribution of Y_t given the past information is asymmetric. Other distributions that also have a stable mode in its location parameter are the MSNBurr distribution (Iriawan 2000; Choir et al. 2019; Pravitasari et al. 2020), the skewed Studen t distribution (Fernández and Steel 1998), and the log F-distribution (Brown et al. 2002).

The Bayesian technique using Markov Chain Monte Carlo (MCMC) is proposed to estimate the model parameters. Among the algorithms in the MCMC, the Gibbs sampling (Geman and Geman 1984) and the Metropolis (Metropolis et al. 1953) algorithms are widely applied and well-known algorithms. However, these algorithms have slow convergence due to inefficiencies in the MCMC processes, especially in the case of models with many correlated parameters (Gelman et al. 2014, p. 269). Furthermore, Neal (2011) has shown that the Hamiltonian Monte Carlo (HMC) algorithm is a more efficient and robust sampler than Metropolis or Gibbs sampling for models with complex posteriors. However, the HMC suffers from a computational burden and the tuning process. The HMC can be tuned in three places (Gelman et al. 2014, p. 303), i.e., the probability distribution for the momentum variables φ, the step size of the leapfrog ε, and the number of leapfrog steps L per iteration. To overcome the challenges related to computation and tuning, the Stan program (Gelman et al. 2014, p. 307; Carpenter et al. 2015, 2017) was developed to automatically apply the HMC. Stan runs HMC using the no-U-turn sampler (NUTS) (Hoffman and Gelman 2014). Al Hakmani and Sheng (2017) used NUTS for the two-parameter mixture IRT (Mix2PL) model and discussed in more detail its performance in estimating model parameters under eight conditions, i.e., two sample sizes per class (250 and 500), two test lengths (20 and 30), and two levels of latent classes (2-class and 3-class). The results indicated that overall, NUTS performs well in retrieving model parameters. Therefore, this research applies the Bayesian method to estimate the parameters of the ZMAR model, using MCMC with the NUTS algorithm, as well as simulation studies to examine different scenarios in order to evaluate whether the proposed mixture model outperforms its counterparts. The models are applied to both the daily IBM common stock closing price from 17 May 1961 to 2 November 1962 (Box et al. 2015, p. 627) and the Brent crude oil price (World Bank 2020). For model selection, we used cross-validation Leave-One-Out (LOO) coupled with the Pareto-smoothed important sampling (PSIS), namely PSIS-LOO. This approach has very efficient computation and was stronger than the Widely Applicable Information Criterion (WAIC) (Vehtari et al. 2017).

The rest of this study is organized as follows. Section 2 describes the definition and properties of Fisher’s z distribution in detail. In Section 3, we introduce the ZMAR model. Section 4 demonstrates the flexibility of the ZMAR model compared with the TMAR and GMAR models using simulated datasets. Section 5 contains the application and comparison of the models using the daily IBM stock prices and the monthly Brent crude oil prices. The conclusion and discussion are given in Section 6.

2. Four-Parameter Fisher’s z Distribution

Let $Y$ be a random variable distributed as an F distribution with $d_1$ and $d_2$degrees of freedom. The density of $Z=\frac{1}{2}\ln Y$ can be defined as

$${\mathsf{\zeta }}_{\left(d_{1k},d_{2k}\right)}\left(z\right)=f_Z\left(z;d_1,d_2\right)=\frac{2{\left(\frac{d_2}{d_1}\right)}^{\frac{1}{2}{d}_2}}{\mathrm{B}\left(\frac{1}{2}{d}_1,\frac{1}{2}{d}_2\right)}\frac{{\mathrm{e}}^{-d_{2}z}}{{\left(1+{\mathrm{e}}^{-2z+\ln \frac{d_2}{d_1}}\right)}^{\left(d_1+d_2\right)/2}},$$

(1)

and the cumulative distribution function (CDF) of $Z$ is expressed as

$${\mathrm{Z}}_{\left(d_{1k},d_{2k}\right)}\left(z\right)={\mathrm{I}}_{z^{*}}\left(\frac{1}{2}{d}_2,\frac{1}{2}{d}_1\right)=\frac{{{\displaystyle \int }}_0^{z^{*}}{t}^{\frac{1}{2}{d}_2-1}{\left(1-t\right)}^{\frac{1}{2}{d}_1-1}dt}{\mathrm{B}\left(\frac{1}{2}{d}_1,\frac{1}{2}{d}_2\right)},$$

(2)

where $\mathrm{e}$ is the exponential constant; $z^{*}=\frac{d_2{\mathrm{e}}^{-2z}}{d_1+d_2{\mathrm{e}}^{-2z}},$ ${\mathrm{I}}_{z^{*}}(.)$ is the incomplete beta function ratio; and $\mathrm{B}(.)$ is the beta function, $-\infty <z<\infty , d_1>0, d_2>0$. Equations (1) and (2) are defined as a probability density function (p.d.f) and a CDF of standardized Fisher’s z distribution, respectively. Let $Z$ be a random variable distributed as a standardized Fisher’s z distribution. Let $\mu$ be a location parameter, and let $\sigma$ be a scale parameter. The density of $X=\sigma Z+\mu$ is (Solikhah et al. 2021)

$$f_X\left(x;d_1,d_2,\mu ,\sigma \right)=\frac{2}{\sigma }\frac{{\left(\frac{d_2}{d_1}\right)}^{\frac{1}{2}{d}_2}}{\mathrm{B}\left(\frac{1}{2}{d}_1,\frac{1}{2}{d}_2\right)}\frac{{\mathrm{e}}^{-d_2\left(\frac{x-\mu }{\sigma }\right)}}{{\left(1+{\mathrm{e}}^{-2\left(\frac{x-\mu }{\sigma }\right)+\ln \frac{d_2}{d_1}}\right)}^{\left(d_1+d_2\right)/2}},$$

(3)

where $-\infty <x<\infty ,-\infty <\mu <\infty ,$and $\sigma >0.$ Equation (3) is defined as a p.d.f of Fisher’s z distribution. It is denoted as $\mathrm{z}\left(d_1,d_2,\mu ,\sigma \right)$. The CDF of the Fisher’s z distribution is expressed as

$$F_X\left(x;d_1,d_2,\mu ,\sigma \right)={\mathrm{I}}_{x^{*}}\left(\frac{1}{2}{d}_2,\frac{1}{2}{d}_1\right)=\frac{1}{\mathrm{B}\left(\frac{1}{2}{d}_1,\frac{1}{2}{d}_2\right)}\underset{0}{\overset{x^{*}}{{\displaystyle \int }}}{t}^{\frac{1}{2}{d}_2-1}{\left(1-t\right)}^{\frac{1}{2}{d}_1-1}dt,$$

(4)

where $x^{*}=\frac{d_2{\mathrm{e}}^{-2\left(\frac{x-\mu }{\sigma }\right)}}{d_1+d_2{\mathrm{e}}^{-2\left(\frac{x-\mu }{\sigma }\right)}}.$ The quantile function (QF) of the Fisher’s z distribution is defined as

$$x_p=\mu +\frac{\sigma }{2}\ln \frac{d_2{\mathrm{I}}^{-1}{}_{x_p}\left(\frac{1}{2}{d}_1,\frac{1}{2}{d}_2\right)}{d_1\left(1-{\mathrm{I}}^{-1}{}_{x_p}\left(\frac{1}{2}{d}_1,\frac{1}{2}{d}_2\right)\right)},$$

(5)

where $\left(d_2{\mathrm{I}}^{-1}{}_{x_p}\left(\frac{1}{2}{d}_1,\frac{1}{2}{d}_2\right)\right)/\left(d_1\left(1-{\mathrm{I}}^{-1}{}_{x_p}\left(\frac{1}{2}{d}_1,\frac{1}{2}{d}_2\right)\right)\right)$ is the QF of the F-distribution and ${\mathrm{I}}^{-1}{}_{x_p}(.)$ is the inversion of the incomplete beta function ratio. Let ${\mathcal{P}}^{-1}{}_{x_p}(.)$ be the inversion of the incomplete gamma function ratio. The QF of the Fisher’s z distribution can be expressed as

$$x_p=\mu +\frac{\sigma }{2}\ln \left(\frac{d_2\left({\mathcal{P}}^{-1}{}_{v_{1p}}\left(d_1/2\right)\right)}{d_1\left({\mathcal{P}}^{-1}{}_{v_{2p}}\left(d_2/2\right)\right)}\right) ,$$

(6)

where $2 {\mathcal{P}}^{-1}{}_{v_{1p}}\left(\frac{1}{2}{d}_1\right)$ and $2 {\mathcal{P}}^{-1}{}_{v_{2p}}\left(\frac{1}{2}{d}_2\right)$ are the QF of the chi-square distribution with $d_1$ and $d_2$ degrees of freedom, respectively. The proofs of Equation (1) up to Equation (6) are postponed to Appendix A. The parameters $d_1$ and $d_2$, known as the shape parameters, are defined for both skewness (symmetrical if $d_1=d_2$, asymmetrical if $d_1\ne d_2$) and fatness of the tails (large $d_1$ and $d_2$ imply thin tails). The Fisher’s z distribution is also always unimodal and has the mode at $x=\mu .$ Furthermore, a change in the value of the parameter $\mu$ only affects the mean of the distribution. It does not affect the variance, skewness, and kurtosis of the distribution. The detailed properties of the Fisher’s z distribution are shown in Appendix B.

A useful tutorial on adding custom function to Stan is provided by Stan Development Team (2018) and Annis et al. (2017). To add a user-defined function, it is first necessary to define a block of function code. The function block must precede all other blocks of Stan code. The code for the random numbers generator function (fisher_z_rng) is shown in Appendix C.1, and the log probability function (fisher_z_lpdf) is shown in Appendix C.2. As an illustration, the p.d.f and CDF of the Fisher’s z distribution with various parameter settings can be seen in Figure 1 and Figure 2, respectively.

3. Fisher’s z Mixture Autoregressive Model

3.1. Model Specification

Let $y_t; t=1,2,\dots ,T$be the real-valued time series of interest; let ${ℱ}_{\left(t-1\right)}$ denote the information set up to time $t-1$, let $F\left(y_t|{ℱ}_{\left(t-1\right)}\right);k=1,2,\cdots ,K$ be the conditional CDF of $Y_t$ given the past information, evaluated at $y_t$; and $K$ is the number of components in the ZMAR model. Let ${\mathrm{Z}}_{\left(d_{1k},d_{2k}\right)}(.)$ be the CDF of the standardized Fisher’s z distribution with $d_{1k}$and $d_{2k}$ shape parameters, given by Equation (2); let ${ϵ}_{k.t}$ be a sequence of independent standardized Fisher’s z random variables such that ${ϵ}_{k.t}$ is independent of $\left\{y_{t-i}, i>0\right\}$; and let ${\sigma }_k$ is a scale parameter of the kth component. The $K$-component ZMAR can be defined as

$$F\left(y_t|{ℱ}_{\left(t-1\right)}\right)={\displaystyle \sum }_{k=1}^K{\eta }_k{\mathrm{Z}}_{\left(d_{1k},d_{2k}\right)}\left(\frac{y_t-{\mu }_{k.t}}{{\sigma }_k}\right),$$

(7)

or

$$y_t=\{\begin{array}{cc}{\mu }_{1.t}+{\sigma }_1{ϵ}_{1.t};& \mathrm{with}\mathrm{probability}{\eta }_1;\\ {\mu }_{2.t}+{\sigma }_2{ϵ}_{2.t};& \mathrm{with}\mathrm{probability}{\eta }_2;\\ ⋮& \\ {\mu }_{K.t}+{\sigma }_K{ϵ}_{K.t};& \mathrm{with}\mathrm{probability}{\eta }_K,\end{array}$$

(8)

with

$${\mu }_{k.t}={\varphi }_{k.0}+{{\displaystyle \sum }}_{i=1}^{p_k}{\varphi }_{k.\mathrm{i}}{y}_{t-\mathrm{i}};k=1,2,\cdots ,K,$$

(9)

where the vector $\mathit{\eta }=\left({\eta }_1, . . . ,{\eta }_K\right)$ is called the weights. $\mathit{\eta }$ takes a value in the unit simplex ${ℰ}_K$, which is a subspace of ${\left({\Re }^{+}\right)}^K$, defined by the following constraint, ${\eta }_k>0$ and ${\displaystyle \sum }_{k=1}^K{\eta }_k=1.$ We use the abbreviation $\mathrm{ZMAR}\left(K;p_1, p_2,\dots , p_K\right)$ for this model, with the parameter $\mathit{\vartheta }=\left({\mathit{\theta }}_1,{\mathit{\theta }}_2,\dots , {\mathit{\theta }}_K,\mathit{\eta }\right);\mathit{\eta }=\left({\eta }_k\right);{\mathit{\theta }}_k=\left(d_{1k},d_{2k},{\sigma }_k,{\varphi }_{k.0},{\mathit{\varphi }}_k\right);{\mathit{\varphi }}_k=\left({\varphi }_{k.1},{\varphi }_{k.2}, \dots ,{\varphi }_{k.p_k}\right);$ to each value $k\in \left\{1,2, . . . ,K\right\};$ taking values in the parameter space ${\Theta }_K={\Theta }^K\times {ℰ}_K$, ${\varphi }_{k.\mathrm{i}}$ denotes the AR coefficient on the kth component and ith lag; $i=1,2, \dots , p_k$; and $p_k$ denotes the autoregressive order of the kth component. Using the parameters ${\mathit{\theta }}_k$, we first define the $K$ auxiliary Fisher’s z AR($p_k$) processes

$$f_k\left(y_t|{ℱ}_{\left(t-1\right)}\right)={\varphi }_{k.0}+{{\displaystyle \sum }}_{i=1}^{p_k}{\varphi }_{k.\mathrm{i}}{y}_{t-\mathrm{i}}+{\sigma }_k{ϵ}_{k.t};k=1,2,\dots ,K,$$

where the AR coefficients ${\mathit{\varphi }}_k$ are assumed to satisfy

$$1-{{\displaystyle \sum }}_{i=1}^{p_k}{\varphi }_{k.\mathrm{i}}{\mathbb{C}}^i\ne 0\mathrm{for}\left|\mathbb{C}\right|\le 1;k=1,\dots ,K.$$

This condition implies that the processes $f_k\left(y_t|{ℱ}_{\left(t-1\right)}\right)$ are stationary and that each component model in (7) or (8) satisfies the usual stationarity condition of the linear AR($p_k$) model.

Suppose $p=\max \left(p_1,p_2, . . . , p_K\right)$. Let a univariate time series $\mathbf{y}=\left(y_1,y_2,\cdots ,y_T\right)$ be influenced by a hidden discrete indicator variables $\mathit{Q}=\left(Q_1,Q_2,\cdots ,Q_T\right)$, where $Q_t$ takes values in the set $\left\{1,2,\dots ,K\right\}$. The probability of sampling from the group labeled $Q_t=k$ is equal to ${\eta }_k$. Suppose that the conditional density of $Y_t$ given $Q_t=k$ is $f_k\left(y_t|{ℱ}_{\left(t-1\right)}\right);\left(k=1,2, \dots , K\right).$Let $\mathit{H}=\left({\mathit{H}}_1{}^{\prime },{\mathit{H}}_2{}^{\prime }, . . . , {\mathit{H}}_T{}^{\prime }\right){}^{\prime }$ be the unobserved random variable, where ${\mathit{H}}_t$ is a $K$-dimensional vector with the values $\mathit{h}=\left({\mathit{h}}_1{}^{\prime },{\mathit{h}}_2{}^{\prime }, . . . , {\mathit{h}}_T{}^{\prime }\right){}^{\prime }$, where $h_{k.t}=1$ if $Q_t=k$ and $h_{k.t}=0$, otherwise. Thus, ${\mathit{H}}_t$ is distributed according to a multinomial distribution consisting of one draw on $K$ categories with probabilities ${\eta }_1, {\eta }_2,\dots , {\eta }_K$ (McLachlan and Peel 2000, p. 7); that is,

$$\mathrm{P}\left({\mathit{H}}_t={\mathit{h}}_t\right)={\eta }_1{}^{h_{1.t}}{\eta }_2{}^{h_{2.t}}\dots {\eta }_K{}^{h_{K.t}}.$$

(10)

Otherwise, this can be written as ${\mathit{H}}_t~{\mathrm{Mult}}_K\left(1,\mathit{\eta }\right)$; where $\mathit{\eta }=\left({\eta }_1, {\eta }_2, . . . ,{\eta }_K\right)$. The conditional likelihood for the $\mathrm{ZMAR}\left(K; p_1, p_2,\dots , p_K\right)$ can be formed by

$$\mathrm{P}\left(y_t|\mathit{\vartheta }\right)={\displaystyle \prod }_{k=1}^K{\displaystyle \prod }_{t=p+1}^T{\left[{\eta }_k{\mathsf{\Delta }}_1{\mathsf{\Delta }}_2\right]}^{h_{k.t}},$$

(11)

where ${\mathsf{\Delta }}_1=\frac{2{\left(d_{2k}/d_{1k}\right)}^{\frac{1}{2}{d}_{2k}}}{{\sigma }_k\mathrm{B}\left(\frac{1}{2}{d}_{1k},\frac{1}{2}{d}_{2k}\right)},$ ${\mathsf{\Delta }}_2=\frac{\exp \left(-d_{2k}\left(a_{k.t}/{\sigma }_k\right)\right)}{{\left[1+\exp \left(-2\left(a_{k.t}/{\sigma }_k\right)+\ln \left[d_{2k}/d_{1k}\right]\right)\right]}^{\left(d_{1k}+d_{2k}\right)/2}},$ and $a_{k.t}=y_t-{\mu }_{k.\mathrm{t}}$.

Let ${\tau }_{t.k}$ be the probability for each t-th observation, $t=\left(p+1\right),\left(p+2\right),\dots ,T$, as members of the k-th component, $k=1,2,\dots ,K$ of a mixture distribution. Suppose ${\mathit{y}}^t=\left(y_1,y_2,\cdots ,y_t\right)$, Bayes’ rule to compute the ${\tau }_{t.k}$ can be expressed as (Frühwirth-Schnatter 2006, p. 26)

$$\mathrm{P}\left(Q_t=k|{\mathit{y}}^t,\mathit{\vartheta }\right)={\tau }_{t.k}=\frac{{\eta }_k{f}_k\left(y_t|{ℱ}_{\left(t-1\right)}\right)}{{{\displaystyle \sum }}_{k=1}^K{\eta }_k{f}_k\left(y_t|{ℱ}_{\left(t-1\right)}\right)}.$$

(12)

Let ${\mathsf{\zeta }}_{\left(d_{1k},d_{2k}\right)}(.)$ be the p.d.f of the standardized Fisher’s z distribution and given by Equation (1). Then, the mixing weights for the ZMAR model can be expressed as

$${\tau }_{t.k}=\frac{{\eta }_k{\sigma }_k{}^{-1}{\mathsf{\zeta }}_{\left(d_{1k},d_{2k}\right)}\left(a_{k.t}/{\sigma }_k\right)}{{{\displaystyle \sum }}_{k=1}^K{\eta }_k{\sigma }_k{}^{-1}{\mathsf{\zeta }}_{\left(d_{1k},d_{2k}\right)}\left(a_{k.t}/{\sigma }_k\right)};k=1,2,\dots ,K.$$

Suppose ${\omega }_{k.t}$ and ${\xi }_k{}^2$ signify the conditional mean and conditional variance of the kth component, which are defined by

$${\omega }_{k.t}=\mathrm{E}\left[Y_t| {\mathit{\theta }}_k\right]={\mu }_{k.t}+\frac{{\sigma }_k}{2}\left(\ln \frac{d_{2k}}{d_{1k}}-\psi \left(\frac{1}{2}{d}_{2k}\right)+\psi \left(\frac{1}{2}{d}_{1k}\right)\right),$$

(13)

and

$${\xi }_k{}^2=\mathrm{Var}\left[Y_t| {\mathit{\theta }}_k\right]={\left(\frac{{\sigma }_k}{2}\right)}^2\left({\psi }^{\prime }\left(\frac{1}{2}{d}_{2k}\right)+{\psi }^{\prime }\left(\frac{1}{2}{d}_{1k}\right)\right),$$

(14)

where $\psi (.)$ is the digamma function and ${\psi }^{\prime }(.)$ is the trigamma function. The conditional mean of $Y_t$ (Frühwirth-Schnatter 2006, p. 10) is obtained as

$${\omega }_t=\mathrm{E}\left[Y_t|{ℱ}_{\left(t-1\right)}\right]={\displaystyle \sum }_{k=1}^K{\eta }_k{\omega }_{k.t},$$

(15)

the conditional variance of $Y_t$ (Frühwirth-Schnatter 2006, p. 11) is obtained as

$${\xi }^2=\mathrm{Var}\left[Y_t|{ℱ}_{\left(t-1\right)}\right]={\displaystyle \sum }_{k=1}^K{\eta }_k\left({\omega }_{k.t}{}^2+{\xi }_k{}^2\right)-{\omega }_t{}^2,$$

(16)

and the higher-order moments around the mean of $Y_t$ (Frühwirth-Schnatter 2006, p. 11) are obtained as

$$\mathrm{E}\left({\left(Y-{\omega }_t\right)}^3|{ℱ}_{\left(t-1\right)}\right)={\displaystyle \sum }_{k=1}^K\left({\left({\omega }_{k.t}-{\omega }_t\right)}^2+3{\xi }_k{}^2\right)\left({\omega }_{k.t}-{\omega }_t\right){\eta }_k,$$

(17)

$$\mathrm{E}\left({\left(Y-{\omega }_t\right)}^4|{ℱ}_{\left(t-1\right)}\right)={\displaystyle \sum }_{k=1}^K\left({\left({\omega }_{k.t}-{\omega }_t\right)}^4+6{\left({\omega }_{k.t}-{\omega }_t\right)}^2{\xi }_k{}^2+3{\xi }_k{}^4\right){\eta }_k.$$

(18)

These expressions apply to any specification of the mixing weights ${\eta }_k$.

3.2. Bayesian Approach for ZMAR Model

In this paper, we apply a Bayesian method to estimate the parameters $\mathit{\vartheta }$. The Bayesian analysis requires the joint posterior density $\pi \left(\mathit{\vartheta }|\mathit{y}\right)$, which is defined by

$$\pi \left(\mathit{\vartheta }|\mathit{y}\right)\propto \mathrm{P}\left(y_t|\mathit{\vartheta }\right)\pi \left(\mathit{\vartheta }\right),$$

(19)

where $\mathrm{P}\left(y_t|\mathit{\vartheta }\right)$ is the conditional likelihood function given by Equation (11) and $\pi \left(\mathit{\vartheta }\right)$ is the prior of the parameter model, which is defined by

$$\pi \left(\mathit{\vartheta }\right)=\pi \left(\mathit{\eta }\right){\displaystyle \prod }_{k=1}^K\pi \left(d_{1k}\right)\pi \left(d_{2k}\right)\pi \left({\sigma }_k\right)\pi \left({\varphi }_{k.0}\right){\displaystyle \prod }_{i=1}^{p_k}\pi \left({\varphi }_{k.i}\right),$$

(20)

where the $\pi \left(\mathit{\eta }\right)$ is prior for the $\mathit{\eta }$ parameter, the $\pi \left(d_{1k}\right), \pi \left(d_{2k}\right), \pi \left({\sigma }_k\right), \pi \left({\varphi }_{k.0}\right)$ and $\pi \left({\varphi }_{k.i}\right)$are prior for the $d_1, d_2, \sigma , {\varphi }_0$ and ${\varphi }_i$ parameters at the kth component; index $i$ denotes the ith lag; $i=1,2,\dots ,p_k$;$k=1,2,\dots ,K$.

Various noninformative prior distributions have been suggested for the prior of AR coefficients, scale parameters, and the selection probabilities in similar models. Huerta and West (1999) analyzed and used the uniform Dirichlet distribution as the prior distribution for the latent variables related to the latent components of an autoregressive model. Gelman (2006) suggested working within the half-t family of prior distributions for variance parameters in the hierarchical modelling, which are more flexible and have better behavior near 0, compared to the inverse-gamma family. Albert and Chib (1993) used the normal distribution as the prior for the autoregressive coefficient in the Markov switching autoregressive model. Based on the findings of the previous studies, we take the singly truncated Student t distribution (positive values only) (Kim 2008) for the priors of the $d_{1k}, d_{2k}$, and ${\sigma }_k$, with the degrees of freedoms $v_{1k}, v_{2k}$,${\nu }_{3k}$, the location parameters $m_{1k},m_{2k},m_{3k}$, and the scale parameters $s_{1k}{}^2, s_{2k}{}^2,s_{3k}{}^2,$ respectively. Therefore, it can be written as $d_{1k}~ t_{{\nu }_{1k}}\left(m_{1k},s_{1k}{}^2\right)I\left(0,\infty \right),$ $d_{2k}~ t_{{\nu }_{2k}}\left(m_{2k},s_{2k}{}^2\right)I\left(0,\infty \right),$and ${\sigma }_k~ t_{{\nu }_{3k}}\left(m_{3k},s_{3k}{}^2\right)I\left(0,\infty \right)$. We take the Dirichlet distribution (Kotz et al. 2000, p. 485) for the prior of the $\mathit{\eta }$ parameter, thus ${\eta }_1,{\eta }_2,\dots ,{\eta }_{K-1} ~ \mathrm{Dir}\left({\delta }_1,{\delta }_2,\dots ,{\delta }_K\right)$. For the priors of the ${\varphi }_{k.0}$ and ${\varphi }_{k.i}$, we take the normal distribution with the location parameters $u_{k.0}, u_{k.i}$ and the scale parameters $g_{k.0}{}^2, g_{k.i}{}^2,$ thus ${\varphi }_{k.0}~N\left(u_{k.0},g_{k.0}{}^2\right)$ and ${\varphi }_{k.i}~N\left(u_{k.i},g_{k.i}{}^2\right);i=1,2,\dots ,p_k.$ Employing the setup of prior distributions, as shown above, the natural logarithm of the joint posterior distribution of the model is given by

$$\ln \pi \left(\mathit{\vartheta }|\mathit{y}\right)\propto {\displaystyle \sum }_{k=1}^K\left(\left({\displaystyle \sum }_{t=p+1}^T{h}_{k.t}\left(\ln {\eta }_k+{\mathsf{\Lambda }}_1+{\mathsf{\Lambda }}_2\right)\right)+{\mathsf{\Lambda }}_3+{\mathsf{\Lambda }}_4+{\mathsf{\Lambda }}_5+{\mathsf{\Lambda }}_6+{\mathsf{\Lambda }}_7+{\displaystyle \sum }_{i=1}^{p_k}{\mathsf{\Lambda }}_8\right),$$

(21)

where ${\mathsf{\Lambda }}_1=d_{2k}\ln \left[d_{2k}/d_{1k}\right]+\ln \left[\mathsf{\Gamma }\left(\frac{1}{2}\left(d_{1k}+d_{2k}\right)\right)\right]-\ln \left[\mathsf{\Gamma }\left(\frac{1}{2}{d}_{1k}\right)\right]-\ln \left[\mathsf{\Gamma }\left(\frac{1}{2}{d}_{2k}\right)\right]-\ln {\sigma }_k,$ ${\mathsf{\Lambda }}_2=-d_{2k}\left(\frac{a_{k.t}}{{\sigma }_k}\right)-\frac{d_{1k}+d_{2k}}{2}\ln \left[1+\text{exp}\left[-2\frac{a_{k.t}}{{\sigma }_k}+\ln \left[d_{2k}/d_{1k}\right]\right]\right]$, ${\mathsf{\Lambda }}_3=\left({\delta }_K-1\right)\ln \left[1-{\displaystyle \sum }_{k=1}^{K-1}{\eta }_k\right]+{\displaystyle \sum }_{k=1}^{K-1}\left({\delta }_k-1\right)\ln {\eta }_k$, ${\mathsf{\Lambda }}_4=-\frac{1}{2}\left({\nu }_{1k}+1\right)\ln \left[1+\frac{1}{{\nu }_{1k}}{\left(\frac{d_{1k}-m_{1k}}{s_{1k}}\right)}^2\right]$, ${\text{ Λ}}_5=-\frac{1}{2}\left({\nu }_{2k}+1\right)\ln \left[1+\frac{1}{{\nu }_{2k}}{\left(\frac{d_{2k}-m_{2k}}{s_{2k}}\right)}^2\right]$, ${\mathsf{\Lambda }}_6=-\frac{1}{2}\left({\nu }_{3k}+1\right)\ln \left[1+\frac{1}{{\nu }_{3k}}{\left(\frac{{\sigma }_k-m_{3k}}{s_{3k}}\right)}^2\right]$, ${\mathsf{\Lambda }}_7=-\frac{1}{2}{\left(\frac{{\varphi }_{k.0}-u_{k.0}}{g_{k.0}}\right)}^2$, and ${\mathsf{\Lambda }}_8=-\frac{1}{2}{\left(\frac{{\varphi }_{k.\mathrm{i}}-u_{k.i}}{g_{k.i}}\right)}^2.$

HMC requires the gradient of the ln-posterior density. In practice, the gradient must be computed analytically (Gelman et al. 2014, p. 301). The gradient of $\ell \left(\mathit{\vartheta }|\mathit{y}\right)$ is

$$\frac{\partial }{\partial \mathit{\vartheta }}\ell \left(\mathit{\vartheta }|\mathit{y}\right)=\frac{\partial }{\partial \mathit{\vartheta }}\ln \pi \left(\mathit{\vartheta }|\mathit{y}\right)$$

$$\frac{\partial }{\partial \mathit{\vartheta }}\ell \left(\mathit{\vartheta }|\mathit{y}\right)=\left(\frac{\partial \ell \left(\mathit{\vartheta }|\mathit{y}\right)}{\partial {\eta }_k},\frac{\partial \ell \left(\mathit{\vartheta }|\mathit{y}\right)}{\partial d_{1k}},\frac{\partial \ell \left(\mathit{\vartheta }|\mathit{y}\right)}{\partial d_{2k}},\frac{\partial \ell \left(\mathit{\vartheta }|\mathit{y}\right)}{\partial {\sigma }_k},\frac{\partial \ell \left(\mathit{\vartheta }|\mathit{y}\right)}{\partial {\varphi }_{k.0}},\frac{\partial \ell \left(\mathit{\vartheta }|\mathit{y}\right)}{\partial {\varphi }_{k.i}}\right)$$

(22)

where $i=1,2,\dots ,p_k;k=1,2,\dots ,K.$ The HMC algorithm for estimating the parameters of the ZMAR model is as follows:

• Determine the initial value of the parameter ${\mathit{\vartheta }}^0$, the diagonal mass matrix $\mathit{M}$, the scale factor of the leapfrog steps $\epsilon$, the number of leapfrog steps $\mathcal{N}$, and the number of iterations $ℛ$.
• For each iteration $r;r=1,2,\dots ,ℛ$ where $ℛ$ represents the number of iterations,
  • Generate the momentum variables $\mathit{\phi }$ with $\mathit{\phi }~\mathrm{Normal}\left(\mathbf{0},\mathit{M}\right)$;
  • For each iteration $n=1,2,\dots ,\mathcal{N}$,

    (1)

    Use the gradient of the ln-posterior density of ${\mathit{\vartheta }}_n$ to make a half-step of ${\mathit{\phi }}_n$

    ${\mathit{\phi }}_{n+0.5}\leftarrow {\mathit{\phi }}_n+\frac{1}{2}\epsilon \frac{\partial \ln \pi \left({\mathit{\vartheta }}_n|\mathit{y}\right)}{\partial {\mathit{\vartheta }}_n},$

    (2)

    Update the vector ${\mathit{\vartheta }}_n$ using the vector ${\mathit{\phi }}_{n+0.5}$

    ${\vartheta }_{\mathit{n}+0.5}\leftarrow {\mathit{\vartheta }}_n+\epsilon {\mathit{M}}^{-1}{\mathit{\phi }}_{n+0.5}$

    (3)

    Update the next half-step for $\mathit{\phi }$

    ${\mathit{\phi }}_{n+1}\leftarrow {\mathit{\phi }}_{n+0.5}+\frac{1}{2}\epsilon \frac{\partial \ln \pi \left({\mathit{\vartheta }}_{n+0.5}|\mathit{y}\right)}{\partial {\mathit{\vartheta }}_{n+0.5}}.$

  • Labels ${\mathit{\vartheta }}^{\left(r-1\right)}$ and ${\mathit{\phi }}^{\left(r-1\right)}$ as the value of the parameter and momentum vectors at the start of the leapfrog process and ${\mathit{\vartheta }}^{\#}, {\mathit{\phi }}^{\#}$ as the value after the $\mathcal{N}$ steps.
  • Compute
    $b=\frac{\pi \left({\mathit{\vartheta }}^{\#}|\mathit{y}\right)f\left({\mathit{\phi }}^{\#}\right)}{\pi \left({\mathit{\vartheta }}^{\left(r-1\right)}|\mathit{y}\right)f\left({\mathit{\phi }}^{\left(r-1\right)}\right)}.$
  • Set ${\mathit{\vartheta }}^{\left(r\right)}$
    ${\mathit{\vartheta }}^{\left(r\right)}=\{\begin{array}{c}{\mathit{\vartheta }}^{\#} \mathrm{with}\mathrm{probability}min\left(b,1\right) \\ {\mathit{\vartheta }}^{\left(r-1\right)}\mathrm{otherwise}.\end{array}$
  • Save ${\mathit{\vartheta }}^{\left(r\right)};$ $r=1,2,\dots ,ℛ$.

The performance of the HMC is very sensitive to two user-defined parameters, i.e., the step size of the leapfrog $\epsilon$ and the number of leapfrog steps $L$. The No-U-Turn Sampler (NUTS) could eliminate the need to set the parameter $L$ and could adapt the step size parameter $\epsilon$ on the fly based on a primal-dual averaging (Hoffman and Gelman 2014). The NUTS algorithm was implemented in C ++ as part of the open-source Bayesian inference package, Stan (Gelman et al. 2014, p. 304; Carpenter et al. 2017). Stan is also a platform for computing log densities and their gradients, so that the densities and gradients are easy to obtain (Carpenter et al. 2015, 2017). Stan can be called from R using the rstan package. An example of the Stan code to fit the ZMAR model can be seen in Section 4.

4. Simulation Studies

A simulation study was carried out to evaluate the performance of the ZMAR model compared to the TMAR and GMAR models. We consider the simulations to accommodate eight scenarios for the conditional density in the first component of the ZMAR model. Furthermore, the conditional densities in the second and third components are specified as the symmetric-fat-tail of the Fisher’s z distributions.

We conducted a Bayesian analysis on the eight simulated datasets, whose datasets were generated by the following steps:

• Step 1: Specify the ZMAR model with three components as $\mathrm{ZMAR}\left(3;1,1,1\right)$ where $f\left(y_t|{ℱ}_{\left(t-1\right)}\right)={\displaystyle \sum }_{k=1}^3{\eta }_k{f}_k\left(y_t|{ℱ}_{\left(t-1\right)}\right)={\displaystyle \sum }_{k=1}^3{\eta }_k\left({\varphi }_{k.0}+{\varphi }_{k.1}{y}_{t-1}+e_{k.t}\right);{\eta }_1={\eta }_2={\eta }_3=\frac{1}{3};{\varphi }_{1.1}=-0.6; {\varphi }_{2.1}=0.2;{\varphi }_{3.1}=0.7$ with $e_{k.t}~\mathrm{z}\left(d_{1.k},d_{2.k},0,{\sigma }_k\right);k=1,2,$3 are the innovations in the first, second, and third components. The scenarios for the simulation are as follows:

  ○

  Scenario 1: represent a mixture of a highly skewed (Bulmer 1967, p. 63) to the left and two symmetrical distributions, where ${\varphi }_{1.0}={\varphi }_{2.0}={\varphi }_{3.0}=0$, excess unconditional kurtosis is 5.73, $e_{1.t}~\mathrm{z}\left(0.2, 10, 0, 5\right),$ $e_{2.t}~\mathrm{z}\left(1, 1, 0, 8\right),$ and $e_{3.t}~\mathrm{z}\left(30, 30, 0, 10\right);$

  ○

  Scenario 2: represent a mixture of a highly skewed to the right and two symmetrical distributions, where ${\varphi }_{1.0}={\varphi }_{2.0}={\varphi }_{3.0}=0,$ excess unconditional kurtosis is 2.23, $e_{1.t}~\mathrm{z}\left(20, 1, 0, 5\right),$ $e_{2.t}~\mathrm{z}\left(1, 1, 0, 8\right),$ and $e_{3.t}~\mathrm{z}\left(30, 30, 0, 10\right);$

  ○

  Scenario 3: represent a mixture of three symmetrical distributions, where ${\varphi }_{1.0}={\varphi }_{2.0}={\varphi }_{3.0}=0,$ excess unconditional kurtosis is 1.67, $e_{1.t}~\mathrm{z}\left(0.5, 0.5, 0, 5\right)$, $e_{2.t}~\mathrm{z}\left(1, 1, 0, 8\right),$ and $e_{3.t}~\mathrm{z}\left(30, 30, 0, 10\right)$;

  ○

  Scenario 4: represent a mixture of moderately skewed (Bulmer 1967, p. 63) distributions, where ${\varphi }_{1.0}={\varphi }_{2.0}={\varphi }_{3.0}=0,$ excess unconditional kurtosis is 2.12, $e_{1.t}~\mathrm{z}\left(3, 10, 0, 8\right),$ $e_{2.t}~\mathrm{z}\left(1, 1, 0, 5\right),$ and $e_{3.t}~\mathrm{z}\left(30, 30, 0, 10\right);$

  ○

  Scenario 5: represent a mixture of three fairly symmetrical distributions (Bulmer 1967, p. 63), where ${\varphi }_{1.0}={\varphi }_{2.0}={\varphi }_{3.0}=0,$ excess unconditional kurtosis is 4.25, $e_{1.t}~\mathrm{z}\left(7, 10, 0, 8\right),$ $e_{2.t}~\mathrm{z}\left(1, 1, 0, 5\right),$ and $e_{3.t}~\mathrm{z}\left(30, 30, 0, 10\right);$

  ○

  Scenario 6: represent a mixture of three symmetrical distributions, where ${\varphi }_{1.0}={\varphi }_{2.0}={\varphi }_{3.0}=0$, excess unconditional kurtosis is 4.47, $e_{1.t}~\mathrm{z}\left(10, 10, 0, 8\right)$, $e_{2.t}~\mathrm{z}\left(1, 1, 0, 5\right),$ and $e_{3.t}~\mathrm{z}\left(30, 30, 0, 10\right)$;

  ○

  Scenario 7: represent a mixture of three symmetrical distributions, where ${\varphi }_{1.0}=-1,$ ${\varphi }_{2.0}={\varphi }_{3.0}=0,$ excess unconditional kurtosis is 1.63, $e_{1.t}~\mathrm{z}\left(0.5, 0.5, 0, 5\right)$, $e_{2.t}~\mathrm{z}\left(1, 1, 0, 8\right),$ and $e_{3.t}~\mathrm{z}\left(30, 30, 0, 10\right)$;

  ○

  Scenario 8: represent a mixture of three symmetrical distributions, where ${\varphi }_{1.0}=-20$, ${\varphi }_{2.0}={\varphi }_{3.0}=0,$ excess unconditional kurtosis is 0.67, $e_{1.t}~\mathrm{z}\left(0.5, 0.5, 0, 5\right)$, $e_{2.t}~\mathrm{z}\left(1, 1, 0, 8\right),$ and $e_{3.t}~\mathrm{z}\left(30, 30, 0, 10\right)$;

The innovations of Scenario 7 and Scenario 8 are the same as those of Scenario 3. The comparison of graph visualizations for the innovations in the first component in Scenario 1 to Scenario 6 represented specifically for the Fisher’s z distribution can be seen in Figure 1 and Figure 2;

• Step 2: Generate $e_{1.t},e_{2.t},$ and $e_{3.t}$;
• Step 3: Generate ${\mathit{H}}_t~{\mathrm{Mult}}_3\left(1,\mathit{\eta }\right)$; where $\mathit{\eta }=\left(\frac{1}{3}, \frac{1}{3}, \frac{1}{3}\right)$;
• Step 4: Compute $f_k\left(y_t|{ℱ}_{\left(t-1\right)}\right)={\varphi }_{k.0}+{\varphi }_{k.1}{y}_{t-1}+e_{k.t}$; $k=1,2,$3;
• Step 5: Compute $y_t=f_1{\left(y_t|{ℱ}_{\left(t-1\right)}\right)}^{h_{1.t}}{f}_2{\left(y_t|{ℱ}_{\left(t-1\right)}\right)}^{h_{2.t}}{f}_3{\left(y_t|{ℱ}_{\left(t-1\right)}\right)}^{h_{3.t}}$.

Furthermore, we generated 600 datasets and fit all of the models for each simulated dataset, respectively, to find the best performance in terms of model comparisons. Figure 3 shows the simulated datasets for Scenario 1 to Scenario 8. We implemented the models using the rstan package (Stan Development Team 2020), the R interface to Stan developed by the Stan Development Team in the R software. Suppose $\mathrm{p}k=p_k;\mathrm{con}\_\mathrm{eta}=\mathrm{c}\left(1,1,1\right);\mathrm{tau}\left[\mathrm{t},\mathrm{k}\right]={\tau }_{t.k};\mathrm{eta}\left[\mathrm{k}\right]={\eta }_k;\mathrm{sigma}\left[\mathrm{k}\right]={\sigma }_k;\mathrm{d}1\left[\mathrm{k}\right]=d_{1k};\mathrm{d}2\left[\mathrm{k}\right]=d_{2k};\mathrm{phi}k\left[1\right]={\varphi }_{k.1};k=1,2,3.$ Here is an example of the Stan code to fit the ZMAR model in Scenario 1:

fitMAR1 = "

functions{

real fisher_z_lpdf(real x,real d1,real d2,real mu,real sigma){

return (log(2)+0.5*d2*(log(d2)-log(d1))-d2*(x-mu)/sigma-log(sigma)-lbeta(0.5*d1,0.5*d2)-

    (d1+d2)/2*log1p_exp((-2*(x-mu)/sigma)+log(d2)-log(d1)));

}

}

data {

int<lower=0> p1;

int<lower=0> p2;

int<lower=0> p3;

int T;

vector[T] y;

vector[3] con_eta;

}

parameters {

simplex[3] eta; //mixing proportions

vector<lower = 0>[3] sigma;

vector<lower = 0>[3] d1;

vector<lower = 0>[3] d2;

real phi1[p1];

real phi2[p2];

real phi3[p3];

}

model {

matrix[T,3] tau;

real lf[T];

eta ~ dirichlet(con_eta);

//priors

sigma[1] ~ student_t(3,5,0.1);

sigma[2] ~ student_t(3,8,0.1);

sigma[3] ~ student_t(3,10,0.1);

d1[1] ~ student_t(3,0.2,0.1);

d1[2] ~ student_t(3,1,0.1);

d1[3] ~ student_t(3,30,0.1);

d2[1] ~ student_t(3,10,0.1);

d2[2] ~ student_t(3,1,0.1);

d2[3] ~ student_t(3,30,0.1);

phi1[1] ~ normal(-0.6,0.1);

phi2[1] ~ normal(0.2,0.1);

phi3[1] ~ normal(0.7,0.1);

//ZMAR model

for(t in 1:T) {

if(t==1) {

tau[t,1] = log(eta[1])+fisher_z_lpdf(y[t]/sigma[1]|d1[1], d2[1],0,1)-log(sigma[1]);

tau[t,2] = log(eta[2])+fisher_z_lpdf(y[t]/sigma[2]|d1[2], d2[2],0,1)-log(sigma[2]);

tau[t,3] = log(eta[3])+fisher_z_lpdf(y[t]/sigma[3]|d1[3], d2[3],0,1)-log(sigma[3]);

} else {

real mu1 = 0;

real mu2 = 0;

real mu3 = 0;

for (i in 1:p1)

mu1 += phi1[i] * y[t-i];

for (i in 1:p2)

mu2 += phi2[i] * y[t-i];

for (i in 1:p3)

mu3 += phi3[i] * y[t-i];

tau[t,1] = log(eta[1])+fisher_z_lpdf((y[t]-mu1)/sigma[1]|d1[1],d2[1],0,1)-log(sigma[1]);

tau[t,2] = log(eta[2])+fisher_z_lpdf((y[t]-mu2)/sigma[2]|d1[2],d2[2],0,1)-log(sigma[2]);

tau[t,3] = log(eta[3])+fisher_z_lpdf((y[t]-mu3)/sigma[3]|d1[3],d2[3],0,1)-log(sigma[3]);

}

lf[t] = log_sum_exp(tau[t,]);

}

target += sum(lf);

} "

The warm-up stage in these simulation studies was set to 1500 iterations, 3 chains with 5000 sampling iterations, and 1 thin. The adapt_delta parameter was set to 0.99, and the max_treedepth was set to 15. For all scenarios, the parameter priors of the ZMAR, TMAR, and GMAR models are shown in Appendix D, and their posterior inferences are presented in Appendix E. There are a variety of convergence diagnoses, such as the potential scale reduction factor $\hat{R}$ (Gelman and Rubin 1992; Susanto et al. 2018; Gelman et al. 2014, p. 285; Vehtari et al. 2020) and the effective sample size n_eff (Gelman et al. 2014, p. 266; Vehtari et al. 2020). If the MCMC chain has reached convergence, the $\hat{R}$ statistic is less than 1.01, and the n_eff statistic is greater than 400 (Vehtari et al. 2020). To compare the performance of the models, we use the PSIS-LOO.

Table 1. Model comparison using pareto-smoothed important sampling leave-one-out cross-validation (PSIS-LOO) in several scenarios.

Scenarios	ZMAR		TMAR		GMAR
	Model	PSIS-LOO	Model	PSIS-LOO	Model	PSIS-LOO
1	ZMAR(3;1,1,1)	4563.6 *	TMAR(3;1,1,1)	4639.1	GMAR(4;1,1,1,1)	4646.6
2	ZMAR(3;1,1,1)	4120.0 *	TMAR(3;1,1,1)	4139.0	GMAR(3;1,1,1)	4130.8
3	ZMAR(3;1,1,1)	4483.3	TMAR(3;1,1,1)	4483.1	GMAR(3;1,1,1)	4481.1 *
4	ZMAR(3;1,1,1)	3639.7 *	TMAR(3;1,1,1)	3652.7	GMAR(3;1,1,1)	3660.9
5	ZMAR(3;1,1,1)	3523.9 *	TMAR(3;1,1,1)	3534.1	GMAR(4;1,1,1,1)	3539.5
6	ZMAR(3;1,1,1)	3490.1 *	TMAR(3;1,1,1)	3496.2	GMAR(3;1,1,1)	3494.6
7	ZMAR(3;1,1,1)	4483.6	TMAR(3;1,1,1)	4481.5	GMAR(3;1,1,1)	4478.0 *
8	ZMAR(3;1,1,1)	4549.1 *	TMAR(3;1,1,1)	4579.4	GMAR(3;1,1,1)	4575.1

Simulation Model: $f\left(y_t|{ℱ}_{\left(t-1\right)}\right)=\frac{1}{3}\left({\varphi }_{1.0}+-0.6y_{t-1}+e_{1.t}\right)+\frac{1}{3}\left({\varphi }_{2.0}+0.2y_{t-1}+e_{2.t}\right)+\frac{1}{3}\left({\varphi }_{3.0}+0.7y_{t-1}+e_{3.t}\right);$ Scenario 1: ${\varphi }_{1.0}={\varphi }_{2.0}={\varphi }_{3.0}=0$, $e_{1.t}~\mathrm{z}\left(0.2,10,0,5\right),$ $e_{2.t}~\mathrm{z}\left(1,1,0,8\right),e_{3.t}~\mathrm{z}\left(30,30,0,10\right);$ Scenario 2: ${\varphi }_{1.0}={\varphi }_{2.0}={\varphi }_{3.0}=0$, $e_{1.t}~\mathrm{z}\left(20,1,0,5\right),$ $e_{2.t}~\mathrm{z}\left(1,1,0,8\right),$ $e_{3.t}~\mathrm{z}\left(30,30,0,10\right)$; Scenario 3: ${\varphi }_{1.0}={\varphi }_{2.0}={\varphi }_{3.0}=0$, $e_{1.t}~\mathrm{z}\left(0.5,0.5,0,5\right),$ $e_{2.t}~\mathrm{z}\left(1,1,0,8\right),$ $e_{3.t}~\mathrm{z}\left(30,30,0,10\right)$; Scenario 4: ${\varphi }_{1.0}={\varphi }_{2.0}={\varphi }_{3.0}=0$, $e_{1.t}~\mathrm{z}\left(3,10,0,8\right),$ $e_{2.t}~\mathrm{z}\left(1,1,0,5\right),$ $e_{3.t}~\mathrm{z}\left(30,30,0,10\right)$; Scenario 5: ${\varphi }_{1.0}={\varphi }_{2.0}={\varphi }_{3.0}=0$, $e_{1.t}~\mathrm{z}\left(7,10,0,8\right),$ $e_{2.t}~\mathrm{z}\left(1,1,0,5\right),e_{3.t}~\mathrm{z}\left(30,30,0,10\right)$; Scenario 6: ${\varphi }_{1.0}={\varphi }_{2.0}={\varphi }_{3.0}=0$, $e_{1.t}~\mathrm{z}\left(10,10,0,8\right),$ $e_{2.t}~\mathrm{z}\left(1,1,0,5\right),e_{3.t}~\mathrm{z}\left(30,30,0,10\right)$; Scenario 7: ${\varphi }_{1.0}=-1$, ${\varphi }_{2.0}={\varphi }_{3.0}=0,$ $e_{1.t}~\mathrm{z}\left(0.5,0.5,0,5\right),$ $e_{2.t}~\mathrm{z}\left(1,1,0,8\right),e_{3.t}~\mathrm{z}\left(30,30,0,10\right)$; Scenario 8: ${\varphi }_{1.0}=-20$, ${\varphi }_{2.0}={\varphi }_{3.0}=0, e_{1.t}~\mathrm{z}\left(0.5,0.5,0,5\right),$ $e_{2.t}~\mathrm{z}\left(1,1,0,8\right),$ $e_{3.t}~\mathrm{z}\left(30,30,0,10\right)$; * indicates the smallest value.

shows the summary simulation result for all scenarios, which indicates that the ZMAR model performs the best when the datasets are generated from the ZMAR model, in which one of the components is asymmetric. When all the components are symmetric, the ZMAR model also performs better than TMAR and GMAR, as long as the excess unconditional kurtosis is large enough or the intercept distances between the components are far apart. However, when all of the mixture components are symmetric, the excess unconditional kurtosis is small and the intercept distances between the components are close enough or the intercepts are the same, then the GMAR model plays the best. Let us now focus on the results of Scenario 3 and Scenario 6.

In the third and sixth scenarios, the datasets are generated from three symmetrical distributions. The two scenarios have the same intercepts and are generated with different unconditional kurtosis. Scenario 3 has a smaller unconditional kurtosis than Scenario 6. The best ZMAR, TMAR, and GMAR models for the two scenarios are ZMAR(3;1,1,1), TMAR(3;1,1,1), and GMAR(3;1,1,1), where the values of PSIS-LOO are 4483.30, 4483.10, and 4481.10, for the third scenario and 3490.10, 3496.20, and 3494.60, for the sixth scenario. Clearly, the PSIS-LOO value for the GMAR model is smaller than for the ZMAR and TMAR models in the third scenario, and the ZMAR model has the smallest PSIS-LOO value in the sixth scenario. When the intercepts in Scenario 3 are varied and determined as in Scenario 7, the GMAR model is also the best. However, when the intercepts in Scenario 3 are varied and determined as in Scenario 8, the ZMAR model is the best. For other scenarios, in datasets generated from asymmetric components (the first, second, fourth, and fifth scenarios), the ZMAR model is the best.

5. Application for Real Data

5.1. IBM Stock Prices

To illustrate the potential of the ZMAR model, we consider the daily IBM common stock closing price from 17 May 1961 to 2 November 1962 (Box et al. 2015, p. 627). This time series has been analyzed by many researchers such as Le et al. (1996) and Wong and Li (2000). Wong and Li (2000) used the EM algorithm to estimate the parameters model and has been identified that the best GMAR model for the series was a GMAR(3;1,1,0).

Figure 4 shows that the IBM stock prices series has a trimodal marginal distribution, where the estimated locations of the modes are at 377.48, 481.78, and 547.90 points. Therefore, we choose the three-component ZMAR model for the differenced series. The orders of the autoregressive components are chosen by the minimum PSIS-LOO. The best three-component ZMAR model is ZMAR(3;0,1,1) without intercept and the value of PSIS-LOO is 2424.40. The warm-up stage for the ZMAR, TMAR, and GMAR models was set to 1500 iterations followed by 3 chains with 5000 sampling iterations and 1 thin; the adapt_delta parameter was set to 0.99, and the max_treedepth was set to 15.

Table 2. Summary of posterior inferences for ZMAR, TMAR, and GMAR models in the IBM stock prices (first-differenced series).

Parameters	ZMAR Model					TMAR Model					GMAR Model
	mean	2.5%	97.5%	n_eff	Rhat	mean	2.5%	97.5%	n_eff	Rhat	mean	2.5%	97.5%	n_eff	Rhat
d1[1]	1.95	1.65	2.26	7185	1	-	-	-	-	-	-	-	-	-	-
d1[2]	1.87	1.65	2.21	6655	1	-	-	-	-	-	-	-	-	-	-
d1[3]	4.91	4.58	5.23	6008	1	-	-	-	-	-	-	-	-	-	-
d2[1]	3.90	3.60	4.23	5905	1	-	-	-	-	-	-	-	-	-	-
d2[2]	6.41	6.09	6.75	2798	1	-	-	-	-	-	-	-	-	-	-
d2[3]	1.70	1.50	1.95	8370	1	-	-	-	-	-	-	-	-	-	-
nu[1]	-	-	-	-	-	12.52	12.21	12.84	7072	1	-	-	-	-	-
nu[2]	-	-	-	-	-	10.77	10.46	11.09	6020	1	-	-	-	-	-
nu[3]	-	-	-	-	-	14.03	13.72	14.36	5384	1	-	-	-	-	-
eta[1]	0.01	0.00	0.05	13945	1	0.58	0.44	0.71	6826	1	0.54	0.41	0.67	5941	1
eta[2]	0.46	0.34	0.58	10248	1	0.40	0.27	0.54	6713	1	0.42	0.29	0.55	5853	1
eta[3]	0.53	0.41	0.64	10235	1	0.02	0.00	0.05	10685	1	0.04	0.01	0.09	9810	1
sigma[1]	28.28	27.97	28.58	4083	1	4.97	4.65	5.21	5145	1	4.82	4.56	5.08	7747	1
sigma[2]	9.77	9.40	10.05	2413	1	5.80	5.50	6.07	7715	1	6.01	5.72	6.29	7369	1
sigma[3]	6.34	6.05	6.63	9774	1	25.02	24.69	25.35	7061	1	19.04	18.71	19.36	4525	1
phi1[1]	-	-	-	-	-	−0.29	−0.41	−0.17	8124	1	−0.32	−0.43	−0.21	6681	1
phi2[1]	0.61	0.47	0.76	11354	1	0.68	0.52	0.84	8061	1	0.67	0.52	0.82	7592	1
phi3[1]	−0.28	−0.39	−0.17	10868	1	-	-	-	-	-	-	-	-	-	-
PSIS-LOO	2424.40					2431.30					2431.70

Note: $\mathrm{eta}\left[1\right]={\eta }_1; \mathrm{eta}\left[2\right]={\eta }_2;\mathrm{eta}\left[3\right]={\eta }_3;\mathrm{sigma}\left[1\right]={\sigma }_1; \mathrm{sigma}\left[2\right]={\sigma }_2;\mathrm{sigma}\left[3\right]={\sigma }_3;\mathrm{d}1\left[1\right]=d_{11};\mathrm{d}1\left[2\right]=d_{12};$ $\mathrm{d}1\left[3\right]=d_{13};\mathrm{d}2\left[1\right]=d_{21};\mathrm{d}2\left[2\right]=d_{22};\mathrm{d}2\left[3\right]=d_{23};\mathrm{nu}\left[1\right]={\nu }_1;\mathrm{nu}\left[2\right]={\nu }_2;\mathrm{nu}\left[3\right]={\nu }_3;\mathrm{phi}1\left[1\right]={\varphi }_{1.1};$ $\mathrm{phi}2\left[1\right]={\varphi }_{2.1};\mathrm{phi}3\left[1\right]={\varphi }_{3.1}.$

shows the summary of posterior inferences for all models. The prior distributions and posterior density plots for the parameters of the ZMAR model are presented in Appendix F.1 and Appendix G.1.1, respectively. For all the parameters of the ZMAR model, the MCMC chain has reached convergence, which is shown by the $\widehat{R}$ statistic being less than 1.01 and the $n_{eff}$ statistic being greater than 400. The three-component ZMAR model for the differenced series was then transformed into a three-component ZMAR model for the original series, namely

$$\begin{array}{c}F\left(y_t|{ℱ}_{\left(t-1\right)}\right)=0.01{\mathrm{Z}}_{\left(1.95,3.90\right)}\left(\frac{y_t-y_{t-1}}{28.28}\right)+0.46{\mathrm{Z}}_{\left(1.87,6.41\right)}\left(\frac{y_t-1.61y_{t-1}+0.61y_{t-2}}{9.77}\right)\\ +0.53{\mathrm{Z}}_{\left(4.91,1.70\right)}\left(\frac{y_t-0.72y_{t-1}-0.28y_{t-2}}{6.34}\right).\end{array}$$

(23)

We compared the ZMAR model with the TMAR and GMAR models. The best three-component TMAR and GMAR models are TMAR(3;1,1,0) and GMAR(3;1,1,0), without intercept. The PSIS-LOO values of the TMAR and GMAR models are 2431.30 and 2431.70, respectively. The prior distributions and posterior density plots for the parameters of the TMAR and GMAR models are presented in Appendix F.1 and Appendix G.1. Here is the result of the summary posterior inferences for all models

The MCMC chains for the TMAR parameter model and the GMAR parameter model have also reached convergence, which is shown by the $\widehat{R}$ statistic being less than 1.01 and the $n_{eff}$ statistic being greater than 400. Let $\mathsf{\Phi }(.)$, and let ${\mathcal{T}}_{{\nu }_{\mathrm{k}}}(.)$ be the CDF of the standard normal distribution and the standardized Student t distribution with ${\nu }_k;k=1,2,3$ degrees of freedom. The three-component TMAR model for the original series, namely

$$\begin{array}{ll}F\left(y_t|{ℱ}_{\left(t-1\right)}\right)& =0.58{\mathcal{T}}_{12.52}\left(\frac{y_t-0.71y_{t-1}-0.29y_{t-2}}{4.97}\right)\\ & +0.40{\mathcal{T}}_{10.77}\left(\frac{y_t-1.68y_{t-1}+0.68y_{t-2}}{5.80}\right)+0.02{\mathcal{T}}_{14.03}\left(\frac{y_t-y_{t-1}}{25.02}\right).\end{array}$$

(24)

and a three-component GMAR for the original series, namely

$$\begin{array}{ll}F\left(y_t|{ℱ}_{\left(t-1\right)}\right)& =0.54 \mathsf{\Phi }\left(\frac{y_t-0.68 y_{t-1}-0.32 y_{t-2}}{4.82}\right)\\ & +0.42 \mathsf{\Phi }\left(\frac{y_t-1.67 y_{t-1}+0.67 y_{t-2}}{6.01}\right)+0.04 \mathsf{\Phi }\left(\frac{y_t-y_{t-1}}{19.04}\right).\end{array}$$

(25)

Therefore, the ZMAR model is preferred over the TMAR and GMAR models, which are indicated by the PSIS-LOO value of the ZMAR model being the smallest.

5.2. Brent Crude Oil Prices

In December 1988, the Organization of Petroleum Exporting Countries (OPEC) decided to adopt Brent as a new benchmark, rather than the value of the Arabian light (Carollo 2012, p. 10). Since then, Brent has been one of the main benchmarks for oil purchases in the world, the other being West Texas Intermediate (WTI). Figure 5a shows the monthly Brent crude oil price from January 1989 to June 2020 (in U.S. dollars per barrel), taken from the World Bank (2020). Figure 5b shows the first-differenced series. Figure 6 shows the marginal distribution of the original series as being trimodal, where the estimated locations of the modes are at 18.57, 61.64, and 110.15 points. Therefore, we also decided to choose a three-component mixture model for the ZMAR, TMAR, and GMAR models applied to the differenced series. The best three-component mixture model estimates for each of the ZMAR, TMAR, and GMAR models are ZMAR (3;1,2,2), TMAR (3;2,1,3), and GMAR (3;3,3,3), without intercept.

The warm-up steps for all the models were set to 1500 iterations with 5000 sampling iterations and 1 thin, followed by 3 chains; the adapt_delta parameters were set to 0.99, and the max_treedepths were set to 15.

Table 3. Summary of posterior inferences for ZMAR, TMAR, and GMAR models in the Brent monthly crude oil prices (first-differenced series).

Parameters	ZMAR Model					TMAR Model					GMAR Model
	mean	2.5%	97.5%	n_eff	Rhat	mean	2.5%	97.5%	n_eff	Rhat	mean	2.5%	97.5%	n_eff	Rhat
d1[1]	13.22	12.89	13.57	1136	1	-	-	-	-	-	-	-	-	-	-
d1[2]	0.99	0.78	1.16	7606	1	-	-	-	-	-	-	-	-	-	-
d1[3]	10.09	9.76	10.43	783	1	-	-	-	-	-	-	-	-	-	-
d2[1]	4.49	4.19	4.83	6163	1	-	-	-	-	-	-	-	-	-	-
d2[2]	1.91	1.64	2.21	8505	1	-	-	-	-	-	-	-	-	-	-
d2[3]	4.37	4.06	4.80	786	1	-	-	-	-	-	-	-	-	-	-
nu[1]	-	-	-	-	-	14.98	14.64	15.30	5742	1	-	-	-	-	-
nu[2]	-	-	-	-	-	12.08	11.77	12.39	6135	1	-	-	-	-	-
nu[3]	-	-	-	-	-	4.34	4.01	4.68	4090	1	-	-	-	-	-
eta[1]	0.40	0.30	0.51	7598	1	0.34	0.25	0.44	7002	1	0.45	0.35	0.55	8395	1
eta[2]	0.39	0.28	0.50	8658	1	0.40	0.29	0.49	7202	1	0.28	0.18	0.39	7317	1
eta[3]	0.21	0.12	0.31	8287	1	0.26	0.17	0.36	7343	1	0.27	0.18	0.36	9667	1
sigma[1]	5.14	4.78	5.40	5164	1	4.93	4.23	5.23	1419	1	4.18	3.96	4.52	5266	1
sigma[2]	2.82	2.57	3.08	9271	1	1.62	1.39	1.86	7816	1	1.66	1.37	1.88	6374	1
sigma[3]	7.12	6.03	7.46	674	1	1.80	1.42	2.03	4462	1	1.92	1.69	2.15	10687	1
phi1[1]	−0.37	−0.49	−0.26	8682	1	0.61	0.45	0.77	8641	1	0.14	0.01	0.28	7866	1
phi1[2]	-	-	-	-	-	−0.42	−0.56	−0.27	7932	1	−0.35	−0.48	−0.22	8284	1
phi1[3]	-	-	-	-	-	-	-	-	-	-	−0.19	−0.31	−0.07	8919	1
phi2[1]	0.68	0.55	0.81	8187	1	−0.28	−0.42	−0.15	7045	1	0.48	0.36	0.60	7799	1
phi2[2]	−0.34	−0.48	−0.21	8060	1	-	-	-	-	-	−0.17	−0.28	−0.06	9453	1
phi2[3]	-	-	-	-	-	-	-	-	-	-	0.38	0.24	0.54	7424	1
phi3[1]	0.69	0.55	0.84	9548	1	0.54	0.44	0.64	7735	1	0.57	0.47	0.67	6397	1
phi3[2]	0.72	0.57	0.86	9333	1	0.86	0.75	0.97	6974	1	0.91	0.81	1.02	7326	1
phi3[3]	-	-	-	-	-	−0.28	−0.41	−0.14	7427	1	−0.29	−0.41	−0.16	7545	1
PSIS-LOO	2024.50					2034.10					2048.70

Note: $\mathrm{eta}\left[1\right]={\eta }_1; \mathrm{eta}\left[2\right]={\eta }_2;\mathrm{eta}\left[3\right]={\eta }_3;\mathrm{sigma}\left[1\right]={\sigma }_1; \mathrm{sigma}\left[2\right]={\sigma }_2;\mathrm{sigma}\left[3\right]={\sigma }_3;\mathrm{d}1\left[1\right]=d_{11};\mathrm{d}1\left[2\right]=d_{12};$ $\mathrm{d}1\left[3\right]=d_{13};\mathrm{d}2\left[1\right]=d_{21};\mathrm{d}2\left[2\right]=d_{22};\mathrm{d}2\left[3\right]=d_{23};\mathrm{nu}\left[1\right]={\nu }_1;\mathrm{nu}\left[2\right]={\nu }_2;\mathrm{nu}\left[3\right]={\nu }_3;\mathrm{phi}1\left[1\right]={\varphi }_{1.1};\mathrm{phi}1\left[2\right]={\varphi }_{1.2};$ $\mathrm{phi}1\left[3\right]={\varphi }_{1.3};\mathrm{phi}2\left[1\right]={\varphi }_{2.1};\mathrm{phi}2\left[2\right]={\varphi }_{2.2};\mathrm{phi}2\left[3\right]={\varphi }_{2.3};\mathrm{phi}3\left[1\right]={\varphi }_{3.1};\mathrm{phi}3\left[2\right]={\varphi }_{3.2};\mathrm{phi}3\left[3\right]={\varphi }_{3.3}.$

shows the summary of posterior inferences for all models. The prior distributions and posterior density plots for the parameters of the ZMAR model, the TMAR model, and the GMAR model are presented, respectively, in Appendix F.2 and Appendix G.2.

For all the parameter models, the MCMC chains reached convergence, which was shown by the $\widehat{R}$ statistic being less than 1.01 and the $n_{eff}$ statistic being greater than 400.

The three-component ZMAR model for the original series, namely

$$\begin{array}{ll}F\left(y_t|{ℱ}_{\left(t-1\right)}\right)& =0.40{\mathrm{Z}}_{\left(13.22,4.49\right)}\left(\frac{y_t-0.63y_{t-1}-0.37y_{t-2}}{5.14}\right)\\ & +0.39{\mathrm{Z}}_{\left(0.99,1.91\right)}\left(\frac{y_t-1.68y_{t-1}+1.02y_{t-2}-0.34y_{t-3}}{2.82}\right)\\ & +0.21{\mathrm{Z}}_{\left(10.09,4.37\right)}\left(\frac{y_t-1.69y_{t-1}-0.03y_{t-2}+0.72y_{t-3}}{7.12}\right).\end{array}$$

(26)

a three-component TMAR model for the original series, namely

$$\begin{array}{ll}F\left(y_t|{ℱ}_{\left(t-1\right)}\right)& =0.34{\mathcal{T}}_{14.98}\left(\frac{y_t-1.61y_{t-1}+1.03y_{t-2}-0.42y_{t-3}}{4.93}\right)\\ & +0.40{\mathcal{T}}_{12.08}\left(\frac{y_t-0.72y_{t-1}-0.28y_{t-2}}{1.62}\right)\\ & +0.26{\mathcal{T}}_{4.34}\left(\frac{y_t-1.54y_{t-1}-0.32y_{t-2}+1.14y_{t-3}-0.28y_{t-4}}{1.80}\right).\end{array}$$

(27)

and a three-component GMAR for the original series, namely

$$\begin{array}{ll}F\left(y_t|{ℱ}_{\left(t-1\right)}\right)& =0.45 \mathsf{\Phi }\left(\frac{y_t-1.14y_{t-1}+0.49y_{t-2}-0.16y_{t-3}-0.19y_{t-4}}{4.18}\right)\\ & +0.28 \mathsf{\Phi }\left(\frac{y_t-1.48y_{t-1}+0.65y_{t-2}-0.55y_{t-3}+0.38y_{t-4}}{1.66}\right)\\ & +0.27 \mathsf{\Phi }\left(\frac{y_t-1.57y_{t-1}-0.34y_{t-2}+1.20y_{t-3}-0.29y_{t-4}}{1.92}\right).\end{array}$$

(28)

The PSIS-LOO values of the ZMAR, TMAR, and GMAR models are 2024.40, 2034.10, and 2048.70, respectively. Therefore, the ZMAR model is preferred over the TMAR and GMAR models, which is indicated by the PSIS-LOO value of the ZMAR model being the smallest.

6. Conclusions

We have discussed the definition and properties of the four-parameter Fisher’s z distribution. The four-parameters of the Fisher’s z distribution are $\mu ,\sigma , d_1,$ and $d_2.$ The $\mu$ is a location parameter, the $\sigma$ is a scale parameter, and the $d_1$,$d_2$ are known as the shape parameters, defined for both skewness (symmetric if $d_1=d_2$, asymmetric if $d_1\ne d_2$) and fatness of the tails (large $d_1$ and $d_2$ imply thin tails). The Fisher’s z distribution is always unimodal and has the mode at $x=\mu .$ The value of $\mu$ only affects the mean of the distribution. It does not affect the variance, skewness, and kurtosis of the distribution. Furthermore, if $d_1=d_2$, then the mean is equal to $\mu$; if $d_1<d_2$, then the mean is less than $\mu$; and if $d_1>d_2$, then the mean is greater than $\mu$. The excess kurtosis value for this distribution is always positive.

We also discussed a new class of nonlinearity in the level (or mode) model for capturing time series with heteroskedasticity and with multimodal conditional distribution, using Fisher’s z distribution as an innovation in the MAR model. The model offers great flexibility that other models, such as the TMAR and GMAR models, do not. The MCMC algorithm, using NUTS, allows for the easy estimation of the parameters in the model. The paper provides a simulation study using eight scenarios to indicate the flexibility and superiority of the ZMAR model compared with the TMAR and GMAR models. The simulation result shows that the ZMAR model is the best for representing the datasets generated from asymmetric components. When all the components are symmetrical, the ZMAR model also performs the best, as long as the excess unconditional kurtosis is large enough or the intercept distances between the components are far apart. However, when the datasets are generated from symmetrical components with small excess unconditional kurtosis and close intercept distances between the components, the GMAR model is the best. Furthermore, we compared the proposed model with the GMAR and TMAR models using two real data, namely the daily IBM stock prices and the monthly Brent crude oil prices. The results show that the proposed model outperforms the existing ones.

Fong et al. (2007) extended univariate the GMAR models to a Gaussian Mixture Vector Autoregressive (GMVAR) model. The ZMAR model can also be extended to a multivariate time-series context. Jones (2002) extended the standard multivariate F distribution to the multivariate skew t distribution and the multivariate Beta distribution. Likewise, the F distributed can also be extended to multivariate Fisher’s z distribution.