NUMERICAL SOLUTION TO AN INTEGRAL EQUATION FOR THE KTH MOMENT FUNCTION OF A GEOMETRIC PROCESS

In this paper, an integral equation for the kth moment function of a geometric process is derived as a generalization of the lower-order moments of the process. We propose a general solution to solve this integral equation by using the numerical method, namely trapezoidal integration rule. The general solution is reduced to the numerical solution of the integral equations which will be given for the third and fourth moment functions to compute the skewness and kurtosis of a geometric process. To illustrate the numerical method, we assume gamma, Weibull and lognormal distributions for the rst interarrival time of the geometric process. 1. INTRODUCTION The geometric process (GP), which is a natural generalization of a renewal process (RP), is an important stochastic monotone model used in many areas of statistics and applied probability, especially for statistical analysis of series of events. Since its introduction by [9], many researcher and authors made a signicance effort on GP by publishing more than 200 papers. For instance, the GP has been used as a model in modelling of an epidemic disease [8], software reliability [16,18], maintenance [23], warranty analysis [5,12] and electricity prices [7]. This process is dened as follows. Denition 1. Let fN(t); t 0g be a counting process (CP) and Xi be the interarrival time between (i 1)th and ith event of this process for i = 1; 2; : : : . The CP fN(t); t 0g is called a GP with the ratio parameter a if there exists a real number a > 0 such that a Xi; i = 1; 2; : : : are independent and identically distributed 2020 Mathematics Subject Classication. Primary 60K99, 65R20; Secondary 60E10 , 62E10. Keywords and phrases. Geometric process, moment functions, skewness and kurtosis, numerical solution. mpekalp@ankara.edu.tr-Corresponding author; aysenurydgd@gmail.com 0000-0002-5183-8394; 0000-0001-6439-8431. c 2021 Ankara University Communications Facu lty of Sciences University of Ankara-Series A1 Mathematics and Statistics 731 732 M. H. PEKALP, A. AYDO 1⁄4 GDU (iid) random variables with a distribution function F , where F is the distribution function of the rst interarrival time X1. The GP is also called quasi-renewal process with ratio parameter = 1=a by [18]. Let fN(t); t 0g be a GP with the ratio parameter a and Fi be the distribution function of Xi for i = 1; 2; : : : . Then, it is easy to verify that Fi (x) = F a x for i = 1; 2; : : : . Further, it can be shown that GP is stochastically increasing if a < 1 and stochastically decreasing if a > 1. When a = 1, GP reduces to a RP. By considering the distribution functions of the interarrival times, one of the important di¤erences between RP and GP can be given as follows. In RP, the distribution function of the interarrival times remains same over is, that is Fi (x) = F (x) for i = 1; 2; : : : . However, in GP, the distribution of the interarrival time Xi does not remain same over is, that is Fi (x) = F (a x) for i = 1; 2; : : : . This provides some advantages to the GP in applications, in particular for reliability mathematics since it can be used as a model for deteriorating systems which may have decreasing working times between failures. In order to understand the place of the GP model in the literature, see the following papers, [1, 19,20,21,22]. Let fN (t) ; t 0g be a GP with ratio parameter a. Set S0 = 0 and Sn = Pn i=1Xi for n = 1; 2; : : : . Thus, Sn is called nth arrival time of the GP. The distribution function of Sn is given by F1 F2 Fn (t), where * denotes the Stieltjes convolution and Fi (t) for i = 1; 2; : : : is the distribution function of Xi. Further, since the events (Sn t) and (N (t) n) are equivalent, the probability distribution of the random variable N (t) is given by P (N (t) = n) = F1 Fn (t) F1 Fn+1 (t) for each xed t 0. Now, let give the following theorem which states the existence of the moments of the GP. The proof of this theorem can be found in many manuscripts and monographs, for example, [6, 11]. Theorem 2. Consider the GP fN (t) ; t 0g with ratio parameter a. If a 1, then Mk (t) = E N (t) < 1 for all t 0 and k 0. If a > 1 and F (t) > 0 for all t > 0, then the moments of N (t) are innite. Let fN (t) ; t 0g be a GP with ratio parameter a. The mean value function of a GP, which is also called the geometric function, is given by M1 (t) = E (N (t)). M1 (t) is the number of the events that have occurred by time t. The geometric function M1 (t) satises the following integral equation. M1 (t) = F (t) + Z t 0 M1(a (t x))dF (x) ; t 0: (1) The second moment function of a GP is given by M2 (t) = E N (t) . [15] show that M2 (t) satises the following integral equation. M2 (t) = 2M1 (t) F (t) + Z t 0 M2 (a (t x))dF (x) ; t 0: (2) NUMERICAL SOLUTION TO AN INTEGRAL EQUATION 733 According to the Theorem 1, for a 1, the geometric function M1 (t) and the second moment function M2 (t) are nite for all t 0. Furthermore, if F is continuous, then the integral equations (1) and (2) can be solved uniquely although M1 (t) and M2 (t) cannot be obtained in analytical forms. In the case of a > 1, M1 (t) and M2 (t) are innite for all t > 0. When the GP model is used as a model for the data sets come from series of event, the distribution function of the rst interarrival time is assumed to be one of four common lifetime distributions as exponential, gamma, Weibull and lognormal. See [13] for the details how to discriminate the lifetime distributions in GP model. Under a lifetime distribution assumption, it is of importance to calculate the moment functions of the GP. Many researchers and authors made some studies on the rst and second moment functions of the GP by considering these lifetime distributions. [6] deal with the boundary problem forM1 (t). [17] propose a numerical method, namely trapezoidal integration rule, for M1 (t) with the help of the integral equation (1). In addition, [4] and [3] suggest power series expansions for M1 (t) depending on the integral equation (1). [2] obtain the numerical calculation and Monte Carlo estimation of the variance function, which is M2 (t) M 1 (t) ; by using the convolution of the distribution functions. Alternative methods for computing M2 (t) are given in [15]. They adapt the Tang and Lams method to M2 (t) and propose a power series expansion forM2 (t) with the help of the integral equation (2). Further, [14] show the asymptotic solution of the integral equation for the second moment function. To the best of our knowledge, in the literature, there is no study about the higher-order moment functions of the GP. However, in order to calculate, for instance, the skewness and kurtosis of the process, the third and fourth moment functions are required. Moreover, to compare the estimators proposed for some parameters of GP and to examine some theoretical properties of the process, higher-order moment functions should be known. The rest of the paper organized as follows. In section 2, rstly, we obtain the integral equations for the third and fourth moment functions of the GP. Then, a generalization of the integral equation for kth moment function of the GP is given. We adapt the Tang and Lams numerical method for the kth moment function of the GP with the help of the integral equation given for this function. Then, we reduce this general approximation to third and fourth moment functions of the GP in Section 3. Illustrative examples are provided in Section 4. Conclusion remarks are presented in Section 5. 2. INTEGRAL EQUATIONS FOR THE MOMENT FUNCTIONS OF THE GP In this section, rstly, integral equations for the third and fourth moment functions of the GP are obtained. Then, in general, we present an integral equation for the kth moment function of the GP. 734 M. H. PEKALP, A. AYDO 1⁄4 GDU Let fN (t) ; t 0g be a GP with ratio parameter a and let us assume that the rst interarrival time X1 follows the distribution function F . The third moment function of the GP is given by M3 (t) = E N (t) ; t 0: Conditioning on the rst interarrival time X1, we have M3 (t) = E N (t) = Z 1 0 E N (t) j X1 = x dF (x): Since E(N (t) j X1 = x) = E(1 +N (a (t x))); x < t and E(N (t) j X1 = x) = 0; x t, we rewrite the equality as follows.


