Fisher’s z Distribution-Based Mixture Autoregressive Model

: 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 classiﬁed 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.


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 where −∞ < x < ∞, −∞ < µ < ∞, and σ > 0. Equation (3) is defined as a p.d.f of Fisher's z distribution. It is denoted as z(d 1 , d 2 , µ, σ). The CDF of the Fisher's z distribution is expressed as where . The quantile function (QF) of the Fisher's z distribution is defined as where d 2 I −1 x p 1 2 d 1 , 1 2 d 2 / d 1 1 − I −1 x p 1 2 d 1 , 1 2 d 2 is the QF of the F-distribution and I −1 x p (.) is the inversion of the incomplete beta function ratio. Let 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 where 2 P −1 v 1p 1 2 d 1 and 2 P −1 v 2p 1 2 d 2 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 = 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 = µ. 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. = + 2 ln 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.   = + 2 ln 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.

Model Specification
Let y t ; 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 y t F (t−1) ; k = 1, 2, · · · , 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 Z (d 1k ,d 2k ) (.) 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 {y t−i , i > 0}; and let σ k is a scale parameter of the kth component. The K-component ZMAR can be defined as or with where the vector η = (η 1 , . . . , η K ) is called the weights. η takes a value in the unit simplex E K , which is a subspace of ( + ) K , defined by the following constraint, η k > 0 and We use the abbreviation ZMAR(K; p 1 , p 2 , . . . , p K ) for this model, with the parameter denotes the AR coefficient on the kth component and ith lag; i = 1, 2, . . . , p k ; and p k denotes the autoregressive order of the kth component. Using the parameters θ k , we first define the K auxiliary Fisher's z AR(p k ) processes where the AR coefficients φ k are assumed to satisfy This condition implies that the processes f k y t F (t−1) are stationary and that each component model in (7) or (8) satisfies the usual stationarity condition of the linear AR(p k ) model.
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 σ k , with the degrees of freedoms v 1k , v 2k ,ν 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 . 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 u k.0 , u k.i and the scale parameters g k.0 2 , g k.i 2 , thus φ k.0 ∼ N u k.0 , g k.0 2 and φ k.i ∼ N u k.i , g k.i 2 ; i = 1, 2, . . . , 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 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 where i = 1, 2, . . . , p k ; 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.
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(Carpenter et al. , 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.

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 3: Step 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 pk = p k ; con_eta = c(1, 1, 1); tau [t, k] 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))); tential scale reduction factor (Gelman and Rubin 1992;Susanto et al. 2018;Gelman et al. 2014, 285;Vehtari et al. 2020) and the effective sample size (Gelman et al. 2014, 266;Vehtari et al. 2020). If the MCMC chain has reached convergence, the statistic is less than 1.01, and the statistic is greater than 400 (Vehtari et al. 2020). To compare the performance of the models, we use the PSIS-LOO.  Table 1 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. mparison using pareto-smoothed important sampling leave-one-out cross-validation (PSIS-LOO) in sev-  In the third and sixth scenarios, the datasets are generated from three symmetrical 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 factorR (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, theR 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 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.

ZMAR
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. Simulation Model:

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 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 Appendices F.1 and G.1.1, respectively. For all the parameters of the ZMAR model, the MCMC chain has reached convergence, which is shown by theR statistic being less than 1.01 and the n e f f 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 6.41 6.09 6.75 We compared the ZMAR model with the TMAR and GMAR models. The best threecomponent 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 Appendices F.1 and 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 theR statistic being less than 1.01 and the n e f f statistic being greater than 400. Let Φ(.), and let T ν k (.) be the CDF of the standard normal distribution and the standardized Student t distribution with ν k ; k = 1, 2, 3 degrees of freedom. The three-component TMAR model for the original series, namely +0.42 Φ y t −1.67 y t−1 +0.67 y t−2 6.01 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.

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. 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 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, Econometrics 2021, 9, x FOR PEER REVIEW in the world, the other being West Texas Intermediate (WTI). Figure 5a shows the m Brent crude oil price from January 1989 to June 2020 (in U.S. dollars per barrel), tak the World Bank (2020). Figure 5b shows the first-differenced series. Figure 6 sh marginal distribution of the original series as being trimodal, where the estimat tions of the modes are at 18.57, 61.64, and 110.15 points. Therefore, we also dec choose a three-component mixture model for the ZMAR, TMAR, and GMAR mo plied to the differenced series. The best three-component mixture model estim each of the ZMAR, TMAR, and GMAR models are ZMAR (3;1,2,2), TMAR (3;2,1 GMAR (3;3,3,3), without intercept.  The warm-up steps for all the models were set to 1500 iterations with 5000 sa iterations and 1 thin, followed by 3 chains; the adapt_delta parameters were set and the max_treedepths were set to 15. Table 3 shows the summary of posterior inf for all models. The prior distributions and posterior density plots for the param 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 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 Appendices F.2 and G.2. For all the parameter models, the MCMC chains reached convergence, which was shown by theR statistic being less than 1.01 and the n e f f statistic being greater than 400.
The three-component ZMAR model for the original series, namely 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.

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 µ, σ, d 1 , and d 2 . The µ is a location parameter, the σ 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 = 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 = µ. The value of µ 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 µ; if d 1 < d 2 , then the mean is less than µ; and if d 1 > d 2 , 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 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.
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.

Informed Consent Statement: Not applicable.
Data Availability Statement: Data available from stated sources.

Conflicts of Interest:
The authors declare no conflict of interest.

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 F distribution to the Fisher's z distribution. Let Y be a random variable distributed as an F distribution with two parameters d 1 and d 2 . The density of Z = 1 2 ln Y is (Fisher 1924;) where −∞ < z < ∞, d 1 > 0, d 2 > 0 and B(.) is the beta function. Interchanging d 1 and d 2 is equivalent to replacing z with −z (Fisher 1924;, thus Equation (A1) can also be defined as If the denominator and numerator of Equation (A2) are divided by d 1 (d 1 +d 2 )/2 then we get Equation (A3) can also be defined as Appendix A.2. Proof of the Equation (2) Let Y be a random variable distributed as an F distribution with d 1 and d 2 degrees of freedom, let I y * (.) be the incomplete beta function ratio, and let B(.) be the beta function. The cumulative distribution function (CDF) of the Y = e 2z be defined as follows (Johnson et al. 1995, vol. 2, p. 327) where y * = d 1 y d 2 +d 1 y . If Z = 1 2 ln Y then Y = e 2z , thus the CDF of Z is where z * = d 1 e 2z d 2 +d 1 e 2z . Equation (A4) can also be defined as where z * = d 2 e −2z d 1 +d 2 e −2z .

Appendix A.3. Proof of the Equation (3)
Let Z be a random variable distributed as a standardized Fisher's z distribution 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 where J(x) is the Jacobian of the transformation and is defined as Therefore, the p.d.f of the Fisher's z distribution can be expressed as Appendix A.4. Proof of the Equation (4) Let Z be a random variable distributed as a standardized Fisher's z distribution 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 Therefore, the CDF of the Fisher's z distribution can be expressed as .
Appendix A.5. Proof of the Equation (5) Let I x * (.) be the incomplete beta function ratio; thus, Equation (4) can be expressed as The value x p is called the p-quantile of the population, if P X ≤ x p = p with 0 ≤ p ≤ 1 (Gilchrist 2000, p. 12). Let I −1 x p (.) be the inversion of the incomplete beta function ratio, then Interchanging d 1 and d 2 is equivalent to replacing x with −x; thus, the QF can also be defined as: Appendix A.6. Proof of the Equation (6) Let us denote a chi-square random variables with d 1 and d 2 degrees of freedom by V 2 1 and V 2 2 , respectively. Let P v * (.) be the incomplete gamma function ratio, and the CDF of V 2 1 and V 2 2 can be defined as Let P −1 v p (.) be the inversion of the incomplete gamma function ratio, and then the QF of the V 2 1 and V 2 2 can be defined as The Beta distribution arises naturally as the distribution of X = V 2 1 / V 2 1 + V 2 2 (Johnson et al. 1995, vol. 2, p. 212); therefore, QF of the Fisher's z distribution can be expressed as

Appendix B.1. Properties of the Fisher's z Distribution
Let X be a random variable distributed as a Fisher's z distribution and let M X (θ) be the Moment Generating Function (MGF) of a random variable X. The MGF of the Fisher's z distribution be expressed as where Γ(.) is a gamma function. Let K X (θ) be the Cumulant Generating Function (CGF) of a random variable X. The CGF of the Fisher's z distribution be given by (A7) The coefficient of θ!/j! in the Taylor expansion of the CGF is the jth cumulant of X and be denoted as κ j . The jth cumulant, therefore, can be obtained by differentiating the expansion j times and evaluating the result at zero.
; j = 1, 2, 3, · · · The first cumulant of the Fisher's z distribution be defined as and the jth cumulant would be defined in Equation (A9).
where ψ(.) is the digamma function, ψ (.) is the trigamma function, ψ (.) 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 X be denoted respectively E(X) and Var(X). The mean of the Fisher's z distribution is given by and the variance is defined as On the basis of Equation (A10), it can be concluded that Let the skewness and the excess kurtosis of a random variable X be denoted, respectively, as γ 1 and γ 2 . The skewness of the Fisher's z distribution is given by and the excess kurtosis is 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 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.   Figure A1. Posterior density plots for some parameters on the ZMAR model in the IBM stock prices (first-differenced series).
Appendix G.1.2. TMAR model Figure A2. Posterior density plots for some parameters on the TMAR model in the IBM stock prices (first-differenced series). Appendix G.1.2. TMAR model Figure A2. Posterior density plots for some parameters on the TMAR model in the IBM stock prices (first-differenced series). Figure A2. Posterior density plots for some parameters on the TMAR model in the IBM stock prices (first-differenced series).