A New INAR(1) Model for Z-Valued Time Series Using the Relative Binomial Thinning Operator

A new first-order integer-valued autoregressive process (INAR(1)) with extended Poisson innovations is introduced based on a signed version of the thinning operator, called relative binomial thinning operator, which can be considered as an extension of standard binomial thinning operator introduced by Steutel, F.W. and van Harn, K. (1979. Discrete analogues of self-decomposability and stability. Ann. Probab. 7: 893–899). It is appropriate for modeling Z-valued time series and either positive or negative correlations. Some properties of the process are established. Conditional least squares, Yule–Walker and conditional maximum likelihood methods are considered for the parameter estimation of the model. Moreover, simulation experiments are carried out to attest to the performance of the estimation methods. The applicability of the proposed model is investigated through a practical data set of the Saudi stock market.


Introduction
Discrete variable time series data are fairly common in practice, which has attracted the attention of many researchers. The early studies about modeling this type of data have been conducted by McKenzie (1985) and Al-Osh and Alzaid (1987), who suggested the first-order integer-valued autoregressive (INAR(1)) models based on the binomial thinning operator introduced by Steutel and van Harn (1979). Afterward, the researchers have introduced alternative INAR(1) process with different marginal and innovation distributions based on this operator. Moreover, various modification of the binomial thinning operator have been proposed for modeling count time series. A good review on INAR models can be found in Weiß (2008) and Scotto et al. (2015).
All the above models mainly are applicable to analyze a time series with non-negative integer-valued, but in practice, one can encounter also integer-valued time series data that include negative values. For example, in stock market, we analyze intra-daily stock prices the changes belongs on a discrete valued set, which can be represented by Z, known also as ticks (the price can go up or down on certain predefined ranges of value). In image analysis, we model the distribution of intensity differences between two neighborhood pixels. Thus, since intensities are discrete actually, we have the difference of two discrete variables. In football, the results of a match can be expressed in Z as the difference in the number of goals, a value often used in sport betting. In medicine, for some diseases we care about the difference in some count variables, such as number of pimples, before and after some treatment is applied. Moreover, in some fields, e.g. meteorology, by sake of simplicity and to facilitate the reading and interpret, data recorded are rounded, which lead to data with Z as support. Furthermore, when we analyze a non-stationary integer-valued time series with non-negative values, we use the differencing operator to achieve stationary. Therefore, we may obtain a time series on Z.
To the best of our knowledge, limited studies have been conducted on modeling time series data defined on the set Z, with an autoregressive structure. Kim and Park (2008) have proposed an integer-valued autoregressive process of order p ≥ 1 with signed binomial thinning operator. Kachour and Yao (2009) have introduced the first-order rounded integer-valued autoregressive process, based on the rounding operator. A more general setup has been introduced since then by Kachour (2014). Recently, Liu et al. (2021) have proposed a semiparametric autoregressive model with a log-concave innovation, based on the model proposed by Kachour (2014). Zhang et al. (2010) have presented an integer-valued autoregressive processes of order p ≥ 1 with signed generalized power series thinning operator. Kachour and Truquet (2011) have proposed an extension of INAR models using a modified version of the generalized thinning operator, which is different from that has been introduced by Kim and Park (2008). Zhang et al. (2012) have introduced a random coefficient process called generalized random coefficient first-order integer-valued autoregressive process with signed thinning operator. In literature, there exists a family of models defined on Z that arise as the difference between two discrete distributions. For example, Freeland (2010) has defined the true integer-valued autoregressive process of order one as the difference of two Poisson INAR(1) processes, Barreto-Souza and Bourguignon (2015) have proposed the skew INAR(1) process, which is defined as the difference of two INAR(1) process with geometric marginal distribution, Bourguignon and Vasconcellos (2016) have defined the new skew integer-valued process as the difference between a Poisson INAR(1) process and a geometric INAR(1) process, and Taveira da Cunha et al. (2018) have proposed a new integer-valued autoregressive process with generalized Poisson difference marginal distributions based on difference of two quasi-binomial thinning operators. It is also important to mention the existence of certain studies that deal with the multivariate Z-valued autoregressive models. For example, Bulla et al. (2011) have introduced the bivariate autoregressive integer-valued time-series models, based on the signed thinning operator. More recently, Chen et al. (2023) have proposed a new bivariate Z-valued autoregressive model, based on the bivariate Skellam distribution.
In this paper we introduce a new integer-valued process, denoted by EP-RBINAR(1), by considering a parametric assumption on the common distribution of the counting sequence of the signed integer-valued autoregressive process proposed by Kachour and Truquet (2011). Thus, the EP-RBINAR(1) can fit integervalued times series with possible negative values. Moreover, similar to an AR(1) process, the autocorrelation function of EP-RBINAR(1) can also have negative values. Indeed, unlike models that arise as the difference between two discrete distributions, EP-RBINAR(1) process can fit Z-valued time series without take into account this specific construction. Moreover, in comparison with the model proposed by Kachour and Yao (2009), the EP-RBINAR(1) doesn't have the issue related to small lack of identifiability on the parameter. Note that, the thinning part of the new process can be considered as the sum of a "discrete random walk", where one can go a one step forward, or one step back, or keep still. Furthermore, innovations of EP-RBINAR(1) process follow the extended Poisson distribution introduced by Bakouch et al. (2016). This distribution is defined on Z, having a dispersion flexibility, and under some assumptions, it can be approximated by the Gaussian distribution.
The paper is structured as follows. EP-RBINAR(1) process is formally defined in Section 2 and some of its properties are outlined. In Section 3, estimation methods for the process parameters are proposed. Section 4 discusses some simulation results for the estimation methods. Moreover, the EP-RBINAR(1) model is applied to a practical data set of the Saudi stock market. Finally, the proofs of all propositions and theorems are contained in the appendix.
New Model for ℤ-Valued Time Series 2 The EP-RBINAR(1) Process