INTRODUCTION
The geometric process (GP), which is a natural generalization of a renewal process (RP), is an important stochastic monotone model used in many areas of statistics and applied probability, especially for statistical analysis of series of events. Since its introduction by [9], many researcher and authors made a signi…cance effort on GP by publishing more than 200 papers. For instance, the GP has been used as a model in modelling of an epidemic disease [8], software reliability [16,18], maintenance [23], warranty analysis [5,12] and electricity prices [7]. This process is de…ned as follows. De…nition 1. Let fN (t); t 0g be a counting process (CP) and X i be the interarrival time between (i 1)th and ith event of this process for i = 1; 2; : : : . The CP fN (t); t 0g is called a GP with the ratio parameter a if there exists a real number a > 0 such that a i 1 X i ; i = 1; 2; : : : are independent and identically distributed The GP is also called quasi-renewal process with ratio parameter = 1=a by [18]. Let fN (t); t 0g be a GP with the ratio parameter a and F i be the distribution function of X i for i = 1; 2; : : : . Then, it is easy to verify that F i (x) = F a i 1 x for i = 1; 2; : : : . Further, it can be shown that GP is stochastically increasing if a < 1 and stochastically decreasing if a > 1. When a = 1, GP reduces to a RP.
By considering the distribution functions of the interarrival times, one of the important di¤erences between RP and GP can be given as follows. In RP, the distribution function of the interarrival times remains same over i's, that is F i (x) = F (x) for i = 1; 2; : : : . However, in GP, the distribution of the interarrival time X i does not remain same over i's, that is F i (x) = F (a i 1 x) for i = 1; 2; : : : . This provides some advantages to the GP in applications, in particular for reliability mathematics since it can be used as a model for deteriorating systems which may have decreasing working times between failures. In order to understand the place of the GP model in the literature, see the following papers, [1,19,20,21,22].
Let fN (t) ; t 0g be a GP with ratio parameter a. Set S 0 = 0 and S n = P n i=1 X i for n = 1; 2; : : : . Thus, S n is called nth arrival time of the GP. The distribution function of S n is given by where * denotes the Stieltjes convolution and F i (t) for i = 1; 2; : : : is the distribution function of X i . Further, since the events (S n t) and (N (t) n) are equivalent, the probability distribution of the random variable N (t) is given by for each …xed t 0. Now, let give the following theorem which states the existence of the moments of the GP. The proof of this theorem can be found in many manuscripts and monographs, for example, [6,11].
Theorem 2. Consider the GP fN (t) ; t 0g with ratio parameter a. If a 1, then M k (t) = E N k (t) < 1 for all t 0 and k 0. If a > 1 and F (t) > 0 for all t > 0, then the moments of N (t) are in…nite.
Let fN (t) ; t 0g be a GP with ratio parameter a. The mean value function of a GP, which is also called the geometric function, is given by M 1 (t) = E (N (t)). M 1 (t) is the number of the events that have occurred by time t. The geometric function M 1 (t) satis…es the following integral equation.
The second moment function of a GP is given by M 2 (t) = E N 2 (t) . [15] show that M 2 (t) satis…es the following integral equation.
According to the Theorem 1, for a 1, the geometric function M 1 (t) and the second moment function M 2 (t) are …nite for all t 0. Furthermore, if F is continuous, then the integral equations (1) and (2) can be solved uniquely although M 1 (t) and M 2 (t) cannot be obtained in analytical forms. In the case of a > 1, M 1 (t) and M 2 (t) are in…nite for all t > 0.
When the GP model is used as a model for the data sets come from series of event, the distribution function of the …rst interarrival time is assumed to be one of four common lifetime distributions as exponential, gamma, Weibull and lognormal. See [13] for the details how to discriminate the lifetime distributions in GP model. Under a lifetime distribution assumption, it is of importance to calculate the moment functions of the GP. Many researchers and authors made some studies on the …rst and second moment functions of the GP by considering these lifetime distributions. [6] deal with the boundary problem for M 1 (t). [17] propose a numerical method, namely trapezoidal integration rule, for M 1 (t) with the help of the integral equation (1). In addition, [4] and [3] suggest power series expansions for M 1 (t) depending on the integral equation (1). [2] obtain the numerical calculation and Monte Carlo estimation of the variance function, which is M 2 (t) M 2 1 (t) ; by using the convolution of the distribution functions. Alternative methods for computing M 2 (t) are given in [15]. They adapt the Tang and Lam's method to M 2 (t) and propose a power series expansion for M 2 (t) with the help of the integral equation (2). Further, [14] show the asymptotic solution of the integral equation for the second moment function. To the best of our knowledge, in the literature, there is no study about the higher-order moment functions of the GP. However, in order to calculate, for instance, the skewness and kurtosis of the process, the third and fourth moment functions are required. Moreover, to compare the estimators proposed for some parameters of GP and to examine some theoretical properties of the process, higher-order moment functions should be known.
The rest of the paper organized as follows. In section 2, …rstly, we obtain the integral equations for the third and fourth moment functions of the GP. Then, a generalization of the integral equation for kth moment function of the GP is given. We adapt the Tang and Lam's numerical method for the kth moment function of the GP with the help of the integral equation given for this function. Then, we reduce this general approximation to third and fourth moment functions of the GP in Section 3. Illustrative examples are provided in Section 4. Conclusion remarks are presented in Section 5.

