Monitoring the Zero-Inflated Time Series Model of Counts with Random Coefficient

In this research, we consider monitoring mean and correlation changes from zero-inflated autocorrelated count data based on the integer-valued time series model with random survival rate. A cumulative sum control chart is constructed due to its efficiency, the corresponding calculation methods of average run length and the standard deviation of the run length are given. Practical guidelines concerning the chart design are investigated. Extensive computations based on designs of experiments are conducted to illustrate the validity of the proposed method. Comparisons with the conventional control charting procedure are also provided. The analysis of the monthly number of drug crimes in the city of Pittsburgh is displayed to illustrate our current method of process monitoring.


Introduction
This work is motivated by an empirical analysis and process control of a monthly drug crime series, which contains excess zeros (over 40%) and shows clear serial dependence (see Section 5 for more details). To solve this problem, an appropriate integer-valued model is selected, further, control charts based on this model are developed. Among kinds of integer-valued models, a specific kind featured with first-order integer-valued autoregressive (INAR(1)) models plays an very important role and has been widely studied in the literature. In reality, serial dependence among the count data have been demonstrated to arise extensively in practice, typical examples are infectious disease counts, defect counts and unemployment counts, etc. These data are important indicators of the epidemic study, quality control and economics analysis, and the process monitoring is essential to detect the shifts in them.
The first INAR(1) model proposed by Al-Osh and Alzaid [1] is in the following form where the binomial thinning operator "• is defined by Steutel and Van Harn [2], α • X t−1 = ∑ X t−1 i=1 Y i , α ∈ (0, 1), {Y i } N is a sequence of independent and identically distributed (i.i.d.) random variables with Bernoulli(α) distribution; and {ε t } N is a sequence of i.i.d. random variables, independent of all {Y i }. The INAR(1) model is currently applied in various kinds of real-world problems because of its good interpretability. As one example, we let X t represent the number of patients of an infectious disease in a community at time t, ε t the number of new patients entering the community at time t, and suppose each patient at time t − 1 survives at time t with survival probability α. As for the crime data, α • X t−1 can be considered as the number of re-offendings provoked by X t−1 with probability α. Depending on the nature of this kind of observed data, the INAR(1) models have been modified and 1, 2, · · · , and p ∈ (0, 1), θ > 0. {ε t } N is independent of the past of the solution {X s ; s < t} and the binary sequence {α t }, parameters are also constrained by the condition p/(β + p(1 − β)) < α < 1.
The ZIGINAR RC (1) process is quite suitable for modelling some real-life phenomena in which counted events may survive or vanish with the random survival probability α t . Such series are studied in section 5 with an example of the counts of the drug crimes, where the re-offending rate may be affected by public security situation and financial situation. The mean, variance, and first-order autocorrelation function of the process are, respectively, Obviously the process is characterized by the property of overdispersion, i.e., the variance greater than the expectation. Figure 1 shows some sample paths of simulated ZIGINAR RC (1) processes for θ = 1, 3, 5; p = 0.2, 0.5; α = 0.5, 0.8 and β = 0.3, 0.8. As we can see, the model has larger process mean with larger θ, and larger percentage of zeros with larger p. Following the Theorem 2.1. in Bakouch, Mohammadpour and Shirozhan [21], the ZIGINAR RC (1) model has a unique, strictly stationary solution given by And the probability mass function of {ε t } N is × [αθ(β + p(1 − β))] j [1 + αθ(β + p(1 − β))] j+1 , j = 0, 1, 2, · · · .
where I (j) = 1 for j = 0 and 0 else. It can be deduced that the innovation series {ε t } N is a mixture of three random variables, a degenerate distribution at 0, Geometric(θ/(1 + θ)) Following Theorem 2.1 in Bakouch, Mohammadpour and Shirozhan [21], the ZIGINAR RC (1) model has a unique, strictly stationary solution given by Furthermore, the probability mass function of {ε t } N is where I (j) = 1 for j = 0 and 0 else. It can be deduced that the innovation series {ε t } N is a mixture of three random variables, a degenerate distribution at 0, Geometric(θ/(1 + θ)) and Geometric(αθ(β + p(1 − β))/(1 + αθ(β + p(1 − β)))) distributions with three different mixing portions. The following form is the transition probability of the process {X t } N 0 Some other important probabilistic properties of the process, like spectral density, multistep conditional mean and variance, extreme order statistics, distributional properties of length of run of zeros, have also been discussed in Bakouch, Mohammadpour and Shirozhan [21]. Furthermore, the unknown parameters of the model could be estimated by conditional least squares or maximum likelihood methods.

