scieee AI-readable full text Open interactive document viewer

Fisher's z distribution-based mixture autoregressive model

Solikhah, Arifatus,Kuswanto, Heri,Iriawan, Nur,Fithriasari, Kartika

Abstract

EconStor is a publication server for scholarly economic literature, provided as a non-commercial public service by the ZBW.

Full text

Solikhah, Arifatus; Kuswanto, Heri; Iriawan, Nur; Fithriasari, Kartika Article Fisher's z distribution-based mixture autoregressive model Econometrics Provided in Cooperation with: MDPI – Multidisciplinary Digital Publishing Institute, Basel Suggested Citation: Solikhah, Arifatus; Kuswanto, Heri; Iriawan, Nur; Fithriasari, Kartika (2021) : Fisher's z distribution-based mixture autoregressive model, Econometrics, ISSN 2225-1146, MDPI, Basel, Vol. 9, Iss. 3, pp. 1-35, https://doi.org/10.3390/econometrics9030027 This Version is available at: https://hdl.handle.net/10419/247617 Standard-Nutzungsbedingungen: Die Dokumente auf EconStor dürfen zu eigenen wissenschaftlichen Zwecken und zum Privatgebrauch gespeichert und kopiert werden. Sie dürfen die Dokumente nicht für öffentliche oder kommerzielle Zwecke vervielfältigen, öffentlich ausstellen, öffentlich zugänglich machen, vertreiben oder anderweitig nutzen. Sofern die Verfasser die Dokumente unter Open-Content-Lizenzen (insbesondere CC-Lizenzen) zur Verfügung gestellt haben sollten, gelten abweichend von diesen Nutzungsbedingungen die in der dort genannten Lizenz gewährten Nutzungsrechte. Terms of use: Documents in EconStor may be saved and copied for your personal and scholarly purposes. You are not to copy documents for public or commercial purposes, to exhibit the documents publicly, to make them publicly available on the internet, or to distribute or otherwise use the documents in public. If the documents have been made available under an Open Content Licence (especially Creative Commons Licences), you may exercise further usage rights as specified in the indicated licence. https://creativecommons.org/licenses/by/4.0/ econometrics Article Fisher’s zDistribution-Based Mixture Autoregressive Model Arifatus Solikhah 1,2 , Heri Kuswanto 1,*, Nur Iriawan 1and Kartika Fithriasari 1   Citation: Solikhah, Arifatus, Heri Kuswanto, Nur Iriawan, and Kartika Fithriasari. 2021. Fisher’s z Distribution-Based Mixture Autoregressive Model. Econometrics 9: 27. https://doi.org/10.3390/ econometrics9030027 Academic Editor: In Choi Received: 10 December 2020 Accepted: 23 June 2021 Published: 29 June 2021 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). 1Department of Statistics, Faculty of Science and Data Analytics, Institut Teknologi Sepuluh Nopember, Surabaya 60111, Indonesia; [email protected] (A.S.); [email protected] (N.I.); [email protected] (K.F.) 2Badan Pusat Statistik (BPS—Statistics Indonesia), Jakarta 10710, Indonesia *Correspondence: [email protected] 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 zautoregressive models with the mixing proportions changing over time. This model can capture time series with both heteroskedasticity and multimodal conditional distribution, using Fisher’s zdistribution 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 zdistribution 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 tMixture 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. Keywords: Fisher’s zdistribution; mixture autoregressive model; the IBM stock prices; the Brent crude oil prices; Bayesian analysis; no-U-turn sampler; Stan program 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 KGaussian 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 Econometrics 2021,9, 27. https://doi.org/10.3390/econometrics9030027 https://www.mdpi.com/journal/econometrics Econometrics 2021,9, 27 2 of 35 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 zMixture Autoregressive (ZMAR) model which assumes that the distribution of innovations belongs to the Fisher’s zdistributions (Solikhah et al. 2021). The ZMAR model consists of a mixture of K-component Fisher’s zautoregressive models, where the numbers of components are based on the number of modes in the marginal density. The Fisher’s zdistribution’s mode is stable in its location parameters, whether it is symmetrical or skewed. Therefore, Fisher’s zuses 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 tdistribution (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 Lper 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 2describes the definition and properties of Fisher’s zdistribution in detail. In Section 3, we introduce the ZMAR model. Section 4demonstrates the flexibility of the ZMAR model compared with the TMAR and Econometrics 2021,9, 27 3 of 35 GMAR models using simulated datasets. Section 5contains 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 zDistribution Let Y be a random variable distributed as an Fdistribution with d1 and d2 degrees of freedom. The density of Z=1 2ln Ycan be defined as ζ(d1k,d2k)(z)=fZ(z;d1,d2)=2d2 d11 2d2 B1 2d1,1 2d2 e−d2z 1+e−2z+ln d2 d1(d1+d2)/2 , (1) and the cumulative distribution function (CDF) of Zis expressed as Z(d1k,d2k)(z)=Iz∗1 2d2,1 2d1=Rz∗ 0t1 2d2−1(1−t)1 2d1−1dt B1 2d1,1 2d2, (2) where e is the exponential constant; z∗=d2e−2z d1+d2e−2z , Iz∗(.) is the incomplete beta function ratio; and B(.) is the beta function, −∞<z<∞,d1>0, d2> 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 zdistribution. Let µ be a location parameter, and let σ be a scale parameter. The density of X=σZ+µis (Solikhah et al. 2021) fX(x;d1,d2,µ,σ)=2 σd2 d11 2d2 B1 2d1,1 2d2 e−d2(x−µ σ) 1+e−2(x−µ σ)+ln d2 d1(d1+d2)/2 , (3) where −∞<x<∞ , −∞<µ<∞ , and σ> 0. Equation (3) is defined as a p.d.f of Fisher’s zdistribution. It is denoted as z(d1,d2,µ,σ) . The CDF of the Fisher’s zdistribution is expressed as FX(x;d1,d2,µ,σ)=Ix∗1 2d2,1 2d1=1 B1 2d1,1 2d2 x∗ Z 0 t1 2d2−1(1−t)1 2d1−1dt, (4) where x∗=d2e−2(x−µ σ) d1+d2e−2(x−µ σ) . The quantile function (QF) of the Fisher’s zdistribution is defined as xp=µ+σ 2ln d2I−1xp1 2d1,1 2d2 d11−I−1xp1 2d1,1 2d2, (5) where d2I−1xp1 2d1,1 2d2/d11−I−1xp1 2d1,1 2d2 is the QF of the F-distribution and I−1xp(.) is the inversion of the incomplete beta function ratio. Let P−1xp(.) be the inversion of the incomplete gamma function ratio. The QF of the Fisher’s zdistribution can be expressed as xp=µ+σ 2ln  d2P−1v1p(d1/2) d1P−1v2p(d2/2) , (6) Econometrics 2021,9, 27 4 of 35 where 2 P−1v1p1 2d1 and 2 P−1v2p1 2d2 are the QF of the chi-square distribution with d1 and d2degrees of freedom, respectively. The proofs of Equation (1) up to Equation (6) are postponed to Appendix A. The parameters d1 and d2 , known as the shape parameters, are defined for both skewness (symmetrical if d1=d2 , asymmetrical if d16=d2 ) and fatness of the tails (large d1 and d2 imply thin tails). The Fisher’s zdistribution is also always unimodal and has the mode at x=µ . Furthermore, a change in the value of the parameter µ 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 zdistribution 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 zdistribution with various parameter settings can be seen in Figures 1and 2, respectively. Econometrics 2021, 9, x FOR PEER REVIEW 4 of 36 =+  2 ln        (   2 ⁄ )           (   2 ⁄ )   , (6) where 2 󰇡󰇢 and 2 󰇡󰇢 are the QF of the chi-square distribution with  and  degrees of freedom, respectively. The proofs of Equation (1) up to Equation (6) are postponed to Appendix A. The parameters  and , known as the shape parameters, are defined for both skewness (symmetrical if =, asymmetrical if ≠) and fatness of the tails (large  and  imply thin tails). The Fisher’s z distribution is also always unimodal and has the mode at =. Furthermore, a change in the value of the parameter  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 Figures 1 and 2, respectively. Figure 1. The p.d.f of the Fisher’s z distribution when =0 at various choices of , , and . Figure 2. The CDF of the Fisher’s z distribution when =0 at various choices of ,  and . Figure 1. The p.d.f of the Fisher’s zdistribution when µ=0 at various choices of d1,d2, and σ. Econometrics 2021, 9, x FOR PEER REVIEW 4 of 36 =+  2 ln        (   2 ⁄ )           (   2 ⁄ )   , (6) where 2 󰇡󰇢 and 2 󰇡󰇢 are the QF of the chi-square distribution with  and  degrees of freedom, respectively. The proofs of Equation (1) up to Equation (6) are postponed to Appendix A. The parameters  and , known as the shape parameters, are defined for both skewness (symmetrical if =, asymmetrical if ≠) and fatness of the tails (large  and  imply thin tails). The Fisher’s z distribution is also always unimodal and has the mode at =. Furthermore, a change in the value of the parameter  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 Figures 1 and 2, respectively. Figure 1. The p.d.f of the Fisher’s z distribution when =0 at various choices of , , and . Figure 2. The CDF of the Fisher’s z distribution when =0 at various choices of ,  and . Figure 2. The CDF of the Fisher’s zdistribution when µ=0 at various choices of d1,d2and σ. 3. Fisher’s zMixture Autoregressive Model 3.1. Model Specification Let yt ; t= 1, 2, . . . , T be the real-valued time series of interest; let F(t−1) denote the information set up to time t− 1, let FytF(t−1) ; k= 1, 2, · · · , K be the conditional CDF of Yt given the past information, evaluated at yt ; and K is the number of components in the Econometrics 2021,9, 27 5 of 35 ZMAR model. Let Z(d1k,d2k)(.)be the CDF of the standardized Fisher’s zdistribution with d1k and d2k shape parameters, given by Equation (2); let ek.t be a sequence of independent standardized Fisher’s zrandom variables such that ek.t is independent of {yt−i,i>0} ; and let σk is a scale parameter of the kth component. The K -component ZMAR can be defined as FytF(t−1)= K ∑ k=1 ηkZ(d1k,d2k)yt−µk.t σk, (7) or yt=             µ1.t+σ1e1.t; with probability η1; µ2.t+σ2e2.t; with probability η2; . . . µK.t+σKeK.t; with probability ηK, (8) with µk.t=φk.0 +∑pk i=1φk.iyt−i;k=1, 2, · · · ,K, (9) where the vector η=(η1, ...,ηK) is called the weights. η takes a value in the unit simplex EK , which is a subspace of (<+)K , defined by the following constraint, ηk> 0 and K ∑ k=1 ηk= 1. We use the abbreviation ZMAR(K;p1,p2, . . . , pK) for this model, with the parameter ϑ=(θ1,θ2, . . . , θK,η) ; η=(ηk) ; θk=(d1k,d2k,σk,φk.0,φk) ; φk=φk.1,φk.2, . . . , φk.pk ; to each value k∈{1, 2, . . . , K} ; taking values in the parameter space ΘK=ΘK× EK , φk.i denotes the AR coefficient on the kth component and ith lag; i= 1, 2, . . . , pk ; and pk denotes the autoregressive order of the kth component. Using the parameters θk , we first define the Kauxiliary Fisher’s zAR(pk) processes fkytF(t−1)=φk.0 +∑pk i=1φk.iyt−i+σkek.t;k=1, 2, . . . , K, where the AR coefficients φkare assumed to satisfy 1−∑pk i=1φk.iCi6=0 for |C|≤1; k=1, . . . , K. This condition implies that the processes fkytF(t−1) are stationary and that each component model in (7) or (8) satisfies the usual stationarity condition of the linear AR( pk ) model. Suppose p=max(p1,p2, ..., pK) . Let a univariate time series y=(y1,y2,· · · ,yT) be influenced by a hidden discrete indicator variables Q=(Q1,Q2,· · · ,QT) , where Qt takes values in the set {1, 2, . . . , K} . The probability of sampling from the group labeled Qt=k is equal to ηk . Suppose that the conditional density of Yt given Qt=k is fkytF(t−1) ; (k=1, 2, . . . , K) . Let H=(H10,H20, ..., HT0)0 be the unobserved random variable, where Ht is a K -dimensional vector with the values h=(h10,h20, ..., hT0)0 , where hk.t= 1 if Qt=k and hk.t= 0, otherwise. Thus, Ht is distributed according to a multinomial distribution consisting of one draw on K categories with probabilities η1,η2, . . . , ηK(McLachlan and Peel 2000, p. 7); that is, P(Ht=ht)=η1h1.tη2h2.t. . . ηKhK.t. (10) Otherwise, this can be written as Ht∼MultK(1, η) ; where η=(η1,η2, ...,ηK) . The conditional likelihood for the ZMAR(K;p1,p2, . . . , pK)can be formed by P(yt|ϑ)= K ∏ k=1 T ∏ t=p+1 [ηk∆1∆2]hk.t, (11) Econometrics 2021,9, 27 6 of 35 where ∆1=2(d2k/d1k)1 2d2k σkB(1 2d1k,1 2d2k),∆2=exp(−d2k(ak.t/σk)) [1+exp(−2(ak.t/σk)+ln[d2k/d1k])](d1k+d2k)/2 , and ak.t=yt−µk.t. Let τt.k be the probability for each t-th observation, t=(p+1) , (p+2) , . . . , T , as members of the k-th component, k= 1, 2, . . . , K of a mixture distribution. Suppose yt=(y1,y2,· · · ,yt) , Bayes’ rule to compute the τt.k can be expressed as (Frühwirth- Schnatter 2006, p. 26) PQt=kyt,ϑ=τt.k=ηkfkytF(t−1) ∑K k=1ηkfkytF(t−1). (12) Let ζ(d1k,d2k)(.) be the p.d.f of the standardized Fisher’s zdistribution and given by Equation (1). Then, the mixing weights for the ZMAR model can be expressed as τt.k=ηkσk−1ζ(d1k,d2k)(ak.t/σk) ∑K k=1ηkσk−1ζ(d1k,d2k)(ak.t/σk);k=1, 2, . . . , K. Suppose ωk.t and ξk2 signify the conditional mean and conditional variance of the kth component, which are defined by ωk.t=E[Yt|θk]=µk.t+σk 2ln d2k d1k −ψ1 2d2k+ψ1 2d1k, (13) and ξk2=Var[Yt|θk]=σk 22ψ01 2d2k+ψ01 2d1k, (14) where ψ(.) is the digamma function and ψ0(.) is the trigamma function. The conditional mean of Yt(Frühwirth-Schnatter 2006, p. 10) is obtained as ωt=EhYt|F(t−1)i= K ∑ k=1 ηkωk.t, (15) the conditional variance of Yt(Frühwirth-Schnatter 2006, p. 11) is obtained as ξ2=VarhYt|F(t−1)i= K ∑ k=1 ηkωk.t2+ξk2−ωt2, (16) and the higher-order moments around the mean of Yt (Frühwirth-Schnatter 2006, p. 11) are obtained as E(Y−ωt)3F(t−1)= K ∑ k=1(ωk.t−ωt)2+3ξk2(ωk.t−ωt)ηk, (17) E(Y−ωt)4F(t−1)= K ∑ k=1(ωk.t−ωt)4+6(ωk.t−ωt)2ξk2+3ξk4ηk. (18) These expressions apply to any specification of the mixing weights ηk. 3.2. Bayesian Approach for ZMAR Model In this paper, we apply a Bayesian method to estimate the parameters ϑ . The Bayesian analysis requires the joint posterior density π(ϑ|y), which is defined by π(ϑ|y)∝P(yt|ϑ)π(ϑ), (19) Econometrics 2021,9, 27 7 of 35 where P(yt|ϑ) is the conditional likelihood function given by Equation (11) and π(ϑ) is the prior of the parameter model, which is defined by π(ϑ)=π(η) K ∏ k=1 π(d1k)π(d2k)π(σk)π(φk.0) pk ∏ i=1 π(φk.i), (20) where the π(η) is prior for the η parameter, the π(d1k) , π(d2k) , π(σk) , π(φk.0) and π(φk.i) are prior for the d1 , d2 , σ , φ0 and φi parameters at the kth component; index i denotes the ith lag; i=1, 2, . . . , pk;k=1, 2, . . . , 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-tfamily 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 tdistribution (positive values only) (Kim 2008) for the priors of the d1k , d2k , and σk , with the degrees of freedoms v1k , v2k , ν3k , the location parameters m1k , m2k , m3k , and the scale parameters s1k2 , s2k2 , s3k2 , respectively. Therefore, it can be written as d1k∼tν1km1k,s1k2I(0, ∞) , d2k∼tν2km2k,s2k2I(0, ∞) , and σk∼tν3km3k,s3k2I(0, ∞) . We take the Dirichlet distribution (Kotz et al. 2000, p. 485) for the prior of the η parameter, thus η1,η2, . . . , ηK−1∼Dir(δ1,δ2, . . . , δK) . For the priors of the φk.0 and φk.i , we take the normal distribution with the location parameters uk.0 , uk.i and the scale parameters gk.02 , gk.i2 , thus φk.0 ∼Nuk.0,gk.02 and φk.i∼Nuk.i,gk.i2 ; i= 1, 2, . . . , pk . Employing the setup of prior distributions, as shown above, the natural logarithm of the joint posterior distribution of the model is given by ln π(ϑ|y)∝ K ∑ k=1 T ∑ t=p+1 hk.t(ln ηk+Λ1+Λ2)!+Λ3+Λ4+Λ5+Λ6+Λ7+ pk ∑ i=1 Λ8!,(21) where Λ1=d2kln[d2k/d1k]+lnhΓ1 2(d1k+d2k)i−lnhΓ1 2d1ki−lnhΓ1 2d2ki−ln σk , Λ2=−d2kak.t σk−d1k+d2k 2lnh1+exph−2ak.t σk+ln[d2k/d1k]ii , Λ3=(δK−1)ln1−K−1 ∑ k=1 ηk +K−1 ∑ k=1 (δk−1)ln ηk , Λ4=−1 2(ν1k+1)ln1+1 ν1kd1k−m1k s1k2 , Λ5=−1 2(ν2k+1) ln1+1 ν2kd2k−m2k s2k2 , Λ6=−1 2(ν3k+1)ln1+1 ν3kσk−m3k s3k2 , Λ7=−1 2φk.0−uk.0 gk.0 2 , and Λ8=−1 2φk.i−uk.i gk.i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 `(ϑ|y)is ∂ ∂ϑ`(ϑ|y)=∂ ∂ϑln π(ϑ|y) ∂ ∂ϑ`(ϑ|y)=∂`(ϑ|y) ∂ηk ,∂`(ϑ|y) ∂d1k ,∂`(ϑ|y) ∂d2k ,∂`(ϑ|y) ∂σk ,∂`(ϑ|y) ∂φk.0 ,∂`(ϑ|y) ∂φk.i(22) where i= 1, 2, . . . , pk ; k= 1, 2, . . . , K . The HMC algorithm for estimating the parameters of the ZMAR model is as follows: 1. Determine the initial value of the parameter ϑ0 , the diagonal mass matrix M , the scale factor of the leapfrog steps ε , the number of leapfrog steps N , and the number of iterations R. 2. For each iteration r;r=1, 2, . . . , Rwhere Rrepresents the number of iterations, Econometrics 2021,9, 27 8 of 35 a. Generate the momentum variables ϕwith ϕ∼Normal(0,M); b. For each iteration n=1, 2, . . . , N, (1) Use the gradient of the ln-posterior density of ϑn to make a half-step of ϕn ϕn+0.5 ←ϕn+1 2ε∂ln π(ϑn|y) ∂ϑn, (2) Update the vector ϑnusing the vector ϕn+0.5 ϑn+0.5 ←ϑn+εM−1ϕn+0.5 (3) Update the next half-step for ϕ ϕn+1←ϕn+0.5 +1 2ε∂ln π(ϑn+0.5|y) ∂ϑn+0.5 . c. Labels ϑ(r−1) and ϕ(r−1) as the value of the parameter and momentum vectors at the start of the leapfrog process and ϑ#,ϕ#as the value after the Nsteps. d. Compute b=π(ϑ#|y)f(ϕ#) π(ϑ(r−1)|y)f(ϕ(r−1)). e. Set ϑ(r) ϑ(r)=ϑ#with probability min(b, 1) ϑ(r−1)otherwise. f. Save ϑ(r);r=1, 2, . . . , R. The performance of the HMC is very sensitive to two user-defined parameters, i.e., the step size of the leapfrog ε 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 ε 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 zdistributions. 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 ZMAR(3; 1, 1, 1) where fytF(t−1)=3 ∑ k=1 ηkfkytF(t−1)=3 ∑ k=1 ηk(φk.0 +φk.1yt−1+ek.t) ; η1=η2=η3= 1 3 ; φ1.1 =− 0.6; φ2.1 = 0.2; φ3.1 = 0.7 with ek.t∼z(d1.k,d2.k, 0, σk) ; 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 φ1.0 =φ2.0 =φ3.0 = 0, excess unconditional kurtosis is 5.73, e1.t∼z(0.2, 10, 0, 5) , e2.t∼z(1, 1, 0, 8) , and e3.t∼z(30, 30, 0, 10); # Scenario 2: represent a mixture of a highly skewed to the right and two symmetrical distributions, where φ1.0 =φ2.0 =φ3.0 = 0, excess unconditional kurtosis is 2.23, e1.t∼z(20, 1, 0, 5),e2.t∼z(1, 1, 0, 8), and e3.t∼z(30, 30, 0, 10); # Scenario 3: represent a mixture of three symmetrical distributions, where φ1.0 = φ2.0 =φ3.0 = 0, excess unconditional kurtosis is 1.67, e1.t∼z(0.5, 0.5, 0, 5) , e2.t∼z(1, 1, 0, 8), and e3.t∼z(30, 30, 0, 10); Econometrics 2021,9, 27 15 of 35 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 3shows 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 Appendices F.2 and G.2. 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: eta[1]=η1 ; eta[2]=η2 ; eta[3]=η3 ; sigma[1]=σ1 ; sigma[2]=σ2 ; sigma[3]=σ3 ; d 1 [1]=d11 ; d 1 [2]=d12 ; d 1 [3]=d13 ; d 2 [1]= d21 ; d 2 [2]=d22 ; d 2 [3]=d23 ; nu[1]=ν1 ; nu[2]=ν2 ; nu[3]=ν3 ; phi 1 [1]=φ1.1 ; phi 1 [2]=φ1.2 ; phi 1 [3]=φ1.3 ; phi 2 [1]=φ2.1 ; phi 2 [2]= φ2.2; phi2[3]=φ2.3; phi3[1]=φ3.1; phi3[2]=φ3.2; phi3[3]=φ3.3. For all the parameter models, the MCMC chains reached convergence, which was shown by the ˆ Rstatistic being less than 1.01 and the ne f f statistic being greater than 400. The three-component ZMAR model for the original series, namely FytF(t−1)=0.40 Z(13.22,4.49)yt−0.63 yt−1−0.37 yt−2 5.14  +0.39 Z(0.99,1.91)yt−1.68 yt−1+1.02 yt−2−0.34yt−3 2.82  +0.21 Z(10.09,4.37)yt−1.69 yt−1−0.03 yt−2+0.72yt−3 7.12 . (26) a three-component TMAR model for the original series, namely FytF(t−1)=0.34 T14.98yt−1.61 yt−1+1.03yt−2−0.42yt−3 4.93  +0.40 T12.08yt−0.72 yt−1−0.28 yt−2 1.62  +0.26 T4.34yt−1.54 yt−1−0.32 yt−2+1.14 yt−3−0.28 yt−4 1.80 . (27) Econometrics 2021,9, 27 16 of 35 and a three-component GMAR for the original series, namely FytF(t−1)=0.45 Φyt−1.14 yt−1+0.49 yt−2−0.16 yt−3−0.19 yt−4 4.18  +0.28 Φyt−1.48 yt−1+0.65 yt−2−0.55 yt−3+0.38 yt−4 1.66  +0.27 Φyt−1.57 yt−1−0.34 yt−2+1.20 yt−3−0.29 yt−4 1.92 . (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 zdistribution are µ , σ , d1 , and d2 . The µ is a location parameter, the σ is a scale parameter, and the d1 , d2 are known as the shape parameters, defined for both skewness (symmetric if d1=d2 , asymmetric if d16=d2 ) and fatness of the tails (large d1 and d2 imply thin tails). The Fisher’s zdistribution is always unimodal and has the mode at x=µ . The value of µ only affects the mean of the distribution. It does not affect the variance, skewness, and kurtosis of the distribution. Furthermore, if d1=d2 , then the mean is equal to µ ; if d1<d2 , then the mean is less than µ ; and if d1>d2 , then the mean is greater than µ . 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 zdistribution 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 Fdistribution to the multivariate skew tdistribution and the multivariate Beta distribution. Likewise, the F distributed can also be extended to multivariate Fisher’s zdistribution. Author Contributions: A.S., H.K., N.I. and K.F. analyzed and designed the research; A.S. collected, analyzed the data, and drafted the paper. All authors critically read and revised the draft and approved the final paper. All authors have read and agreed to the published version of the manuscript. Funding: This research received no external funding. Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable. Data Availability Statement: Data available from stated sources. Conflicts of Interest: The authors declare no conflict of interest. Econometrics 2021,9, 27 17 of 35 Appendix A. Proofs of the Equations (1)–(6) Appendix A.1. Proof of the Equation (1) Solikhah et al. (2021) described transforming a random variable with the Fdistribution to the Fisher’s zdistribution. Let Y be a random variable distributed as an Fdistribution with two parameters d1and d2. The density of Z=1 2ln Yis (Fisher 1924;Aroian 1941) fZ(z;d1,d2)=2d1 1 2d1d2 1 2d2 B1 2d1,1 2d2 ed1z (d1e2z+d2)(d1+d2)/2 , (A1) where −∞<z<∞,d1>0, d2>0 and B(.)is the beta function. Interchanging d1and d2 is equivalent to replacing z with −z (Fisher 1924;Aroian 1941), thus Equation (A1) can also be defined as fZ(z;d1,d2)=2d1 1 2d1d2 1 2d2 B1 2d1,1 2d2 e−d2z (d2e−2z+d1)(d1+d2)/2 . (A2) If the denominator and numerator of Equation (A2) are divided by d1(d1+d2)/2 then we get fZ(z;d1,d2)=2d2 d11 2d2 B1 2d1,1 2d2 e−d2z 1+d2 d1e−2z(d1+d2)/2 . (A3) Equation (A3) can also be defined as fZ(z;d1,d2)=2d2 d11 2d2 B1 2d1,1 2d2 e−d2z 1+e−2z+ln d2 d1(d1+d2)/2  Appendix A.2. Proof of the Equation (2) Let Ybe a random variable distributed as an Fdistribution with d1 and d2 degrees of freedom, let Iy∗(.) be the incomplete beta function ratio, and let B(.) be the beta function. The cumulative distribution function (CDF) of the Y=e2z be defined as follows (Johnson et al. 1995, vol. 2, p. 327) FY(y;d1,d2)=Iy∗1 2d1,1 2d2=Ry∗ 0t1 2d1−1(1−t)1 2d2−1dt B1 2d1,1 2d2, where y∗=d1y d2+d1y. If Z=1 2ln Ythen Y=e2z, thus the CDF of Zis FZ(z;d1,d2)=Iz∗1 2d1,1 2d2=Rz∗ 0t1 2d1−1(1−t)1 2d2−1dt B1 2d1,1 2d2, (A4) where z∗=d1e2z d2+d1e2z. Equation (A4) can also be defined as FZ(z;d1,d2)=Iz∗1 2d2,1 2d1=Rz∗ 0t1 2d2−1(1−t)1 2d1−1dt B1 2d1,1 2d2, where z∗=d2e−2z d1+d2e−2z. Econometrics 2021,9, 27 18 of 35 Appendix A.3. Proof of the Equation (3) Let Z be a random variable distributed as a standardized Fisher’s zdistribution with the p.d.f given by Equation (1), let µ be a location parameter, and let σ be a scale parameter. The density of X=σZ+µcan be defined as fX(x;d1,d2,µ,σ)=fZx−µ σ|J(x)|, (A5) where J(x)is the Jacobian of the transformation and is defined as J(x)=∂z ∂x=1 σ. Therefore, the p.d.f of the Fisher’s zdistribution can be expressed as fX(x;d1,d2,µ,σ)=2 σd2 d11 2d2 B1 2d1,1 2d2 e−d2(x−µ σ) 1+e−2(x−µ σ)+ln d2 d1(d1+d2)/2  Appendix A.4. Proof of the Equation (4) Let Z be a random variable distributed as a standardized Fisher’s zdistribution with the CDF as in Equation (2), let µ be a location parameter, and let σ be a scale parameter. The CDF of X=σZ+µcan be defined as FX(x;d1,d2,µ,σ)=P(X≤x)=P(σZ+µ≤x)=PZ≤x−µ σ, Therefore, the CDF of the Fisher’s zdistribution can be expressed as FX(x;d1,d2,µ,σ)=1 B1 2d1,1 2d2 x∗ Z 0 t1 2d2−1(1−t)1 2d1−1dt, where x∗=d2e−2(x−µ σ) d1+d2e−2(x−µ σ). Appendix A.5. Proof of the Equation (5) Let Ix∗(.) be the incomplete beta function ratio; thus, Equation (4) can be expressed as FX(x;d1,d2,µ,σ)=Ix∗(d2/2, d1/2);x∗=d2e−2(x−µ σ) d1+d2e−2(x−µ σ) The value xp is called the p-quantile of the population, if PX≤xp=p with 0 ≤ p≤ 1 (Gilchrist 2000, p. 12). Let I−1xp(.) be the inversion of the incomplete beta function ratio, then I−1xp(d2/2, d1/2)=d2e−2(xp−µ σ) d1+d2e−2(xp−µ σ); d1+d2e−2(xp−µ σ)I−1xp(d2/2, d1/2)=d2e−2(xp−µ σ); (d1)I−1xp(d2/2, d1/2)=d2e−2(xp−µ σ)−d2e−2(xp−µ σ)I−1xp(d2/2, d1/2); (d1)I−1xp(d2/2, d1/2)=1−I−1xp(d2/2, d1/2)d2e−2(xp−µ σ); d1(I−1xp(d2/2,d1/2)) d2(1−(I−1xp(d2/2,d1/2))) =e−2(xp−µ σ); Econometrics 2021,9, 27 19 of 35 xp=µ−σ 2ln  d1I−1xp(d2/2, d1/2) d21−I−1xp(d2/2, d1/2) ; Interchanging d1 and d2 is equivalent to replacing x with −x ; thus, the QF can also be defined as: xp=µ+σ 2ln  d2I−1xp(d1/2, d2/2) d11−I−1xp(d1/2, d2/2)   Appendix A.6. Proof of the Equation (6) Let us denote a chi-square random variables with d1 and d2 degrees of freedom by V2 1 and V2 2 , respectively. Let Pv∗(.) be the incomplete gamma function ratio, and the CDF of V2 1and V2 2can be defined as FV2 1(v1;d1)=Pv1∗(d1/2);v1∗=v1 2, FV2 2(v2;d2)=Pv2∗(d2/2);v2∗=v2 2. Let P−1vp(.) be the inversion of the incomplete gamma function ratio, and then the QF of the V2 1and V2 2can be defined as v1p=2P−1v1p(d1/2), v2p=2P−1v2p(d2/2), The Beta distribution arises naturally as the distribution of X=V2 1/V2 1+V2 2 (Johnson et al. 1995, vol. 2, p. 212) ; therefore, QF of the Fisher’s zdistribution can be expressed as xp=µ+σ 2ln    d2 2P−1v1p(d1/2) 2P−1v1p(d1/2)+2P−1v2p(d2/2) d11−2P−1v1p(d1/2) 2P−1v1p(d1/2)+2P−1v2p(d2/2)    xp=µ+σ 2ln  d2P−1v1p(d1/2) d1P−1v2p(d2/2)   Appendix B. Properties of the Fisher’s z Distribution and the Proofs Appendix B.1. Properties of the Fisher’s z Distribution Let Xbe a random variable distributed as a Fisher’s zdistribution and let MX(θ) be the Moment Generating Function (MGF) of a random variable X. The MGF of the Fisher’s zdistribution be expressed as MX(θ)=Eeθx=eθµ d2 d1σθ 2!Γ1 2(d2−σθ)Γ1 2(d1+σθ) Γ1 2d1Γ1 2d2, (A6) where Γ(.) is a gamma function. Let KX(θ) be the Cumulant Generating Function (CGF) of a random variable X. The CGF of the Fisher’s zdistribution be given by Econometrics 2021,9, 27 20 of 35 KX(θ)=θµ +θσ 2ln d2 d1+ln Γ1 2(d1+σθ)+ln Γ1 2(d2−σθ) −ln Γ1 2d1−ln Γ1 2d2.(A7) The coefficient of θ ! /j ! in the Taylor expansion of the CGF is the j th cumulant of X and be denoted as κj . The j th cumulant, therefore, can be obtained by differentiating the expansion jtimes and evaluating the result at zero. κj=∂j ∂θjKX(θ)θ=0 ;j=1, 2, 3, · · · The first cumulant of the Fisher’s zdistribution be defined as κ1=µ+σ 2ln d2 d1 −σ 2ψ1 2d2+σ 2ψ1 2d1, (A8) and the jth cumulant would be defined in Equation (A9). κj=(−1)jσ 2jψ(j−1)1 2d2+σ 2jψ(j−1)1 2d1;j=2, 3, 4, · · · , (A9) where ψ(.) is the digamma function, ψ0(.) is the trigamma function, ψ00 (.) is the tetragamma function, and ψ(3)(.) is the pentagamma function. Generally, ψ(j−1)(.) is the (j+1) -gamma function (Johnson et al. 2005, p. 9). Let the mean and the variance of a random variable Xbe denoted respectively E(X) and Var(X) . The mean of the Fisher’s zdistribution is given by E(X)=κ1=µ+σ 2ln d2 d1 −σ 2ψ1 2d2+σ 2ψ1 2d1, (A10) and the variance is defined as Var(X)=κ2=σ 22ψ01 2d2+ψ01 2d1. (A11) On the basis of Equation (A10), it can be concluded that (d1=d2)−→ (E(X)=µ); (d1<d2)−→ (E(X)<µ); (d1>d2)−→ (E(X)>µ). Let the skewness and the excess kurtosis of a random variable Xbe denoted, respectively, as γ1and γ2. The skewness of the Fisher’s zdistribution is given by γ1=ψ00 1 2d1−ψ00 1 2d2 ψ01 2d1+ψ01 2d23 2 , (A12) and the excess kurtosis is γ2=ψ(3)d2 2+ψ(3)d1 2 ψ0d2 2+ψ0d1 22. (A13) On the basis of Equation (A13), the excess kurtosis value for this distribution is always positive, which shows that the distribution has heavier tails than the Gaussian distribution. Furthermore, based on Equation (A10) through Equation (A13), it can be seen that a change Econometrics 2021,9, 27 21 of 35 in the value of the parameter µ only affects the mean of the distribution. It does not affect the variance, skewness, and kurtosis of the distribution. Appendix B.2. Proof of the Properties Appendix B.2.1. Proof of the Equation (A6) If Z is random variables distributed as a standardized Fisher’s z, then the MGF of Z is expressed as (Aroian 1941;Johnson et al. 1995) MZ(θ)= d2 d1θ 2!Γ1 2(d1+θ)Γ1 2(d2−θ) Γ1 2d1Γ1 2d2 If the random variable Zis transformed to X=σZ+µ, then MX(θ)=EeθX=Eeθ(µ+σZ)=eθµEeθ(σZ)=eθµ MZ(θσ) MX(θ)=eθµ d2 d1θσ 2!Γ1 2(d2−θσ)Γ1 2(d1+θσ) Γ1 2d1Γ1 2d2  Appendix B.2.2. Proof of the Equation (A7) The CGF of the random variable X is the natural logarithm of the moment generating function of X(Johnson et al. 2005, p. 54), therefore KX(θ)=ln MX(θ)=ln eθµ d2 d1θσ 2!Γ1 2(d1+θσ)Γ1 2(d2−θσ) Γ1 2d1Γ1 2d2  KX(θ)=θµ +θσ 2lnd2 d1+ln Γ1 2(d1+θσ)+ln Γ1 2(d2−θσ) −ln Γ1 2d1−ln Γ1 2d2  Appendix B.2.3. Proof of the Equations (A8) and (A9) If the random variable Xhas the CGF in the Equation (A7), then ∂ ∂θ KX(θ)=µ+σ 2lnd2 d1−σ 2ψd2−θσ 2+σ 2ψd1+θσ 2 ∂2 ∂θ2KX(θ)=(−1)2σ 22ψ0d2−θσ 2+σ 22ψ0d1+θσ 2 ∂3 ∂θ3KX(θ)=(−1)3σ 23ψ00 d2−θσ 2+σ 23ψ00 d1+θσ 2 . . . ∂j ∂θjKX(θ)=(−1)jσ 2jψ(j−1)d2−θσ 2+σ 2jψ(j−1)d1+θσ 2. The first cumulant be defined as κ1=∂ ∂θ KX(θ)θ=0 =µ+σ 2lnd2 d1−σ 2ψd2−θσ 2+σ 2ψd1+θσ 2θ=0 =µ+σ 2lnd2 d1−σ 2ψd2 2+σ 2ψd1 2  Econometrics 2021,9, 27 22 of 35 and the jth cumulant be defined as κj=∂j ∂θjKX(θ)θ=0 =(−1)jσ 2jψ(j−1)d2−θσ 2+σ 2jψ(j−1)d1+θσ 2θ=0 =(−1)jσ 2jψ(j−1)1 2d2+σ 2jψ(j−1)1 2d1;j=2, 3, 4, . . .  Appendix B.2.4. Proof of the Equation (A10) The mean of the random variable X is the first cumulant (Zelen and Severo 1970), therefore E(X)=κ1 E(X)=µ+σ 2lnd2 d1−σ 2ψd2 2+σ 2ψd1 2  Appendix B.2.5. Proof of the Equation (A11) The variance of the random variable X is the second cumulant (Zelen and Severo 1970), therefore Var(X)=κ2 Var(X)=σ 22ψ0d2 2+σ 22ψ0d1 2=σ 22ψ0d2 2+ψ0d1 2  Appendix B.2.6. Proof of the Equation (A12) The skewness is formed from the second and third cumulants (Zelen and Severo 1970), namely γ1=κ3 (κ2)3/2 =σ 23−ψ00 d2 2+ψ00 d1 2 σ 23ψ0d2 2+ψ0d1 23 2 ; γ1=ψ00 1 2d1−ψ00 1 2d2 ψ01 2d1+ψ01 2d23 2  Appendix B.2.7. Proof of the Equation (A13) The excess kurtosis can be formed from the second and fourth cumulants (Zelen and Severo 1970), namely γ2=κ4 (κ2)2=σ 24ψ(3)d2 2+ψ(3)d1 2 σ 24ψ0d2 2+ψ0d1 22; γ2=ψ(3)d2 2+ψ(3)d1 2 ψ0d2 2+ψ0d1 22  Econometrics 2021,9, 27 23 of 35 Appendix C. Adding the Fisher’s z Distribution Functions in Stan Appendix C.1. Random Numbers Generator Function We can add the random numbers generator function of Fisher’s zdistribution (fisher_z_rng) in Stan using the following code, functions{ real fisher_z_rng(real d1, real d2, real mu, real sigma){ return(mu+sigma*0.5*log((chi_square_rng(d1)*d2)/(chi_square_rng(d2)*d1))); } } where chi_square_rng(d1) and chi_square_rng(d2) are the chi-square random numbers generator with d1and d2degrees of freedoms. Appendix C.2. Log Probability Density Function We can also add the log probability function of the Fisher’s zdistribution (fisher_z_lpdf) in Stan using the following code, 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))); } } where lbeta(0.5*d1,0.5*d2) is the natural logarithm of the beta function applied to 1 2d1 and 1 2d2 , and log1p_exp((-2*(x-mu)/sigma)+log(d2)-log(d1))) is the natural logarithm of one plus the natural exponentiation of −2x−µ σ+ln d2−ln d1. Appendix D. Priors of Parameters on the ZMAR, TMAR, and GMAR Models in the Simulation Study Appendix D.1. Scenario 1 Table A1. Priors of parameters on the ZMAR, TMAR, and GMAR models in the Scenario 1. ZMAR TMAR GMAR d11 ∼t3(0.2, 0.1)I(0, ∞); d21 ∼t3(10, 0.1)I(0, ∞); σ1∼t3(5, 0.1)I(0, ∞) φ1.1 ∼N(−0.6, 0.1) d12 ∼t3(1, 0.1)I(0, ∞); d22 ∼t3(1, 0.1)I(0, ∞); σ2∼t3(8, 0.1)I(0, ∞) φ2.1 ∼N(0.2, 0.1) d13 ∼t3(30, 0.1)I(0, ∞); d23 ∼t3(30, 0.1)I(0, ∞); σ3∼t3(10, 0.1)I(0, ∞) φ3.1 ∼N(0.7, 0.1) η∼Dir(1, 1, 1) ν1∼t3(2.1, 0.1)I(0, ∞); σ1∼t3(14.6, 0.1)I(0, ∞) φ1.1 ∼N(−0.62, 0.1) ν2∼t3(3.94, 0.1)I(0, ∞); σ2∼t3(8.1, 0.1)I(0, ∞) φ2.1 ∼N(0.2, 0.1) ν3∼t3(9.66, 0.1)I(0, ∞); σ3∼t3(1.62, 0.1)I(0, ∞) φ3.1 ∼N(0.84, 0.1) η∼Dir(1, 1, 1) σ1∼t3(23, 0.1)I(0, ∞) φ1.1 ∼N(−0.07, 0.1) σ2∼t3(33, 0.1)I(0, ∞) φ2.1 ∼N(−0.47, 0.1) σ3∼t3(11, 0.1)I(0, ∞) φ3.1 ∼N(0.6, 0.1) σ4∼t3(2, 0.1)I(0, ∞) φ4.1 ∼N(0.8, 0.1) η∼Dir(1, 1, 1, 1) Econometrics 2021,9, 27 24 of 35 Appendix D.2. Scenario 2 Table A2. Priors of parameters on the ZMAR, TMAR, and GMAR models in the Scenario 2. ZMAR TMAR GMAR d11 ∼t3(20, 0.1)I(0, ∞); d21 ∼t3(1, 0.1)I(0, ∞); σ1∼t3(5, 0.1)I(0, ∞) φ1.1 ∼N(−0.6, 0.1) d12 ∼t3(1, 0.1)I(0, ∞); d22 ∼t3(1, 0.1)I(0, ∞); σ2∼t3(8, 0.1)I(0, ∞) φ2.1 ∼N(0.2, 0.1) d13 ∼t3(30, 0.1)I(0, ∞); d23 ∼t3(30, 0.1)I(0, ∞); σ3∼t3(10, 0.1)I(0, ∞) φ3.1 ∼N(0.7, 0.1) η∼Dir(1, 1, 1) ν1∼t3(2.19, 0.1)I(0, ∞); σ1∼t3(8, 0.1)I(0, ∞) φ1.1 ∼N(−0.62, 0.1) ν2∼t3(3.94, 0.1)I(0, ∞); σ2∼t3(8.12, 0.1)I(0, ∞) φ2.1 ∼N(0.2, 0.1) ν3∼t3(9.66, 0.1)I(0, ∞); σ3∼t3(1.62, 0.1)I(0, ∞) φ3.1 ∼N(0.84, 0.1) η∼Dir(1, 1, 1) σ1∼t3(5.01, 0.1)I(0, ∞) φ1.1 ∼N(−0.59, 0.1) σ2∼t3(2.69, 0.1)I(0, ∞) φ2.1 ∼N(0.64, 0.1) σ3∼t3(14.49, 0.1)I(0, ∞) φ3.1 ∼N(0.57, 0.1) η∼Dir(1, 1, 1) Appendix D.3. Scenario 3 Table A3. Priors of parameters on the ZMAR, TMAR, and GMAR models in the Scenario 3. ZMAR TMAR GMAR d11 ∼t3(0.5, 0.1)I(0, ∞); d21 ∼t3(0.5, 0.1)I(0, ∞); σ1∼t3(5, 0.1)I(0, ∞) φ1.1 ∼N(−0.6, 0.1) d12 ∼t3(1, 0.1)I(0, ∞); d22 ∼t3(1, 0.1)I(0, ∞); σ2∼t3(8, 0.1)I(0, ∞) φ2.1 ∼N(0.2, 0.1) d13 ∼t3(30, 0.1)I(0, ∞); d23 ∼t3(30, 0.1)I(0, ∞); σ3∼t3(10, 0.1)I(0, ∞) φ3.1 ∼N(0.7, 0.1) η∼Dir(1, 1, 1) ν1∼t3(2.1, 0.1)I(0, ∞); σ1∼t3(14.6, 0.1)I(0, ∞) φ1.1 ∼N(−0.62, 0.1) ν2∼t3(3.94, 0.1)I(0, ∞); σ2∼t3(8.1, 0.1)I(0, ∞) φ2.1 ∼N(0.2, 0.1) ν3∼t3(9.66, 0.1)I(0, ∞); σ3∼t3(1.62, 0.1)I(0, ∞) φ3.1 ∼N(0.84, 0.1) η∼Dir(1, 1, 1) σ1∼t3(10, 0.1)I(0, ∞) φ1.1 ∼N(−0.6, 0.1) σ2∼t3(8, 0.1)I(0, ∞) φ2.1 ∼N(0.24, 0.1) σ3∼t3(2, 0.1)I(0, ∞) φ3.1 ∼N(0.74, 0.1) η∼Dir(1, 1, 1) Appendix D.4. Scenario 4 Table A4. Priors of parameters on the ZMAR, TMAR, and GMAR models in the Scenario 4. ZMAR TMAR GMAR d11 ∼t3(3, 0.1)I(0, ∞); d21 ∼t3(10, 0.1)I(0, ∞); σ1∼t3(8, 0.1)I(0, ∞) φ1.1 ∼N(−0.6, 0.1) d12 ∼t3(1, 0.1)I(0, ∞); d22 ∼t3(1, 0.1)I(0, ∞); σ2∼t3(5, 0.1)I(0, ∞) φ2.1 ∼N(0.2, 0.1) d13 ∼t3(30, 0.1)I(0, ∞); d23 ∼t3(30, 0.1)I(0, ∞); σ3∼t3(10, 0.1)I(0, ∞) φ3.1 ∼N(0.7, 0.1) η∼Dir(1, 1, 1) ν1∼t3(2.19, 0.1)I(0, ∞); σ1∼t3(3, 0.1)I(0, ∞) φ1.1 ∼N(−0.62, 0.1) ν2∼t3(3.94, 0.1)I(0, ∞); σ2∼t3(6.8, 0.1)I(0, ∞) φ2.1 ∼N(0.2, 0.1) ν3∼t3(9.66, 0.1)I(0, ∞); σ3∼t3(1.62, 0.1)I(0, ∞) φ3.1 ∼N(0.84, 0.1) η∼Dir(1, 1, 1) σ1∼t3(10, 0.1)I(0, ∞) φ1.1 ∼N(−0.6, 0.1) σ2∼t3(11.95, 0.1)I(0, ∞) φ2.1 ∼N(0.24, 0.1) σ3∼t3(31, 0.1)I(0, ∞) φ3.1 ∼N(0.74, 0.1) η∼Dir(1, 1, 1) Econometrics 2021,9, 27 31 of 35 Appendix G.1.3. GMAR Model Econometrics 2021, 9, x FOR PEER REVIEW 32 of 36 Appendix G.1.3. GMAR Model Figure A3. Posterior density plots for some parameters on the GMAR model in the IBM stock prices (first-differenced series). Appendix G.2. Brent Crude Oil Prices (First-Differenced Series) Appendix G.2.1. ZMAR Model Figure A3. Posterior density plots for some parameters on the GMAR model in the IBM stock prices (first-differenced series). Appendix G.2. Brent Crude Oil Prices (First-Differenced Series) Appendix G.2.1. ZMAR Model Econometrics 2021, 9, x FOR PEER REVIEW 32 of 36 Appendix G.1.3. GMAR Model Figure A3. Posterior density plots for some parameters on the GMAR model in the IBM stock prices (first-differenced series). Appendix G.2. Brent Crude Oil Prices (First-Differenced Series) Appendix G.2.1. ZMAR Model Figure A4. Cont. Econometrics 2021,9, 27 32 of 35 Econometrics 2021, 9, x FOR PEER REVIEW 33 of 36 Figure A4. Posterior density plots for some parameters on the ZMAR model in the Brent crude oil prices (first-differenced series). Appendix G.2.2. TMAR Model Figure A4. Posterior density plots for some parameters on the ZMAR model in the Brent crude oil prices (first-differenced series). Appendix G.2.2. TMAR Model Econometrics 2021, 9, x FOR PEER REVIEW 33 of 36 Figure A4. Posterior density plots for some parameters on the ZMAR model in the Brent crude oil prices (first-differenced series). Appendix G.2.2. TMAR Model Econometrics 2021, 9, x FOR PEER REVIEW 34 of 36 Figure A5. Posterior density plots for some parameters on the TMAR model in the Brent crude oil prices (first-differenced series). Appendix G.2.3. GMAR Model Figure A5. Posterior density plots for some parameters on the TMAR model in the Brent crude oil prices (first-differenced series). Econometrics 2021,9, 27 33 of 35 Appendix G.2.3. GMAR Model Econometrics 2021, 9, x FOR PEER REVIEW 34 of 36 Figure A5. Posterior density plots for some parameters on the TMAR model in the Brent crude oil prices (first-differenced series). Appendix G.2.3. GMAR Model Econometrics 2021, 9, x FOR PEER REVIEW 35 of 36 Figure A6. Posterior density plots for some parameters on the GMAR model in the Brent crude oil prices (first-differenced series). References (Al Hakmani and Sheng 2017) Al Hakmani, Rehab, and Yanyan Sheng. 2017. NUTS for Mixture IRT Models. Paper presented at the Annual Meeting of the Psychometric Society, Zurich, Switzerland, July 18–21, pp. 25–37. (Albert and Chib 1993) Albert, James H., and Siddhartha Chib. 1993. Bayes Inference via Gibbs Sampling of Autoregressive Time Series Subject to Markov Mean and Variance Shifts. Journal of Business & Economic Statistics 11: 1–15. (Annis et al. 2017) Annis, Jeffrey, Brent J. Miller, and Thomas J. Palmeri. 2017. Bayesian Inference with Stan: A Tutorial on Adding Custom Distributions. Behavior Research Methods 49: 863–86. (Aroian 1941) Aroian, Leo A. 1941. A Study of RA Fisher’s z Distribution and the Related F Distribution. The Annals of Mathematical Statistics 12: 429–48. (Azzalini 2014) Azzalini, Adelchi. 2014. The Skew-Normal and Related Families. Cambridge: Cambridge University Press. (Box et al. 2015) Box, George E. P., Gwilym M. Jenkins, Gregory C. Reinsel, and Greta M. Ljung. 2015. Time Series Analysis: Forecasting and Control. Hoboken: John Wiley & Sons. (Brown et al. 2002) Brown, Barry W., Floyd M. Spears, and Lawrence B. Levy. 2002. The Log F: A Distribution for All Seasons. Computational Statistics 17: 47–58. (Bruce and Martin 1989) Bruce, Andrew G., and R. Douglas Martin. 1989. Leave-K-Out Diagnostics for Time Series. Journal of the Royal Statistical Society: Series B Methodological 51: 363–401. (Bulmer 1967) Bulmer, Michael George1967. Principles of Statistics. Cambridge: MIT Press. (Carollo 2012) Carollo, Salvatore. 2012. Understanding Oil Prices: A Guide to What Drives the Price of Oil in Today’s Markets. Hoboken: John Wiley & Sons. (Carpenter et al. 2015) Carpenter, Bob, Matthew D. Hoffman, Marcus Brubaker, Daniel Lee, Peter Li, and Michael Betancourt. 2015. The Stan Math Library: Reverse-Mode Automatic Differentiation in C++. arXiv arXiv:1509.07164. (Carpenter et al. 2017) Carpenter, Bob, Andrew Gelman, Matthew D. Hoffman, Daniel Lee, Ben Goodrich, Michael Betancourt, Marcus Brubaker, Jiqiang Guo, Peter Li, and Allen Riddell. 2017. Stan: A Probabilistic Programming Language. Journal of Statistical Software 76: 1–32. (Choir et al. 2019) Choir, Achmad Syahrul, Nur Iriawan, Brodjol Sutijo Suprih Ulama, and Mohammad Dokhi. 2019. MSEPBurr Distribution: Properties and Parameter Estimation. Pakistan Journal of Statistics and Operation Research 15: 179–93. (Fernández and Steel 1998) Fernández, Carmen, and Mark F. J. Steel. 1998. On Bayesian Modeling of Fat Tails and Skewness. Journal of the American Statistical Association 93: 359–71. (Fisher 1924) Fisher, Ronald Aylmer. 1924. On a Distribution Yielding the Error Functions of Several Well Known Statistics. Paper presented at the International Congress of Mathematics, Toronto, ON, Canada, April 11–16, pp. 805–13. (Fong et al. 2007) Fong, Pak Wing, Wai Keung, Li, C. W. Yau, and Chun Shan Wong. 2007. On a Mixture Vector Autoregressive Model. Canadian Journal of Statistics 35: 135–50. (Frühwirth-Schnatter 2006) Frühwirth-Schnatter, Sylvia. 2006. Finite Mixture and Markov Switching Models. Cham: Springer. (Gelman 2006) Gelman, Andrew. 2006. Prior Distributions for Variance Parameters in Hierarchical Models (Comment on Article by Browne and Draper). Bayesian Analysis 1: 515–34. (Gelman and Rubin 1992) Gelman, Andrew, and Donald B. Rubin. 1992. Inference from Iterative Simulation Using Multiple Sequences. Statistical Science 7: 457–72. (Gelman et al. 2014) Gelman, Andrew, John B. Carlin, Hal S. Stern, David B. Dunson, Aki Vehtari, and Donald B. Rubin. 2014. Bayesian Data Analysis, 3rd ed. Boca Raton: CRC Press. Figure A6. Posterior density plots for some parameters on the GMAR model in the Brent crude oil prices (first-differenced series). References Al Hakmani, Rehab, and Yanyan Sheng. 2017. NUTS for Mixture IRT Models. Paper presented at the Annual Meeting of the Psychometric Society, Zurich, Switzerland, July 18–21; pp. 25–37. Albert, James H., and Siddhartha Chib. 1993. Bayes Inference via Gibbs Sampling of Autoregressive Time Series Subject to Markov Mean and Variance Shifts. Journal of Business & Economic Statistics 11: 1–15. Econometrics 2021,9, 27 34 of 35 Annis, Jeffrey, Brent J. Miller, and Thomas J. Palmeri. 2017. Bayesian Inference with Stan: A Tutorial on Adding Custom Distributions. Behavior Research Methods 49: 863–86. [CrossRef] [PubMed] Aroian, Leo A. 1941. A Study of RA Fisher’s zDistribution and the Related F Distribution. The Annals of Mathematical Statistics 12: 429–48. [CrossRef] Azzalini, Adelchi. 2014. The Skew-Normal and Related Families. Cambridge: Cambridge University Press. Box, George E. P., Gwilym M. Jenkins, Gregory C. Reinsel, and Greta M. Ljung. 2015. Time Series Analysis: Forecasting and Control. Hoboken: John Wiley & Sons. Brown, Barry W., Floyd M. Spears, and Lawrence B. Levy. 2002. The Log F: A Distribution for All Seasons. Computational Statistics 17: 47–58. [CrossRef] Bruce, Andrew G., and R. Douglas Martin. 1989. Leave-K-Out Diagnostics for Time Series. Journal of the Royal Statistical Society: Series B Methodological 51: 363–401. Bulmer, Michael George. 1967. Principles of Statistics. Cambridge: MIT Press. Carollo, Salvatore. 2012. Understanding Oil Prices: A Guide to What Drives the Price of Oil in Today’s Markets. Hoboken: John Wiley & Sons. Carpenter, Bob, Matthew D. Hoffman, Marcus Brubaker, Daniel Lee, Peter Li, and Michael Betancourt. 2015. The Stan Math Library: Reverse-Mode Automatic Differentiation in C++. arXiv, arXiv:1509.07164. Carpenter, Bob, Andrew Gelman, Matthew D. Hoffman, Daniel Lee, Ben Goodrich, Michael Betancourt, Marcus Brubaker, Jiqiang Guo, Peter Li, and Allen Riddell. 2017. Stan: A Probabilistic Programming Language. Journal of Statistical Software 76: 1–32. [CrossRef] Choir, Achmad Syahrul, Nur Iriawan, Brodjol Sutijo Suprih Ulama, and Mohammad Dokhi. 2019. MSEPBurr Distribution: Properties and Parameter Estimation. Pakistan Journal of Statistics and Operation Research 15: 179–93. [CrossRef] Fernández, Carmen, and Mark F. J. Steel. 1998. On Bayesian Modeling of Fat Tails and Skewness. Journal of the American Statistical Association 93: 359–71. Fisher, Ronald Aylmer. 1924. On a Distribution Yielding the Error Functions of Several Well Known Statistics. Paper presented at the International Congress of Mathematics, Toronto, ON, Canada, April 11–16; pp. 805–13. Fong, Pak Wing, Wai Keung Li, C. W. Yau, and Chun Shan Wong. 2007. On a Mixture Vector Autoregressive Model. Canadian Journal of Statistics 35: 135–50. [CrossRef] Frühwirth-Schnatter, Sylvia. 2006. Finite Mixture and Markov Switching Models. Cham: Springer. Gelman, Andrew. 2006. Prior Distributions for Variance Parameters in Hierarchical Models (Comment on Article by Browne and Draper). Bayesian Analysis 1: 515–34. [CrossRef] Gelman, Andrew, and Donald B. Rubin. 1992. Inference from Iterative Simulation Using Multiple Sequences. Statistical Science 7: 457–72. [CrossRef] Gelman, Andrew, John B. Carlin, Hal S. Stern, David B. Dunson, Aki Vehtari, and Donald B. Rubin. 2014. Bayesian Data Analysis, 3rd ed. Boca Raton: CRC Press. Geman, Stuart, and Donald Geman. 1984. Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images. IEEE Transactions on Pattern Analysis and Machine Intelligence 6: 721–41. [CrossRef] [PubMed] Gilchrist, Warren. 2000. Statistical Modelling with Quantile Functions. Boca Raton: CRC Press. Hoffman, Matthew D., and Andrew Gelman. 2014. The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo. Journal of Machine Learning Research 15: 1593–623. Huerta, Gabriel, and Mike West. 1999. Priors and Component Structures in Autoregressive Time Series Models. Journal of the Royal Statistical Society: Series B Statistical Methodology 61: 881–99. [CrossRef] Iriawan, Nur. 2000. Computationally Intensive Approaches to Inference in Neo-Normal Linear Models. Bentley: Curtin University of Technology. Johnson, Norman L., Samuel Kotz, and Narayanaswamy Balakrishnan. 1995. Continuous Univariate Distributions, 2nd ed. Hoboken: John Wiley & Sons, vol. 2. Johnson, Norman L., Adrienne W. Kemp, and Samuel Kotz. 2005. Univariate Discrete Distributions, 3rd ed. Hoboken: John Wiley & Sons. Jones, Chris. 2002. Multivariate t and Beta Distributions Associated with the Multivariate F Distribution. Metrika 54: 215–31. [CrossRef] Kim, Hea-Jung. 2008. Moments of Truncated Student-t Distribution. Journal of the Korean Statistical Society 37: 81–87. [CrossRef] Kotz, Samuel, Narayanaswamy Balakrishnan, and Norman L. Johnson. 2000. Continuous Multivariate Distributions, Volume 1: Models and Applications. Hoboken: John Wiley & Sons. Le, Nhu D., R. Douglas Martin, and Adrian E. Raftery. 1996. Modeling Flat Stretches, Bursts Outliers in Time Series Using Mixture Transition Distribution Models. Journal of the American Statistical Association 91: 1504–15. Maleki, Mohsen, Arezo Hajrajabi, and Reinaldo B. Arellano-Valle. 2020. Symmetrical and Asymmetrical Mixture Autoregressive Processes. Brazilian Journal of Probability and Statistics 34: 273–90. [CrossRef] Martin, R. Douglas, and Victor J. Yohai. 1986. Influence Functionals for Time Series. The Annals of Statistics 14: 781–818. [CrossRef] McLachlan, Geoffrey J., and David Peel. 2000. Finite Mixture Models. Hoboken: John Wiley & Sons. Metropolis, Nicholas, Arianna W. Rosenbluth, Marshall N. Rosenbluth, Augusta H. Teller, and Edward Teller. 1953. Equation of State Calculations by Fast Computing Machines. The Journal of Chemical Physics 21: 1087–92. [CrossRef] Neal, Radford M. 2011. MCMC Using Hamiltonian Dynamics. In Handbook of Markov Chain Monte Carlo. Edited by Steve Brooks, Andrew Gelman, Galin L. Jones and Xiao-Li Meng. Boca Raton: CRC Press, pp. 116–62. Econometrics 2021,9, 27 35 of 35 Nguyen, Hien D., Geoffrey J. McLachlan, Jeremy F.P. Ullmann, and Andrew L. Janke. 2016. Laplace Mixture Autoregressive Models. Statistics & Probability Letters 110: 18–24. Pravitasari, Anindya Apriliyanti, Nur Iriawan, Kartika Fithriasari, Santi Wulan Purnami, and Widiana Ferriastuti. 2020. A Bayesian Neo-Normal Mixture Model (Nenomimo) for MRI-Based Brain Tumor Segmentation. Applied Sciences 10: 4892. [CrossRef] Solikhah, Arifatus, Heri Kuswanto, Nur Iriawan, Kartika Fithriasari, and Achmad Syahrul Choir. 2021. Extending Runjags: A Tutorial on Adding Fisher’sz Distribution to Runjags. AIP Conference Proceedings 2329: 060005. Stan Development Team. 2018. Stan User’s Guide, Version 2.18.0. Available online: https://mc-stan.org/docs/2_18/stan-users-guide/ index.html (accessed on 21 October 2020). Stan Development Team. 2020. RStan: The R Interface to Stan. Available online: https://cran.r-project.org/web/packages/rstan/ vignettes/rstan.html (accessed on 21 October 2020). Susanto, Irwan, Nur Iriawan, Heri Kuswanto, Kartika Fithriasari, Brodjol Sutija Suprih Ulama, Wahyuni Suryaningtyas, and Anindya Apriliyanti Pravitasari. 2018. On The Markov Chain Monte Carlo Convergence Diagnostic of Bayesian Finite Mixture Model for Income Distribution. Journal of Physics: Conference Series 1090: 012014. [CrossRef] Vehtari, Aki, Andrew Gelman, and Jonah Gabry. 2017. Practical Bayesian Model Evaluation Using Leave-One-out Cross-Validation and WAIC. Statistics and Computing 27: 1413–32. [CrossRef] Vehtari, Aki, Andrew Gelman, Daniel Simpson, Bob Carpenter, and Paul Christian Bürkner. 2020. Rank-Normalization, Folding, and Localization: An Improved Rhat for Assessing Convergence of MCMC. Bayesian Analysis 1: 1–28. [CrossRef] Wong, Chun Shan, and Wai Keung Li. 2000. On a Mixture Autoregressive Model. Journal of the Royal Statistical Society: Series B Statistical Methodology 62: 95–115. [CrossRef] Wong, Chun Shan, and Wai Keung Li. 2001. On a Logistic Mixture Autoregressive Model. Biometrika 88: 833–46. [CrossRef] Wong, Chun Shan, Wai-Sum Chan, and P. L. Kam. 2009. A Student T-Mixture Autoregressive Model with Applications to Heavy-Tailed Financial Data. Biometrika 96: 751–60. [CrossRef] World Bank. 2020. World Bank Commodity Price Data (The Pink Sheet). Available online: http://pubdocs.worldbank.org/en/561011 486076393416/CMO-Historical-Data-Monthly.xlsx (accessed on 21 October 2020). Zelen, Marvin, and Norman C. Severo. 1970. Probability Functions. In Handbook of Mathematical Functions; Edited by Milton Abramowitz and Irene A. Stegun. Washington, DC: U.S. Government Printing Office, pp. 925–95.