Bayesian autoregressive adaptive refined descriptive sampling algorithm in the Monte Carlo simulation

This paper deals with the Monte Carlo Simulation in a Bayesian framework. It shows the importance of the use of Monte Carlo experiments through refined descriptive sampling within the autoregressive model , where and the errors are independent random variables following an exponential distribution of parameter θ. To achieve this, a Bayesian Autoregressive Adaptive Refined Descriptive Sampling (B2ARDS) algorithm is proposed to estimate the parameters ρ and θ of such a model by a Bayesian method. We have used the same prior as the one already used by some authors, and computed their properties when the Normality error assumption is released to an exponential distribution. The results show that B2ARDS algorithm provides accurate and efficient point estimates.


Introduction
Time series are often used to model data in different areas of life. We can think of economic forecasts, the evolution of a biological system, time-stamped data from social media, E-commerce and sensor systems. The AutoRegressive (AR) model is the most widely used class of time series models, and particularly the first order AutoRegressive process denoted as AR(1).
Several works related to time series models have been treated by the Bayesian approach. Box et al. (1976), Monahan (1984) and Marriot and Smith (1992) discussed the Bayesian approach to analyse the AR models. Diaz and Farah (1981) devoted their work to a Bayesian technique for computing posterior analysis of an AR process in any arbitrary order. Amaral Turkmann (1990) considered a Bayesian study of an AR(1) model with exponential white noise based on a noninformative prior. Ghosh and Heo (2003) introduced a comparative study to some selected noninformative priors for the AR(1) model. Pereira and Amaral-Turkman (2004) developed a Bayesian analysis of a threshold AR model with exponential white noise.
Several Bayesian methods for estimating autoregressive parameters have been proposed by different authors; for instance, Albert and Chib (1993) and Barnett et al. (1996) proposed a Bayesian method via Gibbs sampling to estimate the parameters in an AR model. Ibazizen and Fellag (2003) proposed a Bayesian estimation of an AR(1) process with exponential white noise by using a noninformative prior, whereas Suparman and Rusiman (2018) proposed a hierarchical Bayesian estimation for stationary AR models using a reversible jump MCMC algorithm. However, Kumar and Agiwal (2019) proposed a Bayesian estimation for a full shifted Panel AR(1) time series model. Bayesian inference is based on the posterior distribution, which can be interpreted as the result of the combination of two information sources: the information provided by the parameter which is contained in the prior distribution, and the information provided by the data which formalizes the likelihood function through the use of a given statistical model.
Posterior inference can be undertaken by applying classical Monte Carlo (MC) methods, which provide iterative algorithms that can generate approximate samples from the posterior distribution. These methods fall into two categories: the Markov Chain Monte Carlo (MCMC) methods such as Metropolis-Hastings and Gibbs algorithms (Smith & Robert, 1993), and the Monte Carlo (MC) method also known as Simple Random Sampling (SRS), for Table 1. Subsets of prime numbers size for a Gumbel distribution with mean E(X) = 0.577, p 1 = 7, p 2 = 11 and p 3 = 13. instance, Quasi Monte Carlo (Keller et al., 2008), Latin hypercube sampling (LHS) (Loh, 1996;Mckay et al., 1979;Stein, 1987), Descriptive Sampling (DS) (Saliby, 1990) as well as Refined Descriptive Sampling (RDS) (Tari & Dahmani, 2006). In the non-Bayesian estimation, the SRS is a widely used method in MC simulation. But unfortunately, the sampling method generates estimation errors, so it is not accurate enough. However, the RDS method (Aloui et al., 2015) has shown its effectiveness for several models and in various fields. For instance, in the field of medicine, Ourbih-Tari and Azzal (2017) proposed an RDS algorithm for a nonparametric estimation of the survival function using Kaplan Meier and Fleming Harrington estimators in a case study. In the field of engineering, Baghdali-Ourbih et al. (2017) have estimated the performance measures of the aggregate quarry of Bejaia city in Algeria, using a discrete-event simulation model. In the field of M/G/1 queuing models, M/M/1 retrial queues and M/G/1 retrial queues, simulation results using RDS outperform those obtained using the MC method (Idjis et al., 2017;Tamiti et al., 2018). Also, a comparison between RDS and LHS carried out on M/G/1 queues showed the efficiency of RDS over LHS (Boubalou et al., 2019). Since it is important to choose a sampling scheme that allows for a better parameter estimate while keeping the number of repetitions to a minimum, the RDS method is chosen for further work. This paper will therefore explore the RDS method in a Bayesian sense by considering the AR(1) model, which has already been studied by Amaral Turkmann (1990) and, later on, by Ibazizen and Fellag (2003). For this purpose, we assume the noninformative prior given in Ibazizen and Fellag (2003) for the AR(1) process with exponential white noise. Then, we develop an algorithm for the estimation of AR(1) model parameters in a Bayesian context by using the most accurate sampling method as the RDS method. The proposed algorithm is called Bayesian Autoregressive Adaptive Refined Descriptive Sampling (B2ARDS).
Some details of the RDS method are given in Section 2. Then, in Section 3, the methodology from Bayesian perspective is explained. Section 4 is concerned with the Bayesian estimates of the AR(1) model parameters and their variances using the RDS method through the B2ARDS algorithm. Section 5 presents simulation experiments performed by using the R Package and corresponding results are summarized in Tables 1 and 2, and Figures 1 and 2. Finally, Section 6 concludes this paper.