Monitoring Procedure
In this section, we present a CUSUM chart for monitoring the ZIGINAR RC (1) process. As this process is used to fit the number of crimes, an increase in the process mean usually means a deteriorating public security environment, and an increase in the process correlation usually means more re-offendings. Thus, our purpose is to detect the increasing of both mean shifts and correlation shifts in the ZIGINAR RC (1) process. According to the model properties, the process mean is affected by the parameters θ and p, the correlation is affected by the parameters α and β. Let θ 0 , p 0 , α 0 and β 0 (θ 1 , p 1 , α 1 and β 1 ) denote the in-control (out-of-control) parameters of the processes, and µ 0 , σ 0 , ρ 0 (µ 1 , σ 1 , ρ 1 ) be the corresponding in-control (out-of-control) process mean, standard deviation and first-order correlation.
Scheme (The ZIGINAR RC (1) CUSUM chart). Let {X t } N 0 be a stationary ZIGINAR RC (1) process, the CUSUM statistics C t is defined as: where k is a positive integer constant representing the reference (k µ 0 ). This chart is said to be out-of-control when C t falls outside the control limit h (h ∈ N), that is, C t > h.
The initial value of the CUSUM statistics is set equal to the integer constant c 0 , i.e., C 0 = c 0 with c 0 < h. The performance evaluation of this chart is accomplished based on the average run length (ARL) measures, which is defined as the average number of points to be plotted on the chart until the first out-of-control signal triggers. As {X t , C t } t∈N 0 of the ZIGINAR RC (1) process is a bivariate Markov chain, the Markov chain approach proposed by Brook and Evans [39] is adapted to evaluate the exact ARLs. Though this method has been described in detail in the relevant literature by Weiß [22], Weiß and Testik [23] and Weiß and Testik [24], we briefly introduce this method here for completeness. The reachable control region (CR) of {X t , C t } N 0 is given by Obviously CR has a finite number of elements and could be ordered in a certain manner. The transition probability matrix of {X t , C t } N 0 is Q (p(n, j|m, i)) (n,j),(m,i)∈CR , p(n, j|m, i) P(X t = n, C t = j|X t−1 = m, C t−1 = i) = I (j−max(0,n−k+i)) P(X t = n|X t−1 = m).
The conditional probability that the run length of {X t , C t } N 0 equals r is defined by where (m, i) ∈ CR. Let the vector µ (k) denote the k-th factorial moments that (u (k) ) m,i ∑ ∞ r=1 r (k) p m,i (r) where k 1 and r (k) r × (r − 1) × · · · × (r − k + 1). Then The ARL is obtained as For simplicity we do not repeat the proof methods, see Weiß [22], Weiß and Testik [23] and Weiß and Testik [24] for more details. It is expected that an efficient chart possesses a large in-control ARL (denoted as ARL 0 ) and a small out-of-control ARL. Along with the ARL, we also assess the performance of the charts through the standard deviation of run length (SDRL) suggested by Weiß [22]. The SDRL of the ZIGINAR RC (1) CUSUM chart could again be computed efficiently by applying the Markov chain method. The second order factorial moments u (2) can be determined recursively from the relation (I − Q)u (2) = 2Qu (1) . Then the SDRL is To implement the proposed monitoring scheme, the chart design pairs (h, k) need to be designed in advance. Generally, a fixed ARL 0 value is set to be the target value, and (h, k) is set accordingly. Some guidelines for the choices of them will be given in the next section.