INTEGRAL EQUATIONS FOR THE MOMENT FUNCTIONS OF THE GP
In this section, …rstly, integral equations for the third and fourth moment functions of the GP are obtained. Then, in general, we present an integral equation for the kth moment function of the GP.
Let fN (t) ; t 0g be a GP with ratio parameter a and let us assume that the …rst interarrival time X 1 follows the distribution function F . The third moment function of the GP is given by M 3 (t) = E N 3 (t) ; t 0: Conditioning on the …rst interarrival time X 1 , we have we rewrite the equality as follows.
Using the integral equations given in (1) and (2), we have Simplifying the expression, we obtain The fourth moment function of the GP is de…ned by M 4 (t) = E N 4 (t) ; t 0: Using similar arguments in the derivation of the integral equation for M 3 (t), we have Following theorem states the kth moment function of a GP model.
Theorem 3. For any k 2 N, the kth moment function M k (t) of the GP is given by where M 0 (t) = F (t).
Proof. The proof of (5) is given by using the mathematical induction. It is obvious that (5) holds for k = 1 when we consider equation (1). Now, let us assume that (5) holds for any integer k and show that it also holds for integer k + 1. The (k + 1)th moment function of the GP is de…ned by M k+1 (t) = E N k+1 (t) ; t 0: Conditioning on the …rst interarrival time X 1 , we have The following equation can be written by separating the (k + 1)th term.
Since we assume that (5) holds for any integer k, we have When we rearrange the terms, we obtain Hence, the proof is completed.
According to the Theorem 1, for a 1, the function M k (t) is …nite for all t 0. Furthermore, if F is continuous, then the integral equation (5) can be solved uniquely although this function does not have an analytical form. We discuss this problem in the next section. In the case of a > 1, M k (t) is in…nite for all t > 0.