The use of RDS
DS is known to have two problems. From a theoretical perspective, DS can be biased. More practically, its strict operation requires a prior knowledge of the required sample size. Based on this, RDS has been proposed to make DS safe, efficient and convenient. RDS makes DS safe by reducing substantially the risk of sampling bias, efficient by producing estimates with lower variances and convenient by removing the need to determine in advance the sample size. Not only this, but RDS beats DS in several comparisons and such is the case, on a flow shop system (Tari & Dahmani, 2005a) and a production system (Tari & Dahmani, 2005b). So, RDS is the preferred method when it comes to performing MC simulation.

Description of the method
RDS is based on a block of sub-samples of randomly generated prime number sizes, so it distributes RDS sample values from sub-samples at the request of the simulation. We stop the process of the generation of sample values when the simulation terminates, say when m prime numbers p j , j = 1, 2, . . . , m, have been used, which derives m sub-runs.
The choice of primes in the RDS method as sub-sample sizes is motivated by the definition of primes since we must avoid integers having multiples, whether they are even or odd numbers. Indeed, if the underlying frequency is different from the sampling frequency (the product of primes) or a multiple of it, then the bias of the estimates will be insignificant as demonstrated in Tari and Dahmani (2006). This should ensure that the method reduces the sampling bias.
Remark: In the case of M replicated simulation runs, the RDS method will distribute M blocks defined by m 1 , m 2 , . . . , m M primes which are not the same for all replicated simulation runs.

RDS sample values
Formally, in an RDS run defined by a sample of size n = m j=1 p j , the sample values drawn from the input realvalued random variable X having F as the cumulative distribution are given by the following formulas: where for any j = 1, 2, . . . , m, the values of each sub-sample (R j i , i = 1, 2, . . . , p j ) are randomly selected without replacement from the regular numbers of the following sub-set and such r j i is the midpoint of the p j subintervals in which the interval [0, 1] is subdivided. The primes p j , j = 1, 2, . . . , m, are randomly generated between 7 and a given value fixed by the user.

The estimator of RDS
The j-th sub-run generating p j values of the input random variable X leads to the following estimate of the parameter θθ Therefore, in a given run, the use of refined descriptive sample of size n = m j=1 p j leads to the following sampling estimate of θ defined by the average of those estimates observed on different sub-runs: where h is the simulation function.

Example
We illustrate the method by taking the Gumbel distribution and the prime numbers p 1 = 7, p 2 = 11 and p 3 = 13. The values of the variable are simulated by the formula x i = − ln(− ln(r i )), i = 1, 2, . . . , p j and j = 1, 2, . . . , m. When m = 3, the obtained sub-sets and sub-samples are given in Table 1 using the formulas (2) and (1), respectively, together with the observed mean (θ j ) in each sub-run using the formula (3), while the overall mean computed at the end of the run is given in Table 2 using the formula (4).