Computation Results
In this section, we evaluate the ZIGINAR RC (1) CUSUM chart performance basing on extensive numerical experiments and presume that the parameters in this model have already been known. In practice, the in-control parameters need to be estimated from the data, as shown in the next section. We search for possible chart designs (integer (h, k) pairs) in order to adjust the ARL 0 close to the target value. Here the target ARL 0 value is set to be 370, which is commonly used in the statistical process monitoring domain. Meanwhile, the values of ARL and SDRL are calculated accurately by the Markov chain method, and we only show the results with two decimal places for simplicity. We first compute ARL 0  Table 1, three important conclusions can be derived. First, it can be observed that when c 0 takes smaller value, the deviation of ARL 0 and SDRL 0 is small. When c 0 takes larger value, there might be a situation where the value of SDRL 0 is significantly greater than the value of ARL 0 (for example, ARL 0 = 409.42, SDRL 0 = 442.13 under (θ 0 , p 0 , α 0 , β 0 , c 0 ) = (1, 0.3, 0.5, 0.8, 6)). Thus, we assume that c 0 = 0 in the following studies to get better robust. Second, as the differences of the values between ARL and SDRL are small when c 0 = 0, we only use ARL as the measure in the following computations to save space. Last, the parameter θ 0 shows a great influence on the selection of control designs (h, k), with a larger θ 0 comes a larger pair of (h, k). Table 1. ARL 0 and SDRL 0 of the CUSUM chart for various θ 0 , p 0 , α 0 , β 0 and c 0 . Due to its simplicity, the conventional Shewhart chart is very popular in monitoring the process shifts. The upper limit for the Shewhart chart is denoted as UCL. For observations, when the value of the process {X t } N 0 exceeds the threshold value UCL (X t > UCL), a fault is declared. Figures 2 and 3 investigate the CUSUM method preliminarily by comparing it with the Shewhart method. In both of these figures, we assume that the in-control parameters are θ 0 = 2, p 0 = 0.2, α 0 = 0.5 and β 0 = 0.5, which are selected based on the real drug crime data in Section 5. According to these parameters, the CUSUM chart designs can be determined, respectively, as h = 31, k = 2 (corresponding ARL 0 = 383.74); h = 19, k = 3 (ARL 0 = 396.12); h = 14, k = 4 (ARL 0 = 373.27); h = 11, k = 5 (ARL 0 = 370.77); h = 9, k = 6 (ARL 0 = 394.03). Furthermore, the Shewhart chart limit UCL = 13 (ARL 0 =381.31) can be used. It should be noted that two types of changes are considered in Figure 2, which both lead to the upward mean shifts. The first type of changes occurs only in the parameter θ, with other parameters invariant, the results are listed in Figure 2a CUSUM charts under most shifts, while the Shewhart chart performs worst among them. For the upward correlation shift scenarios, ARL values under two types of the parameter changes are displayed in Figure 3. The first one considers changes only in the parameter α, and the second one considers changes only in the parameter β. For each scenario, the Shewhart chart performs increase ratio of ARL with the increase of first-order correlation ρ X , and the CUSUM chart has the better behaviour in the figure. In a comprehensive view, the conventional Shewhart chart is insensitive for upward mean shifts caused by changes in parameter p, and fails to detect shifts in the correlation. While the proposed CUSUM chart could overcome these limitations and has superiorities in various coefficient shifts compared with the Shewhart chart. From the figures, we can also conclude that the smaller the value of k, the more sensitive the CUSUM chart is. As the constraint k µ 0 is required to make the chart reasonable, it is natural to recommend k = µ 0 (the smallest integer no less than µ 0 ), then we aim to select the value of h such that ARL 0 is close to 370. Now the computations of the CUSUM chart are extended to general cases with designs of experiments as follow. Shewhart chart performs increase ratio of ARL with the increase of first-order correlation ρ X , and the CUSUM chart has the better behaviour in the figure. In a comprehensive view, the conventional Shewhart chart is insensitive for upward mean shifts caused by changes in parameter p, and fails to detect shifts in the correlation. While the proposed CUSUM chart could overcome these limitations and has superiorities in various coefficient shifts compared with the Shewhart chart. From the figures, we can also conclude that the smaller the value of k, the more sensitive the CUSUM chart is. As the constraint k µ 0 is required to make the chart reasonable, it is natural to recommend k = µ 0 (the smallest integer no less than µ 0 ), then we aim to select the value of h such that ARL 0 is close to 370. Now the computations of the CUSUM chart are extended to general cases with designs of experiments as follow.  In Tables 2,3,4, we focus on situations that there are increasing shifts in process mean, and the correlation remains the same. Each in-control parameter has three levels: We consider the case when the changes only occur in θ, this is the most common case. The out-ofcontrol process mean is µ 1 = µ 0 + δσ 0 , the shift size δ considers potential values in set {0.5, 1, 1.5, 6}. The usual relative deviation (in %) in ARL is defined as dev (%) = 100%× (ARL-ARL 0 )/ARL 0 (Weiß and Testik [24]). From Tables 2,3,4, we can conclude that the ZIGINAR RC (1) CUSUM chart performs quite well in detecting upward mean shifts for all scenarios. For the small shift size of δ = 0.5, the CUSUM chart is efficient with the Shewhart chart performs increase ratio of ARL with the increase of first-order correlation ρ X , and the CUSUM chart has the better behaviour in the figure. In a comprehensive view, the conventional Shewhart chart is insensitive for upward mean shifts caused by changes in parameter p, and fails to detect shifts in the correlation. While the proposed CUSUM chart could overcome these limitations and has superiorities in various coefficient shifts compared with the Shewhart chart. From the figures, we can also conclude that the smaller the value of k, the more sensitive the CUSUM chart is. As the constraint k µ 0 is required to make the chart reasonable, it is natural to recommend k = µ 0 (the smallest integer no less than µ 0 ), then we aim to select the value of h such that ARL 0 is close to 370. Now the computations of the CUSUM chart are extended to general cases with designs of experiments as follow.  In Tables 2,3,4, we focus on situations that there are increasing shifts in process mean, and the correlation remains the same. Each in-control parameter has three levels:  (ARL-ARL 0 )/ARL 0 (Weiß and Testik [24]). From Tables 2-4, we can conclude that the ZIGINAR RC (1) CUSUM chart performs quite well in detecting upward mean shifts for all scenarios. For the small shift size of δ = 0.5, the CUSUM chart is efficient with the minimum 86.38% drop of ARL and the maximum 91.23% drop of ARL. Take larger δ for another illustration (δ = 1), the drop of ARL is at least 92.81%, and up to 96.03% at most. It can also be obtained that δ has to be at least 6 to get an immediate signal with the out-of-control ARL closer to 1. In addition, extensive computation results show that the in-control parameters (θ 0 , p 0 , α 0 , β 0 ) have little effect on the better performance of the CUSUM chart to detect the mean shifts. Table 2. ARL profiles of the CUSUM chart versus mean shifts under θ 0 = 1.