Definition
In this section, we first consider a special version of the signed thinning operator, originally proposed by Latour and Truquet (2008). It is referred to "relative binomial thinning operator" and is an extension of the classical Steutel and van Harn operator to Z-valued random variables. Then, we propose a new signed integer-valued autoregressive process based on this operator with extended Poisson innovations.
Definition 1. (Relative binomial thinning operator) Let {Y i } i∈N be a sequence of independent identically distributed (i.i.d) integer-valued random variables with the common distribution F, and independent of an integer-valued random variable X. The relative binomial thinning operator, denoted by F∘, is defined as where for any integer x ≠ 0, sign(x) = 1 if x > 0 and −1 if x < 0, and F is defined as follow with 0 ≤ α ≤ 1 (F is called a relative Bernoulli distribution).
Remark 1. For classical Steutel and van Harn operator, X is positive and {Y i } i∈N is a sequence of Bernoulli variables (i.e. Y i Ω ( ) = {0, 1}), where "Y i = 1" represents "success" (survival) and "Y i = 0" represents "failure" (death). However, for the relative binomial thinning operator, Y i can be seen as a "discrete random walk", where "Y i = 1" can be interpreted as a one step forward, "Y i = −1" as one step back, and "Y i = 0" as keep still.
Remark 2. Note that originally Chesneau and Kachour (2012) introduced the relative binomial distribution, without giving it this name. Moreover, one can see that . Moreover, let l be a positive integer, such that l > 2. Thus, one can deduce that In the following lemma, some useful basic properties of the relative binomial thinning operator are presented.
Lemma 1. Let x ∈ Z \ {0}, and set F∘X|(X = x) = F∘x, then the relative binomial thinning operator in Definition 1 has the following properties: (i) ∀y ∈ Z, where I A (y) denotes the indicator function, which has the value 1 when y takes values in A and has the value 0, otherwise and Using the relative binomial thinning operator in Definition 1, we now introduce a new process to be used for modeling Z-valued time series.
Definition 2. A sequence {X t } t∈Z is said to be a EP-RBINAR(1) (first-order Extended-Poisson Relative Binomial Intergred-valued Autoregressive) process if it has the following representation: where F is defined in (2) and {ϵ t } t∈Z is a sequence of i.i.d random variables, called innovations, following an extended Poisson distribution introduced by Bakouch et al. (2016). It is denoted by ϵ t ∼ E-Po(p, λ) and defined by the following probability mass function (pmf): where λ > 0 and 0 ≤ p ≤ 1. All {Y i } i∈N in counting sequence are independent of ϵ t . Moreover, ϵ t and the sigma-algebra σ(X t−1 , X t−2 , …) are supposed to be mutually independent.
New Model for ℤ-Valued Time Series Remark 3. Given that innovations (i.e. {ϵ t } t∈Z ) have a known distribution, our model defined in (3), can be considered as a special case of P-SINAR(1) process introduced by Chesneau and Kachour (2012), which is a particular case of SINAR(1) process presented by Kachour and Truquet (2011).
Remark 4. The use of E-Po(p, λ) distribution guarantees innovations with possible negative values and having more dispersion flexibility (i.e. the dispersion of E-Po(p, λ) distribution depends on the value of p).
Remark 5. From Bakouch et al. (2016), the corresponding formulas for mean, variance and probability generating function (pgf) of E-Po(p, λ) distribution with the pmf (4) are: Obviously, if ϵ t ∼ E-Po(p, λ), then |ϵ t |∼ Po(λ). As a result, we have Moreover, for any positive integer l > 2, we have Remark 6. Under the previous parametric assumptions, the EP-RBINAR(1) process can be seen as an extension on Z of the Poisson INAR(1) process, originally proposed Al-Osh and Alzaid (1987).
Remark 7. Based on its construction, one can deduce that, under specific parameters values, the EP-RBINAR(1) process can provide zero-inflated integer (positive and negative) observations. For example, suppose that (by sake of simplicity) that the innovation distribution is symmetric (i.e. p = 0.5) and concentrated on zero (i.e. λ < 1), and α is chosen so that zero be the mode of distribution F. Thus, in this case, it is very likely that zero is the most represented value of the process. Figure 1 shows plot and distribution of n = 1000 observations simulated from EP-RBINAR(1) process, where actual values are α = 0.45, p = 0.5, and λ = 0.3. One can see that the empirical frequency associated with zero equals 60 % and there is an almost balance between the distribution of positive and negative values.