Monte Carlo methods from a Bayesian perspective
We consider a random vector Y having a probability density function f (y | θ), where θ is a vector of unknown parameters. The Bayesian statistic considers θ as a random vector of density π(θ) called prior law of θ.
In a Bayesian sense, any inference is based on the computation of the posterior distribution which is the conditional distribution of the random vector θ knowing y and is denoted by its probability density π(θ | y). Bayesian estimator of the parameter θ is constructed from the posterior distribution π(θ | y) by minimizing an appropriate loss function. Most Bayesian estimators have the form of E π [h(θ )], where E π denotes the mathematical expectation of the posterior distribution π(θ | y).
A major difficulty in an estimation by the Bayesian approach, is the explicit computation of the posterior distribution as soon as the number of parameters exceeds two. Also computing E π [h(θ )] can be very difficult even if the law is known. Thus, several approximation methods such as the MC and MCMC methods are proposed in the literature. The latter allows for the estimation of moments of any order of random variables, too complex to be studied analytically, using iterative generation algorithms.
The Monte Carlo method is a numerical method of solving problems using random numbers. It is one of the most popular mathematical tools with a very wide scope, like differential equation integration, fluid mechanics equations, financial mathematics, particles transport and matrix inversion.
In a Bayesian sense, the Monte Carlo method consists in estimating where θ 1 , θ 2 , . . . , θ n are independent and identically distributed i.i.d. random variables according to π(θ | y). The strong law of large numbers ensures thatÎ h converges almost surely to I h when n tends to infinity.

Presentation of the model
Consider the first order autoregressive process defined by where 0 < ρ < 1 and the Y t 's are i.i.d. random variables with the exponential distribution exp(θ ) of parameter θ having a density f (y) = θ exp(−θ y)I (0, ∞) (y), θ > 0.
X 1 is assumed to be distributed according to exp((1 − ρ)θ) such that the process is mean stationary. The likelihood function based on the observations The Maximum Likelihood Estimators (MLE) of ρ and θ are introduced by Andel (1988) , respectively as follows: The model given by (5) was processed using a Bayesian approach by Amaral Turkmann (1990). The author's study is based on a noninformative prior defined as follows: The author derived the Bayes estimators of the independent parameters ρ and θ as follows: where r = 1 − ρ 0 ( S nx ) such that 0 < r < 1. Ibazizen and Fellag (2003) proposed a Bayesian analysis of the model (5) by using a more general noninformative prior defined by which includes the formula (9) as a special case (β = 1). For θ, the prior (7) is proportional to the scale parameter 1/θ and for ρ, it is a special case of the beta distribution. The authors have derived Bayesian estimators given below, but their formulas are complicated for computation. So, simulation experiments were carried out and the obtained results were compared to the MLE.
is the Gauss hypergeometric function of parameters (a, b, c). We note that to compute the latter, the Pochhammer coefficients (t, m) for t = a, b and c are defined as follows where and B represent symbols of the usual functions Gamma and Beta, respectively. This paper focuses on the Bayesian estimation of ρ and θ using the most accurate RDS method to generate numbers for the simulation, where the prior distribution of ρ and θ is the one given in (10) because it is the most general noninformative prior.

Bayesian estimation of ρ and θ using the RDS method
Suppose that our data consists of the first n observations x = (x 1 , x 2 , . . . , x n ).
By using the prior in (10) and the Bayes theorem, we get the following posterior: which is proportional to Knowing that the parameters ρ and θ are independent, the conditional posterior distributions of ρ and θ are respectively and The use of RDS method requires the computation of the inverse cumulative distribution functions F −1 and H −1 , so we put The inversion method leads to the following ρ and θ values of the conditional posterior distributions which will be used for the generation of observations: and where the variables R 1 and R 2 following U[0, 1] are generated using the RDS algorithm given in subsection 4.3 (see Appendix 1 for the generation of the ρ values and Appendix 2 for the generation of the θ values).