Process Parameters
ARL (dev (%) )  Table 3. ARL profiles of the CUSUM chart versus mean shifts under θ 0 = 3.  The computation study in Tables 5 and 6 concerns the upward shifts in the process correlation. Two levels are accepted for each in-control parameter: θ 0 = {1, 5}; p 0 = {0.1, 0.2}; α 0 = {0.5, 0.6}; β 0 = {0.7, 0.8}. Two types of out-of-control pattern are considered here for comprehensive investigation, the first type is that the upward changes only exist in α (shown in Table 5), and the second type is that the downward changes only exist in β (shown Table 6). In Table 5, the shifts in the magnitude δ α = α 1 − α 0 are from the set {0.1, 0.2, 0.3}. The results imply that the performance of CUSUM chart fluctuates greatly in detecting correlation shifts caused only by α. To be specific, when δ α is 0.3, dev (%) ranges from −6.81% to −23.68%. Another finding based on the design of experiments in Table 5 is that both a smaller θ 0 and a smaller β 0 could slightly improve detection efficiency, while p 0 and α 0 could not. In Table 6, the shifts magnitude δ β = β 1 − β 0 are from the set {−0.1, −0.2, −0.3}. From Table 6, we can see the CUSUM chart is more efficient in detecting the correlation shifts caused by β. As the absolute value of δ β gets bigger, the decreasing proportion of ARL gradually increased. When δ β = −0.3, dev (%) ranges from −21.3% to −40.39%. Meanwhile, we can further conclude that a smaller θ 0 and a larger α 0 often lead to better chart performance, and p 0 , β 0 have little influence. Based on all the analysis above in this paragraph, we can further conclude that a smaller θ 0 and a larger value of initial correlation ρ 0 are helpful to detect the correlation shifts. Furthermore, that we cannot get an immediate signal when only correlation shifts occur. Table 5. ARL profiles of the CUSUM chart versus correlation shifts caused by α.  Table 6. ARL profiles of the CUSUM chart versus correlation shifts caused by β.

