Determining the Order of a Moving Average Model of Time Series Using Reversible Jump MCMC: A Comparison between Laplacian and Gaussian Noises

Moving average (MA) is a time series model often used for pattern forecasting and recognition. It contains a noise that is often assumed to have a Gaussian distribution. However, in various applications, noise often does not have this distribution. This paper suggests using Laplacian noise in the MA model, instead. The comparison of Gaussian and Laplacian noises was also investigated to ascertain the right noise for the model. Moreover, the Bayesian method was used to estimate the parameters, such as the order and coefficient of the model, as well as noise variance. The posterior distribution has a complex form because the parameters are concerened with the combination of spaces of different dimensions. Therefore, to overcome this problem, the Markov Chain Monte Carlo (MCMC) reversible jump algorithm is adopted. A simulation study was conducted to evaluate its performance. After it has worked properly, it was applied to model human heart rate data. The results showed that the MCMC algorithm can estimate the parameters of the MA model. This was developed using Laplace distributed noise. Moreover, when compared with the Gaussian, the Laplacian noise resulted in a higher order model and produced a smaller variance.


Introduction
Moving average (MA) is often used to model data in various fields of life. Studies on its application is seen in several literatures, such as [1], [2], [3], [4], and [5]. In fact, Reghunath et al. in [1] used MA to analyze water resource data. The model was used to determine the long-term trends in groundwater fluctuations. To eliminate seasonal effects on this fluctuation, the 12th order MA was used to model water source data. Furthermore, Akrami et al. in [2] used it to analyze rainfall data. To observe long and short-term trends, several MA models with different orders were examined. Silva de Souza et al. in [4] used it for profitability checks in technical analysis on the stock market. An automated trading system was developed based on moving averages from previous prices. Meanwhile, Gautam and Abhishekh in [5] used the model to analyze the data of fuzzy time series by reducing fluctuations.
Researchers like in [6] and [7] used the Gaussian distributed noise. Meanwhile, Middleton et al. in [6] used the type that contained Gaussian to model the noise channel. Safi in [7] examined the heteroscedastic autocorrelation function of residue in MA. Therefore, in determining this functions, the model was assumed to have Gaussian noise.
In various applications, mathematical models often 722 Determining the Order of a Moving Average Model of Time Series Using Reversible Jump MCMC: A Comparison between Laplacian and Gaussian Noises have non-Gaussian noise, such as [8], [9], [10], and [11]. Jureckova et. al. in [8] examined non-parametric tests in linear autoregressive models that use Pareto noise. Also, Zhang in [9] examined the autoregressive type that uses G-GARCH as noise. Suparman in [10] developed a stationary model that has exponentially distributed noise. In addition, an inversible moving range type was developed, which has exponentially distributed noise. In this autoregressive and moving average model, the order is assumed to be unknown, hence the posterior distribution has a complicated shape. To determine the Bayesian estimator, the research adopted the MCMC reversible jump algorithm.
Conversely, Laplace distributed noise was also investigated by several researchers, such as [12] and [13]. Minchole et. al. in [12] detected that changes in body position from the ECG use Laplacian noise. Also, Miertoiu and Dumitrescu in [13] investigated the signal representation that rarely uses Laplacian noise.
However, MA with Laplacian has not been studied. Therefore, this research proposed the development of the model with Laplacian, which is then used for the human heart rate data. The Laplacian was also compared with Gaussian to determine the effect of noise on MA models.

Background and Methodology
The Bayesian method was used to estimate the MA model parameters, such as its order, coefficients, and noise. In this method, the parameters were treated as random variables that are assumed to have a certain distribution, known as "prior." This was combined with the likelihood function for data, which produced a posterior distribution. In this case, the order of the MA model is also a parameter. Therefore, the shape of posterior distribution becomes very complex, which makes it difficult to explicitly determine the Bayes estimator. To determine this, the Markov chain Monte Carlo (MCMC) reversible jump algorithm was adopted [14]. The basic idea of this was to create a Markov chain. This chain was designed, hence its limit distribution approached the posterior for the parameters. Also, the algorithm uses 3 transformations, such as birth and death of the order, as well as coefficient change. The performance of this algorithm was evaluated using simulation studies and it was applied to model human heart rate data.

Result and Discussion
This section explained the functions of data, posterior and prior distribution, as well as MCMC reversible jump in determining Bayes estimators.