Adaptive RDS algorithm to AR(1) model from the Bayesian perspective by using the R language
In this section, we present a new Bayesian Adaptive AR(1) RDS (B2ARDS) algorithm to estimate the parameters of the model given in (2).

Data structure
For each input random variable, we define a record with the following structure. N: an integer defining the number of replications. n: an integer defining the sample size. p: a prime number defining the size of the sub-sample. P: an array of generated primes.  = ρ), n, innov = rexp(n, rate = θ)).
(d) Compute the values x and S, using formula (6) and nx − ρS.
(e) Before each simulation run, generate a vector of prime numbers smaller than the sample size n and store them in a vector P. (f) Randomize the vector P. (g) Compute the sum of primes successively and stop when the sum is n and collect the number of primes m.
(h) Initialization of the sub-run defined by the prime p. (i) For each prime number p j , j = 1, . . . , m, we do the following sub-run.
(f 1 ) Compute the regular numbers r i , i = 1, . . . , p j using the formula (2) and store them in an array r.
(f 2 ) Randomize the r i by the command: > R = sample(r) and store them in an array R.
. . , p j and store them respectively in arrays ρ and θ. (f 4 ) Compute the simulation estimatesρ j andθ j in the j-th sub-run such that

Simulation study
In the following, we present a simulation study using the R language. We generate samples of n observations x 1 , x 2 , . . . , x n from the AR(1) model given in (5) Tables 1  and 2, together with ρ 0 and θ 0 , computed using formulas (7) and (8), respectively (see Appendix 3 for details). At the end of the simulation, we collect the sequences ofρ j andθ j estimates resulting from step (f 4 ) of the proposed algorithm. For N = 10,000 replications, sequences ofρ j estimates for initial values of ρ = 0.1, 0.3, 0.6 are represented in Figure 1 while Figure 2 represents sequences ofθ j estimates for the same ρ values as in Figure 1.

Interpretation of tables
According to the simulation results given in Tables 3 and 4, we have very small variances (given in parentheses). So, we can say that the sequences of RDS values ρ i and θ i , generated by the proposed B2ARDS algorithm are homogeneous. Also, in each caseρ RDS is too close to the true value of ρ while, for some values of β,ρ RDS is closer to the true value than ρ 0 . It can be seen in Figure 1 that when ρ is close to 0, the best values ofρ RDS are obtained when β ∈ [2.5, 16]; while, when ρ is close to 1, the best values ofρ RDS are obtained for high values of β (β ∈ [50, 111]). So, there are some values for ρ and β for which the resulting Bayesian estimate is ideal in accuracy.
Forθ B RDS , the results in Table 4 and Figure 2 also prove that, for each case, the Bayesian estimate of θ is very close to the true value of θ and for some cases (n = 20 and ρ = 0.6), it is a better estimate than θ 0 . We notice here that the best values ofθ B RDS are obtained when ρ is close to 1.

Interpretation of figures
In each figure, the red line represents the true value of the parameter, while the green line represents the value of the resultingρ RDS andθ RDS estimators for the j-th sub-run. We notice that the points which represent the estimateŝ ρ j andθ j are concentrated in a given band for all figures and nearly all points are concentrated around their mean, which represents in a Bayesian sense the resultingρ RDS andθ RDS estimators of the parameters ρ and θ , respectively.  We can see from all figures that the resulting estimate is not only very close to the true value of the parameter in each case, but it coincides with it in some cases (Figure 1 (c) and Figure 2 (c)).

Conclusion
The proposed algorithm B2ARDS has allowed us to generate refined descriptive samples from input random variables for Bayesian estimation of AR(1) parameter in Monte Carlo simulation using the R software. As a consequence and apart from its already proved efficiency in the classical sense, we argue that the RDS method performs well even in a Bayesian context and allows the estimation of the AR(1) parameter using the Bayesian approach.

Disclosure statement
No potential conflict of interest was reported by the author(s).