Analyses of Drug Crime Count Time Series
In this section, we present a case study of crime count data in Pittsburgh. The data set contains multiple crime types, such as arson, drink-driving, robbery and so on. Monitoring of crime data is needed not only for early warnings of the organised crime, but also for assessments of the social security environment. For the crime data, the readers can download it from the Forecasting Principles site (http://www.forecastingprinciples.com, accessed on 20 March 2021), or email to the corresponding author to access. The subset we analyse is a monthly drug use count data collected from the 56th police car beat, which contains 144 observations from January 1990 to December 2001. There are 67 zeros in this drug use data (the proportion up to 46.53%), which have the greatest proportion among the other values for the data series. The sample mean, variance and first-order autocorrelation of the data are 1.7153, 6.4289 and 0.3886, respectively, which show strong overdispersion and autocorrelation. The sample path and the histogram of the series are in Figure 4. The histograms of estimated ZIG distribution, estimated Geometric distribution and estimated Poisson distribution are also given in Figure 4b, which indicate that the ZIG marginal is the most appropriate to describe the data. The sample autocorrelation function (ACF) and the sample partial autocorrelation function (PACF) in Figure 5 reveal that the series most likely comes from an AR-type process of order 3. While our intention is to illustrate the implementation of the proposed control chart, we will employ the first-order INAR models that are widely studied and applied in the literature. The consideration of more complex models will be left for future study.   In this section, we present a case study of crime count data in Pittsburgh. The data set contains multiple crime types, such as arson, drink-driving, robbery and so on. Monitoring of crime data is needed not only for early warnings of the organised crime, but also for assessments of the social security environment. For the crime data, the readers can download it from the Forecasting Principles site (http://www.forecastingprinciples.com), or email to the corresponding author to access. The subset we analyse is a monthly drug use count data collected from the 56st police car beat, which contains 144 observations from January 1990 to December 2001. There are 67 zeros in this drug use data (the proportion up to 46.53%), which have the greatest proportion among the other values for the data series. The sample mean, variance and first-order autocorrelation of the data are 1.7153, 6.4289 and 0.3886 respectively, which show strong overdispersion and autocorrelation. The sample path and the histogram of the series are in Figure 4. The histograms of estimated ZIG distribution, estimated Geometric distribution and estimated Poisson distribution are also given in Figure 4 (b), which indicate that the ZIG marginal is the most appropriate to describe the data. The sample autocorrelation function (ACF) and the sample partial autocorrelation function (PACF) in Figure 5 reveal that the series most likely comes from an AR-type process of order 3. While our intention is to illustrate the implementation of the proposed control chart, we will employ the first-order INAR models that are widely studied and applied in the literature. The consideration of more complex models will be left for future study.  Except for the ZIGINAR RC (1) model, some competitive models are also applied to the time series, such as Poisson INAR(1) (Al-Osh and Alzaid [1]), GINAR(1) (Alzaid and Al-Osh [7]), ZINAR(1) (Jazi, Jones and Lai [9]), ZMGINAR(1) (Barreto-Souza [11]), NGINAR(1) (Ristić, Bakouch and Nastić [12]), ZIMINAR(1) (Li, Wang and Zhang [14]). The Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC) are suggested to evaluate these models. Numerical results in Table 7 show that the ZIGINAR RC (1) process has the best overall performance, compared with its competitors. Therefore, we assume that the drug crime data is from the ZIGINAR RC (1) model and the estimated parameters in Table 7 are used in the process control procedures. Based on the computation results in Section 4, the CUSUM chart with designs h = 34, k = 2 (corresponding ARL 0 = 364.44) is the best choice, which is shown in Figure 6a. For comparison, we also present CUSUM control charts with designs h=15, k = 4 (ARL 0 = 358.40) in Figure 6b, designs h = 12, k = 5 (ARL 0 = 372.28) in Figure 6c, and the Shewhart chart with control limit UCL = 13 (ARL 0 = 340.25) in Figure 6d. We observe that all the CUSUM charts give out-of-control signals, while there is no outliers in the Shewhart chart. Because the Shewhart chart has been proved to be less effective than the CUSUM chart, the drug crime data set seems to be out-of-control with increasing mean shifts or increasing correlation shifts, and some investigation should be done for further explanation. The CUSUM control charts with three designs also display different detection efficiencies. The CUSUM chart with k = 2 first signals at t = 131 following with continuous alarms as t increases. The signals of the CUSUM chart with k = 4 are first given at t = 133, then go back below the control limit over a period of time, and come again at t = 141. While the outlier of the CUSUM chart with k = 5 occurs at t = 144. The analysis above proves again that the CUSUM chart design k = µ 0 is the most effective in practice.  that the drug crime data is from the ZIGINAR RC (1) model and the estimated parameters in Table 7 are used in the process control procedures. Based on the computation results in Section 4, the CUSUM chart with designs h = 34, k = 2 (corresponding ARL 0 = 364.44) is the best choice, which is shown in Figure 6 (a). For comparison, we also present CUSUM control charts with designs h = 15, k = 4 (ARL 0 = 358.40) in Figure 6 (b), designs h = 12, k = 5 (ARL 0 = 372.28) in Figure 6 (c), and the Shewhart chart with control limit UCL = 13 (ARL 0 = 340.25) in Figure 6 (d). We observe that all the CUSUM charts give out-of-control signals, while there is no outliers in the Shewhart chart. Because the Shewhart chart has been proved to be less effective than the CUSUM chart, the drug crime data set seems to be out-of-control with increasing mean shifts or increasing correlation shifts, and some investigation should be done for further explanation. The CUSUM control charts with three designs also display different detection efficiencies. The CUSUM chart with k = 2 first signals at t = 131 following with continuous alarms as t increases. The signals of the CUSUM chart with k = 4 are first given at t = 133, then go back below the control limit over a period of time, and come again at t = 141. While the outlier of the CUSUM chart with k = 5 occurs at t = 144. The analysis above proves again that the CUSUM chart design k = µ 0 is the most effective in practice.

Conclusions
In this paper, we have made contributions on monitoring the zero-inflated autocorrelated count data which can be described by an INAR(1) process with random coefficient. ARL and SDRL are adopted to be the measures and calculated by the Markov chain ap-

Conclusions
In this paper, we have made contributions on monitoring the zero-inflated autocorrelated count data which can be described by an INAR(1) process with random coefficient. ARL and SDRL are adopted to be the measures and calculated by the Markov chain approach. The design parameter k in the CUSUM chart is proved to have great influence on monitoring efficiency, and the smallest integer no less than µ 0 is the recommended value of k. The proposed ZIGINAR RC (1)CUSUM chart is proved to have superiorities in various coefficient shifts compared with the conventional Shewhart chart. Computation results also show that the CUSUM chart performs quite well in detecting upward mean shifts, and shows fluctuation in detecting upward correlation shifts. Based on the design of experiment, we also find that a larger value of initial correlation ρ 0 is helpful to detect the correlation shifts. An immediate signal occurs after large upward mean shifts, while does not occur when only correlation shifts exist.
There are some possible topics for our future research. First, we can consider a different monitoring scheme for the ZIGINAR RC (1) process and conduct a comparison study. Second, we can explore the monitoring of p-th random coefficient INAR model, which is suitable for the count data with high order dependence. Third, we can study a multivariate INAR(1) process with random coefficient to continuously monitor the serial correlated counts that we are interested in.