NUMERICAL SOLUTION
In this section, we give a method based on the trapezoidal integration rule for the numerical solution of the integral equation (5). This solution is obtained by recursive calculations with respect to k. Now, let us remind the trapezoidal integration rule as follows.
According to the trapezoidal integration rule, an integral R b a g (t) dt can be calculated numerically as where ft 0 ; t 1 ; : : : ; t n g is a partition of the interval [a; b] such that a = t 0 < t 1 < < t n = b, t i = a + ih, i = 0; 1; : : : ; n and h = b a n . The approximation gives more precise results as the number of the partition increases.
Let fN (t); t 0g be a GP with ratio parameter a < 1. Assume that the …rst interarrival time X 1 has probability density function f . Then, the integral equation (5) can be written as If we substitute s = a(t x), we have Assume that T > 0; t 2 [0; T ] and f (0) = 0. Let ft 0 ; t 1 ; : : : ; t n g be a partition of the interval [0; T ] such that 0 = t 0 < t 1 < < t n = T . Take the step width h = T n and t i = ih for i = 0; 1; : : : ; n. By (7), we have where b:c is the greatest integer function. Since at i does not have to belong to this partition, the interval [0; at i ] is divided into two subsets as [0; t baic ] and [t baic ; at i ]. Now, let us de…ne I 1 and I 2 as Considering the trapezoidal integration rule with the partition f0; h; 2h; : : : ; [ai] hg of the interval [0; t baic ], we obtain where M k (0) = 0. To calculate the integral I 2 , it can be again used the trapezoidal integration rule with two points (the values of the lower and upper bounds) on the interval [t baic ; at i ]. Hence, By using the equations (9) and (10) into the equation (8), we have Denote the approximate value of M k (t i ) asM k (t i ). Then, for any given k 1 and i = 0; 1; : : : ; n, the values ofM k (t i ) can be recursively calculated as whereM k (0) = M k (0) = 0: Let us reduce this general solution to the numerical solution of integral equations (3) and (4). By (11), the values ofM 3 (t i ) andM 4 (t i ) can be calculated as respectively, where M 1 (t i ) and M 2 (t i ) can be approximately calculated with the help of the formula given in (11) as Note that these formulas in (12) and (13) are given previously by [17] and [15], respectively.