Properties
The stationarity of the proposed EP-RBINAR(1) process is established in the following theorem.
Theorem 1. Suppose that α ∈ 0, 1 ] [\ { 1 2 } and 0 < p < 1. Hence, one can deduce that the process X t defined in (3) has a unique stationary solution. Moreover, for any positive integer l > 2, we have In the sequel, under stationary conditions and using Lemma 1, we derive some properties of the EP-RBINAR(1) process.
. Under stationary conditions presented in Theorem 1, one can deduce that The probability generating function of the stationary EP-RBINAR(1) process is obtained as and Φ ϵt (s) is the pgf of {ϵ t } given in (5).
Remark 8. The stationary conditions associated with the EP-RBINAR(1) process are similar to that of a classical AR(1) process. Indeed, one can see that where ξ t is a stationary process defined by Moreover, one can see, based on properties of signed thinning operator that ξ = (ξ t ) is an uncorrelated process. Thus, the EP-RBINAR(1) process has the same autocorrelation structure as a classic AR(1) process. In other words, the autocorrelation function of the EP-RBINAR(1) process is obtained as Therefore, the spectral density function of the model, for ω ∈ −π, +π ( ); is obtained as Under stationary conditions, the marginal probability function of the EP-RBINAR(1) process is given in the following proposition.
Proposition 2. Let {X t } be an EP-RBINAR(1) process, with 0 < α < 1, α ≠ 1 2 , and 0 < p < 1. Thus, for all a ∈ Z, we have where The EP-RBINAR(1) process's one-step transition probability, denoted by and for i ≠ 0 where P(ϵ t = j) is the pmf of {ϵ t } defined by (4). Moreover, using the first-order dependence of the process, the joint probability can be obtained as The k step-ahead conditional mean (usually used for time series forecast issues) of the process is derived in the following proposition.

Parameters Estimation
Let (X 1 , …, X m ) t be a vector of observations from our model defined in (3), and denote the unknown parameter vector. In order to estimate θ, we propose three estimation methods, namely, Yule-Walker (YW), conditional least square (CLS), and conditional maximum likelihood (CML).

Yule-Walker Estimation
In order to obtain estimations in this method, we solve a set of equations that are resulted from equating the theoretical and empirical aspects at the same time. Thus, the YW estimation of α is obtained bŷ and the YW estimation of p and λ are obtained by solving the following equationŝ Remark 10. A special case of our model, when we suppose that α = p. In this case, the YW estimation of α (and p) isα YW =ρ (1) + 1 2 , and the YW estimation of λ is given by
Remark 12. A special case of the proposed model, when we suppose that α = p. In this case, the CLS estimation of α (and p) is given aŝ and the CLS estimation of λ is given bŷ

Conditional Maximum Likelihood Estimation
From the joint probability function (9), the conditional log-likelihood function for the EP-RBINAR(1) model can be written as where P(ϵ t = j) is the pmf of {ϵ t } defined by (4) and I A (X s ) denotes the indicator function, which has the value 1 when X s takes values in A and has the value 0, otherwise.
The conditional maximum likelihood (CML) estimatorθ CML = (α CML ,p CML ,λ CML ) ′ of θ = (α, p, λ)′ is defined as the value of θ that maximize the conditional log-likelihood function in (14). Since equating the first-order conditional log-likelihood derivatives to zero leads us to a complicated system of equations, the CML estimates are achieved using numerical methods.

Empirical Study
This section includes two subsections. In the first part, the performance of estimation methods that we used for parameters of the proposed process is compared through a simulation study. To ensure the applicability of the proposed model, the second part is devoted to analyzing a practical data set of the Saudi stock market.

The Case When α = p
In this section, we aim to test the efficiency of the parameter's estimation discussed in the previous section. For the sake of simplicity, we suppose that α = p. Thus, we simulate (using the R programming language) 1000 paths of length 100, 500, 1000 and 10,000. These paths are simulated using Equation (3) with three sets of parameters: (a) α = p = 0.1 and λ = 1; (b) α = p = 0.6 and λ = 2; (c) α = p = 0.9 and λ = 3. Mean values of YW, CLS, and CML estimates for each set of parameters are given in Table 1. The standards deviations of the estimates are stated in brackets under the estimated values.
Remark 13. In the case when α = p, one can use results presented in Remarks 10 and 12 to calculate YW and CLS estimates. However, CML estimates are calculate by using numerical methods. Explicitly, we use the nlm function (Non-Linear Minimization, from package "stats", software R), to find values that maximize the conditional log-likelihood function, denoted by (14) (where we consider that α = p in this equation).
In general way, from the obtained results, one can deduce that YW, CLS and CML methods provide estimates that are quite close to the actual values. Moreover, the performances of YW and CLS methods are very similar. For these methods, one can New Model for ℤ-Valued Time Series see that when the length of series is over 500, they converge quickly to the actual parameter values. However, results of CML method are more efficient even for a small length of series. For all the proposed methods, we can notice that the standard deviation of the estimates decreases as the length series increases. Furthermore, in most of studied cases, the CLM method provides the lower standard deviation of estimates. However, for YW and CLS methods, it is important to point the convergence rate is less quickly for λ. Indeed, this can be explained by the fact that, for both methods, the estimate of λ depends on the estimated value of α. Thus, for these methods, we find that a large current value of λ implies the standard deviation of the estimates is high, especially when the size of the series is small. This remark is illustrated thanks to the simulation results obtained for parameters set (c) (for details, see Table 1). The boxplots for the parameter set (c) are given in Figures 2-4   methods, respectively. One can see that, contrary to the other methods, the estimates from the CML method have few or no outliers.

for CLS, YW and CML
On the other hand, one can see, based on the simulation results of YW and CLS methods, obtained for parameter set (b) when n = 100, that the standard deviation of λ estimates is very high for both methods (see Table 1). In fact, these results are not surprising cause that α (= p) actual value is close to 0.5 (which is an excluded value of the parameter). So, with a small size of series the estimate value of α may be close to 0.5, which can "explode" the estimate value of λ. However, this issue is not present when we use the CML method (in this case, even with a small length series, the standard deviation is weak, see Table 1). Finally, fitting to the Gaussian distribution is illustrated in Figures 5 and 6 for YW, CLS and CML estimators (with the parameter Figure 5: Normal Q-Q plots for errors (α YW − 0.9), (α CLS − 0.9) and (α CML − 0.9), when n = 10,000. The horizontal axis indicates theoretical quantiles and the vertical axis indicates sample quantiles.

The Case When α ≠ p
As we mentioned in the previous section, the performance of the YW and CLS methods are quite close. Moreover, we have also seen that if α ≠ p, then the estimates of these methods do not have an explicit form. Thus, in this section, we will compare the performances resulting from the YW and CML methods based on one set of parameters: α = 0.75, p = 0.4, and λ = 2. We simulate (using the R programming language) 10,000 paths of length 100, 250, and 1000. These paths are simulated using Equation (3) with the above set of parameters. To calculate the CML estimates we use the nlm function (Non-Linear Minimization, from package "stats", software R) to find values that maximize the conditional log-likelihood function, denoted by (14). On the other hand, the YW estimates of α is given by (11). Once we estimate the value of parameter α (based on (11)) and substitute this value into Equation (12), we get the system of two equations with two unknown variables λ and p. To find YW estimates of these parameters, this system has to be solved numerically. Thus, the numerical procedure will be conducted in programming language R, where we will use the BB package. Mean values of YW and CML estimates for each parameters are given in Table 2. The standards deviations of the estimates are stated in brackets under the estimated values. Once again, we can see that the precision of these estimates (from both methods) increases when the size n increases. Remark 14. In order to study the compatibility of the proposed model and the estimation results, we generated, based on Equation (1), n = 1000 observations of the EP-RBINAR(1) process, with α = 0.75, p = 0.4, and λ = 2. Basic descriptive statistics associated with the simulated data are represented in the following Using the simulated data and the procedure presented in the above section, we compute the CML estimates of the process parameters:α = 0.7431100,p = 0.3854399 andλ = 1.9496662. Now, using these estimates and based on the properties of the EP-RBINAR(1) process (see Proposition 1 and Remark 8), we calculate the indicators presented in the above New Model for ℤ-Valued Time Series

A Practical Data Example
Here, we present an application of the model EP-RBINAR(1) based on real-world data from Saudi stock market. In 2007, the minimum amount of change was 0.25 SR (Saudi Riyal) for all stocks. The daily close price as number of ticks (ticks = close price × 4) in 2007 for Saudi Telecommunication Company (STC) stock is considered. Note that these data were originally introduced by Alzaid and Omair (2014) as an application of a new integer-valued model with a Poisson difference marginal distribution. This model can be used as a tool to model non-stationary count data. Basic descriptive statistics concerning these data are presented in Table 3. After, examination of the time series plot of the data, authors (Alzaid and Omair (2014)) notice a non-stationarity in the mean (this issue was confirmed based on a sustained large autocorrelation function (ACF) and exceptionally large first lag partial autocorrelation function (PACF)). Thus, authors considered that differencing is needed. Explicitly, if {X t } is the process behind the above data, then authors propose to study the lag one differenced process, denoted by {Y t }, where Y t = X t − X t−1 . Indeed, one can consider that data observed from process {Y t } as the difference between two consecutive daily ticks (associated with the close price). Thus, a data with negative values represents the number of ticks "lost" at the close between two consecutive days (i.e. a drop in price), while a positive value reveals the number of ticks "gained" at the close between two consecutive days (i.e. a price increase). Finally, if the observed data is zero, this means that the number of ticks at closing has remained "stable" between two consecutive days (i.e. the price has not changed).
The time series plot of the differenced stock is illustrated in Figure 7 and basic descriptive statistics concerning the differenced data are presented in Table 4. Thus, Figure 7 shows the stationarity in the mean is now verified. The ACF and PACF for differenced stock are shown in Figure 8. These figures show that the lag one correlation is positive and significant, which implies that the differenced process have the same autocorrelation structure of an AR(1) process. Moreover, after the lag one differenced stage, the data is integer with negative values.
Remark 15. Bourguignon and Vasconcellos (2016) have used data, dating from 2012, providing also from Saudi Telekom. The choice of these data is consistent with the constructed as being the difference between a Poisson INAR(1) process and a geometric INAR(1) process). At the first glance, it is quite possible to consider the model proposed by Bourguignon and Vasconcellos (2016) to fit the differenced data of our application. However, we cannot say that the differenced data used in our   application comes from the difference between two sets of independent positive integer-valued data and therefore it goes against the philosophy of the model proposed by Bourguignon and Vasconcellos (2016). However, since our model belongs to SINAR(1) process family, it can be considered to fit data used by Bourguignon and Vasconcellos (2016) without take into account how data have been constructed.
Remark 16. Based on properties of the process and using the estimates of the unknown parameters, one can deduce that the mean is equal to 0.3311299 (which is not close enough to 0.004048583 (mean of the data) however, one can consider that both values are near to zero), the standard deviation equals 4.246681 (which is close to 4.715049, the standard deviation calculated from the data), and the lag-1 autocorrelation is equal to 0.289729 (which is close to 0.2194211, lag-1 autocorrelation calculated from the data). These first results support our choice of the EP-RBINAR(1) model to fit the data.
Moreover, based on the conditional one-step ahead mean, the estimated residuals are computed asê Remark 17. In general,ŷ t are real-valued, a mapping into the discrete support of the series is obtained by rounding to the nearest integer.
To verify the adequacy of the considered model, the ACF and the PACF of the estimated residuals are plotted in Figure 9. Clearly, this graph shows that the residuals can be considered as observations from a white noise. Moreover, in Figure 10, we make a comparison between the empirical frequencies ofê t and their associated theoretical probabilities resulting from E-Po(0.5345621, 3.402455). This comparison shows that the frequencies and probabilities associated with the observed modalities are relatively close. This is one of reasons to justify our choice to use the EP-RBINAR(1) model to fit the data.
Remark 18. Based on the conditional one-step ahead standard deviation, the estimated standardized Pearson residuals are computed aŝ For details concerning the standardized Pearson residuals and its role to check model adequacy for count time series, see Weiß et al. (2019). Note that the empirical mean, median, and standard deviation, of the estimated standardized Pearson residuals, are, respectively, equals: −0.05694, 0, and 1.127245. Thus, one can see that the empirical mean is very close to zero and empirical standard deviation is close to one.  Moreover, the ACF, the PACF, and the fitting to the Gaussian distribution (Normal Q-Q plot) of the estimated standardized Pearson residuals are plotted in Figure 11. These results can be considered as additional arguments whose objective is to check the model adequacy to fit the data.
Remark 19. As mentioned before, the data processed in this section have been studied by Alzaid and Omair (2014). Indeed, authors proposed their model, called PDINAR(1) process, to fit the data. However, this process has two variants, denoted by PDINAR + (1) and PDINAR − (1), depending on the sign of the correlation. In practice, to fix which variant should be used to fit the data, it is necessary to have a preliminary idea of the sign ofρ(1). Thus, since the empirical first-order autocorrelation is positive, Alzaid and Omair (2014) considered the PDINAR + (1) to study the data. However, for the use of our model, we don't need to know the sign ofρ(1). On the other hand, Alzaid and Omair (2014) show the that one-step ahead least squares predictions can be calculated bỹ y t = 0.2201y t−1 + 10.4088 − 10.3953.
These predicted values need also to by rounding to the nearest integer to be adequate with the discrete support of the data. From Figure 12, one can deduce that the predictions (based on conditional one-step ahead mean) results provided from PDINAR + (1) process and those from the EP-RBINAR(1) process, are very similar. This implies that both models have a prediction performances very close. Indeed, the Root Mean Squared Error (RMSE) equals 4.604699 (resp. 4.590552), for forecasting values (based on conditional mean) provided from PDINAR + (1) (resp. EP-RBINAR(1)) Figure 11: From top to bottom: ACF, PACF, and normal Q-Q plot of estimated standardized Pearson residuals of the differenced STC stock data set.
process. However, one can see also that for both models there is a rather discrepancy between the forecast and the actual data.
Remark 20. As mentioned above, conditional one-step ahead mean leads to realvalued forecasts, which is incoherent with the support of the integer-valued time series. Recently, coherent techniques have been introduced in order to provide integer-valued forecasts. One of the most used technique is the conditional median (for details, see Homburg et al. (2019)), which is obtained by calculating the conditional probability of each possible integer value, then selecting a forecast value with a cumulative conditional probability greater than 0.5. From Figure 13, one can see that discrepancy between conditional median forecasts and actual data have been

Concluding Remarks
In this paper, we introduce a new stationary first-order integer-valued autoregressive process, denoted by EP-RBINAR(1), where the thinning part of the process can be considered as the sum of a random walk. The new model has several advantages: possible negative values for time series, possible negative values for autocorrelation, and an innovation structure that can be interpreted as an extension on Z of the Poisson distribution. Moreover, under specific parameters values, EP-RBINAR(1) can provide zero-inflated integer (positive and negative) observations. The main properties of the model are derived. Then we considered the problem of parameter estimation and derived properties of the Yule-Walker, conditional least square, and conditional maximum likelihood estimators. An example of application to a real-world data set, based on number of ticks (associated with close price) for Saudi Telecommunication Company (STC) stock (in 2007), illustrates the importance and potentiality of the new model. We are convinced that this new process may attract many other wider applications in time series analysis. As part of future research, it would be of interest to extend the proposed process, studying the process of order p.