Likelihood Function
For example, 1 , … , are data and n is the amount.
Therefore, it is said to be modeled by the order MA, written as ( ), when it satisfies the equation The magnitudes of 1 , … , express the coefficients of the ( ) model. Meanwhile, the variables of ( = 1, … . , ) express noise. This was assumed to be a Laplace distribution with parameter. Therefore, the probability function for is written as follows Furthermore, the probability function for can be determined using the variable transformation between and . The function of can be stated as follows For example, = ( 1 , … , ) and ( ) = � 1 , … , �. Hence, the function can be stated as follows The amount of is calculated using the following equation for = + 1, … , and = 0 for t = 1, 2, … ,q. The ( ) model is inversible when the coefficients of 1 , … , fulfill the inversibility area where For large order values, the area is very difficult to identify. Therefore, a transformation is needed to overcome this problem. For example, 1 , 2 , ⋯ , express the inverse partial autocorrelation functions that correspond to the is the inversibility area [15]. Therefore, this transformation can easily be used to identify areas of inversibility even for large orders. The ( ) model is inversible only when� 1 , 2 , ⋯ , � ∈ (−1,1) .Also, the function for can be rewritten as follows where −1 is the inverse transformation of .

Prior Distribution
The MA parameter is seen as a variable with prior distribution. The distribution for this model order is binomial with and parameters.
The value of is limited to order 10. Meanwhile, the binomial distribution was chosen because it is a conjugate. It means that if the beta distribution is used as a prior, a posterior of binomial likelihood will also be the beta distribution. For example, ( ) = � 1 , 2 , ⋯ , �. The prior or ( ) given the value of order is a uniform distribution at intervals (−1, 1) Finally, the prior distribution for is the Gamma inverse with parameters and where the value of is set to egual 1. The prior distribution for MA model order contains hyperparameter and that for contains . The hyperprior distribution for is uniform at intervals (0.1). Meanwhile, is of Jeffreys ( ) ∝ 1 . Therefore, the combined distribution of prior and hyperprior can be written as follows

Posterior Distribution
By using the Bayes theorem, the posterior distribution for � , ( ) , , , � can be stated as follows � , ( ) , , , Parameter space is a combination of those with different dimensions. This makes it difficult to analytically determine the Bayes estimator. Therefore, the MCMC reversible jump algorithm was used to achieve this.

Reversible Jump MCMC
Simulation of posterior distribution for parameters was performed in 2 stages, which are conditional simulation for( , , ) when given � , ( ) � and for � , ( ) � when given ( , , ). Meanwhile, the conditional distribution for ( , , ) when given � , ( ) � is a multiplication results of the Gamma inversion, Binomial, and the Gamma distribution. Therefore, the conditional simulation for ( , , ) when given � , ( ) � can be performed as follows: However, conditional distribution for � , ( ) � with ( , , ) has a complex form. Therefore, the simulation for � , ( ) � with ( , , ) is performed using the MCMC algorithm. This uses 3 types of transformation, such as coefficient change, birth, and death of the order [10].

Coefficient Change
For example, = � , ( ) � is the old Markov Chain and * = � * , * ( * ) � is the new. The change in the coefficient of MA model does not alter the order value. Meanwhile, the alteration from to * was performed in two steps. Firstly, value * =q was taken, secondly, ∈ {1, … , } was chosen and * = where ~(−1,1) were defined. The ratio between the likelihood function ( | * ) and ( | ) can be stated as follows: The ratio between prior distribution for * and for is Meanwhile, the ratio between conditional prior for * when given * and conditional prior distribution for with is ( * | * ) ( | ) = 1. Also, the ratio between instrument distribution ( * , ) and ( , * ) is: ( , * ) is the acceptance probability for the coefficient change, which is written as follows:

Birth of the order
For example, = � , ( ) � is the old Markov chain and * = � * , * ( * ) � is the new. The birth of MA alters the order and coefficient value of the model. Meanwhile, change from to * was performed in two steps. Firstly, * = + 1 was taken, secondly, * = ( , ) where ~(−1,1) was defined. The ratio between the function ( | * ) and ( | ) can be stated as follows: The ratio between the prior distribution for * and p for is Also, the ratio between conditional prior distribution for * when given * , and for with is ( * ) ( ) = 1. This can be written as: The ratio between the distribution of the instrument ( * , ) and ( , * ) depending on the value of when < 0 is: Meanwhile, when > 0 , the ratio between the distribution of instrument ( * , ) and ( , * ) is: ( , * ) is the acceptance probability for birth of order, which can be written as follows:

Death of Order
The death of the MA model order was different from the transformation. For example, = � + 1, ( +1) � is the old markov chain and * = ( , ) is the new. The birth changes the order and coefficient value of the model. Also, alterations from to * are made in two steps. Firstly, * = was taken, and secondly, * = ∖ � +1 � was defined.
( , * ) is the acceptance probability for the death of the order, which can be written as follows:

Simulation
The performance of MCMC reversible jump algorithm was tested using simulation study. This addresses two MA models. The first assumes that noise is Laplacian distributed. Meanwhile, the second model assumes that noise is Gaussian.

Laplacian Noise
For the first simulation, 250 synthesis time series were made according to equation (1) with Laplacian distributed noise. This was given in equation (2). Meanwhile, the MA model parameters are stated in Table 1. The maximum order is limited to order 10, = 10. The resulting time series is presented in Figure 1. Furthermore, the parameters were estimated based on the synthesis time series. The MCMC reversible jump algorithm was used to estimate the parameters. Also, it was run in 20,000 iterations with a burn-in period of 5000. The histogram for the order is presented in Figure 2.  Figure 2 showed that the maximum frequency was reached in order 3. Therefore, the estimator for is � = 3. Given � = 3, the results of the model coefficients are presented in Table 2. Comparing Table 1 and 2, the estimated values of the parameters are close to the actual values. The simulation showed that the algorithm can estimate the parameters that have Laplacian distributed noise.

Gaussian Noise
For the second simulation, 250 synthesis time series were made according to equation (1) but the noise is Gaussian distributed, which is obtained by the following The MA model parameters are shown in Table 3. As in the first simulation, the order is limited to order 10, = 10.
The resulting time series is presented in Figure 3. Furthermore, the parameters are also estimated based on the synthesis time series. The MCMC was also used to estimate the parameters. As in the first simulation, the algorithm was also run in 20,000 iterations with a burn-in period of 5000. The histogram for the order is presented in Figure 4.  Figure 4 showed that the maximum frequency was reached in order 2. Therefore, the estimator for is � = 2. Given � = 2, the results of the model coefficients are presented in Table 4. As in the first simulation, comparing the Table 3 and 4, the estimated values of the parameters are close to the actual. This simulation study showed that the MCMC algorithm can estimate the parameters of the MA model that have Gaussian distributed noise.

Heart Rate Data
Heart rate is the number of cardiac cycles per minute. The heart is an organ that responds in the same manner as other excitable tissues. In fact, it hypertrophy and becomes stronger when a person exercises. In the absence of exercise, this organ continues to pump blood throughout the body to facilitate tissue repair and recovery. Studies on heart rate can be found in a variety of literatures, such as [16] and [17], and [18].
In this paper, the authors use MA to model human heart rate data. The number of beats per minute was recorded for 100 minutes. The recording data are presented in Figure 5. In this first case, this rate was modeled by MA with Laplacian distributed noise. Meanwhile, the MCMC algorithm was implemented to estimate the parameters with Laplacian noise. The algorithm was run in 100000 iterations with a burn-in period of 10000. Figure 6 showed a histogram for the order.   Figure 6 showed that the maximum frequency was reached in order 10. Therefore, the estimator for is � = 10. Given q = 2, the results of the model coefficients are presented in Table 5. In this second case, the rate was modeled with the Gaussian distributed noise. The MCMC was implemented to estimate the parameters. Furthermore, the algorithm was run in 100000 iterations with a burn-in period of 10000. Figure 6 showed a histogram for the order.  Figure 7 showed that the maximum frequency was reached in order 10. Therefore, the estimator for is � = 7. Given � = 7, the results of the model coefficient are presented in Table 6. When the model with Laplacian noise was compared to the MA with Gaussian, the Laplacian noise resulted in a higher order model. Conversely, it produced a model with smaller noise variance.
As in [19], research related to the MA model can be extended into a piecewise constant MA model. The MA model is a special case of the piecewise constant MA model when the number of MA models is only one.

Conclusions
This research resulted in the development of MA, which is the model with Laplacian noise. Also, the Bayesian approach was used to estimate the model, and the order was assumed as a parameter. Furthermore, the MCMC reversible Jump algorithm was applied to generate a Markov chain. This chain was designed, and the limit distribution approached the posterior for the parameter. This was then used to determine the Bayes estimator.
Simulation study showed that the MCMC algorithm can estimate the order and coefficient of the MA model, and noise variance. Also, the human heart rate was modeled by MA with Gaussian and Laplacian noises. This was implemented to ascertain a suitable MA for the heart rate. Laplacian noise produced a model that has a higher order than Gaussian noise. The advantage of this is that those with Laplacian has a smaller noise variance than a model with Gaussian.