ILLUSTRATIVE EXAMPLES
In this section, we consider gamma, Weibull and lognormal distributions for the …rst interarrival time of the GP to illustrate the numerical method developed in previous section. As indicated in [10], the ratio parameter a satis…es the condition 0:95 a 1:05 for many real data sets …tted by the GP. Further, in the applications of the GP, we mostly encounter with values of a which is less than 1. For this reason, the ratio parameter of the GP is taken as a = 0:95 in each example. It is worth to 740 M . H. PEKALP, A. AYDO ¼ GDU noting that similar results are obtained for the di¤erent values of a. The value of T is chosen to be at least 10E (X 1 ).
In each example, we calculate the approximate values of the skewness and kurtosis of the GP model as ; t 0 respectively.

Example 3. (Lognormal distribution)
Let fN (t); t 0g be a GP with ratio parameter a = 0:95 and assume that the …rst interarrival time X 1 has lognormal distribution LN(0; 1). Taking T = 18 > 10E (X 1 ) = 10e 1=2 , the interval [0; 18] is divided into n = 1800 subintervals with the equal width h = 0:01. Thus, we obtain the approximate values of M 3 (t), M 4 (t) ; S(t) and K(t) in the Table 3 below. It can be concluded from Tables 1-3 that the shape of the distribution of the GP converges to the normal distribution when the value of t gets closer to the mean of the …rst interarrival time.

CONCLUSIONS
Integral equations satis…ed by the third and fourth moment functions of a GP are derived. Further, we present an integral equation for the kth moment function as a generalization of the lower-order moment functions of the GP. In general manner, a numerical method established for solving the integral equation given for the kth moment function is presented. Then, we reduce this general structure to the solutions of the third and fourth integral equations to obtain their solutions. By using the solutions of the integral equation (3) and (4), skewness and kurtosis of the GP model are calculated for some lifetime distributions. According to the numerical calculations, as t gets closer to E(X 1 ), the shape of the distribution of the GP converges to the normal distribution. Note that more precise results can be obtained by taking smaller step width in numerical calculation of the integral equations since the accuracy of the approximation depends on the step width.

Authors Contributions Statement
The authors jointly worked on the results and they read and approved the …nal version of the manuscript.

Declaration of Competing Interests
The authors declare that they have no competing interest.