Control charts for monitoring a Poisson hidden Markov process

Monitoring stochastic processes with control charts is the main field of application in statistical process control. For a Poisson hidden Markov model (HMM) as the underlying process, we investigate a Shewhart individuals chart, an ordinary Cumulative Sum (CUSUM) chart, and two different types of log‐likelihood ratio (log‐LR) CUSUM charts. We evaluate and compare the charts' performance by their average run length, computed either by utilizing the Markov chain approach or by simulations. Our performance evaluation includes various out‐of‐control scenarios as well as different levels of dependence within the HMM. It turns out that the ordinary CUSUM chart shows the best overall performance, whereas the other charts' performance strongly depend on the particular out‐of‐control scenario and autocorrelation level, respectively. For illustration, we apply the HMM and the considered charts to a data set about weekly sales counts.


| INTRODUCTION
Hidden Markov models (HMMs) are commonly used in practice when being concerned with serially dependent processes ðX t Þ N 0 , where N 0 = f0, 1, …g, especially when the process is discrete-valued, such as a count process or a categorical process, see the literature. 1,2 The basic idea of HMMs is the presence of a sequence ðQ t Þ N 0 of latent states (out of finitely many ones), which themselves follow a Markov chain (MC) model; the observations X t are then emitted depending on the current state Q t . More details on HMMs are presented in Section 2. In the present work, we solely focus on the case of ðX t Þ N 0 being a count process; so the range of X t is equal to the set N 0 of nonnegative integers or a subset thereof. But it would also be possible to adapt our methods to HMMs with a different type of range for the observations X t . More precisely, we assume that the process is stationary if running in control and follows a certain HMM for counts. But it may happen that the process leaves this model after some time (it thus runs out of control), and the aim is to detect such a process change as early as possible.
In order to monitor such stochastic processes, the discipline of statistical process control (SPC) provides tools for detecting changes in the process compared to the given in-control model, namely, control charts. The most basic types of control charts are Shewhart control charts, named after Walter A. Shewhart  in view of his pioneering work in this area. They declare the monitored process to be in control as long as the values of its plotted statistics, for example, the individual observations X t themselves, lie between two predefined control limits. Otherwise, one concludes that the process might be out of control. Shewhart control charts are very useful to detect assignable causes resulting, for example, in large shifts in the process mean. At the same time, Shewhart charts exhibit a considerable disadvantage: if small sustained shifts are of interest, these charts are known to be very ineffective, because they ignore the information embodied by earlier data (i. e., they lack memory). When more sensitivity to smaller shifts is preferred, Cumulative Sum (CUSUM) control charts constitute a better choice, as they accumulate information from previous samples, see Montgomery 3 for details. For ease of presentation, we concentrate on the case of using upper-sided charts for monitoring the considered count process ðX t Þ N 0 , but our approaches presented below are easily adapted to lower-sided or two-sided charts as well, by adapting the idea of Yontay et al. 4 As mentioned before, HMMs have been widely applied in practice, see Zucchini et al 1 for references, but there seems to be only little work on SPC methods applied to HMMs. Alshraideh and Runger 5 propose to monitor the residuals obtained from the considered HMM by using a Shewhart chart. Fuh 6 derives asymptotic optimality for some CUSUM-type sequential procedures, and Fuh and Mei 7 study approaches related to the generalized likelihood ratio (GLR) scheme for a two-state HMM. Actually, one more frequently finds applications of HMMs in a surveillance context, see the references in Sparks. 8 Such an example is given by Rafei et al, 9 who monitor weekly counts of tuberculosis cases. In contrast to typical SPC applications, where an in-control model (and a corresponding control chart design) is developed for the "common phase" in order to detect any out-of-control situation that violates the in-control assumptions, Rafei et al 9 propose a two-state HMM with two states representing the common phase (in-control) and a particular type of uncommon phase (out-of-control), that is, in-control and out-of-control are part of one model. An analogous definition of the considered HMM can be found in Simões et al, 10 where diesel engines are monitored and five hidden states are defined that correspond to a degradation of the system. Two of these hidden states are interpreted as fault states which require an immediate maintenance action. Similar HMM designs are commonly used in machining and maintenance applications, where the hidden states express increasing tool wear, see, for example, other works. [11][12][13] In summary, an HMM could be used within control charting in two ways. Either it serves as an explicit change point model or it is the appropriate statistical model for the in-control process. In the first case, however, we assume to know the probability law for reaching the change point while utilizing very simple probability models for the in-control and one or a few more out-of-control phases. In the second case, we make use of the full flexibility of the HMM to describe in-control processes exhibiting particular patterns. Here, we utilize the HMM in its primary sense (or as a parsimonious time series model) for characterizing the in-control process and its potential deviations.
As an illustrative data example, we consider a time series of counts of weekly sales of a soap product, 14 where the hidden states can be interpreted as different (in-control) levels of demand. The monitoring of demand is important, for example, within an inventory management system. The aim is to detect changes compared to the in-control HMM such as shifts within the demand states or distributional changes between the demand states, to be able to act accordingly in presence of changed sale activities. Further details on these data are provided in Section 2, where we also review the definition and important properties of HMMs. Section 3 then presents the control charts for monitoring HMMs considered in this work, namely, the ordinary c and CUSUM chart as well as two different types of log-LR CUSUM charts. The ARL performances of these charts, for diverse out-of-control scenarios, are analyzed in Section 4. Section 5 then picks up the data example from Section 2 and illustrates the application of our control charts in practice. Section 6 discusses the topic of state-dependent control charts, and Section 7 concludes the article.

| HIDDEN MARKOV MODELS AND APPLICATIONS
An HMM is defined to be a bivariate process ðX t , Q t Þ N 0 , where the X t are the observable random variables (counts in our work), and the Q t are the hidden states (latent states) with a finite qualitative range Q, see the literature 1,2 for comprehensive treatments. For simplicity, we label these states by integer numbers, that is, Q = f0, …, d Q g with some d Q 2 N, where d Q is typically a small number. The hidden states are assumed to follow a homogeneous MC. Given the state process, the observation process is generated serially independently, with its probability mass function (PMF) depending solely on the current state Q t .
To put it in a nutshell, a (basic) HMM for ðX t , Q t Þ N 0 satisfies the following three conditions: the observation equation PðQ t = qjQ t − 1 = r, …Þ = PðQ t = qjQ t − 1 = rÞ = a qjr , 8q, r 2 Q: ð2:3Þ In addition, the initial distribution of Q 0 is hereafter assumed to satisfy the invariance equation Aπ = π = ðπ 0 , …, π d Q Þ > , where A = (a qjr ) q,r denotes the transition matrix of the hidden states, and π a vector of marginal probabilities. In other words, the MC ðQ t Þ N 0 given by (2.3), and thus, the whole HMM ðX t ,Q t Þ N 0 , is assumed to be stationary with P(Q t = q) = π q for all t 2 N 0 and all q 2 Q. For a vivid description of an HMM's data-generating mechanism, see Alshraideh and Runger. 5, Section 2 In our performance investigations in Section 4, we ensure this stationarity by defining ðQ t Þ N 0 as a first-order discrete autoregressive (DAR(1)) process, 2 that is, it follows the recursion: ð2:4Þ where the α t are independent and identically distributed (i. i. d.) Bernoulli random variables with "success probability" ϕ, and where the innovations ϵ t are i. i. d. according to π with range Q. Both processes, ðα t Þ N and ðϵ t Þ N , are assumed to be independent of each other and of (Q s ) s < t . This DAR(1) model leads to a parsimoniously parametrized MC model with marginal distribution π and transition matrix The parameter ϕ controls the extent of serial dependence, with larger values corresponding to stronger dependence. We restrict ourselves to time-homogeneous state-dependent distributions; therefore (2.1) extends to P(X t = x j Q t = q) =: p(xjq) for all t 2 N 0 . In applications, mainly parametric distributions are assumed for p(� jq). For illustration, we only consider the Poisson HMM in the sequel, which means that the state-dependent distributions of X t are assumed to be Poisson distributions. So being in state q 2 Q at time t, X t follows Poi(λ q ) with some λ q > 0: Without loss of generality but to simplify the interpretation, we assume that , but the unconditional distribution of X t is overdispersed as a Poisson mixture. It should be noted that any other count model could be used for the HMM as well, also different models for different states. Recently, Adam et al 15 even developed an estimation approach for a nonparametric count HMM. A brief summary of further stochastic properties as well as approaches for the parameter estimation, forecasting, and decoding of the hidden states is provided in Appendix A. More detailed discussions can be found in previous studies. 1,2 These books also provide R codes for doing these computations. In addition, ready software implementations are available for common types of computational software, for example, the "HiddenMarkov" package in R, 16 the command HiddenMarkovProcess in Wolfram Mathematica™, or the hmm… commands in MATLAB's Statistics and Machine Learning Toolbox™. At this point, a practitioner might ask how to recognize at all that an HMM could be used for count time series modeling. In our opinion, one should recall the famous words "All models are wrong but some are useful" 17 here. HMMs are very flexible and can be adapted to many types of data and dependence structure. In fact, Zucchini et al 1 even recommend them as "general-purpose models for time series" (p. 5). In this spirit, Alshraideh and Runger 5 use HMMs as a competitor to the ARIMA models because of their simplicity. Besides their flexibility in capturing quite different time series properties, two further advantages should be stressed: the hidden states are often well interpretable, and a lot of software implementations for HMMs are readily available, also see the previous discussion. So there are several pragmatic reasons for using HMMs in the context of count time series modeling. Nevertheless, any application of an HMM should be complemented by a careful analysis of model adequacy, see, for example, in Weiß. 2, Sections 2.4 and 5.2

| A numerical example
To illustrate the HMM's structure as well as common types of statistical inference for HMMs, we first present a simple numerical example (later in Section 2.2, we also discuss a real-data example). The numerical example is motivated by a manufacturing situation with two hidden states (i. e., d Q = 1). For instance, we may count the number of tool changes at a machine within specified time intervals (the HMM's observations, X t ). The two possible hidden states are caused by the incoming raw material (two quality levels, good vs. medium), and the state-dependent Poisson distribution (2.6) has a lower mean for the good-quality batches than for the medium-quality ones. This setup is similar to the coal mill example in Kisi� c et al, 13 where the wear of the grinding table depends on the quality of the processed coal. Another example for a two-state HMM would be a production environment with two load conditions (normal load vs. overload). The number X t of produced items in the tth time interval has a larger mean under normal-load operation, because overload operation leads to a throttling of the machine output. Similar to the material example, we cannot directly observe the actual load condition, so this constitutes the HMM's hidden state Q t .
We assume the following parametrization of the two-state Poisson HMM: Note that the transition matrix A corresponds to a DAR(1) model with ϕ = 0.7, recall (2.5). The hidden MC features inertia in the sense that it tends to stay in its present state with a probability of 0.8 and 0.9, respectively. The overall probability for being in state "0" is 1/3, and the corresponding state-dependent Poisson distribution has the mean 2 (whereas state "1" goes along with mean 5). Using the formulae given in Appendix A, we compute the marginal mean as 4 and the variance as 6, so the HMM exhibits 50% overdispersion. The autocorrelation function (ACF) takes the values 0.233, 0.163, 0.114,… at lags 1, 2, 3, … For illustration, we simulated a time series of length T = 100 from this HMM, see the left plot in Figure 1. There, we also show the true values of the hidden states as empty circles, plotted at the value λ q if the respective state equals q. Then, we fitted a two-state Poisson HMM to the data via maximum likelihood (ML) estimation (see Appendix A for details), leading to the estimatesÂ Next, the fitted model was applied to uncover the hidden states. Here, one distinguishes between global decoding of all hidden states, q 1 , … , q 100 , and local decoding, where only a single hidden state Q t is determined. The result of a global decoding (using the Viterbi algorithm as described in Appendix A) is shown by grey squares in Figure 1 (plotted as resulting meanλq t ). Nine of the hidden states were decoded incorrectly, most of them between t = 54, …, 60. Finally, the fitted model was used for computing the h-step-ahead forecasting distributions P(X 100 + h = x j x 1 , … , x 100 ); the ones for the forecast horizons h = 1,5 are plotted in the right part of Figure 1. Defining the respective point forecast as the count value maximizing the forecast PMF (mode forecast value), the point forecast for time 101 is 4, and the one for time 105 is 3.

| Application to sales counts time series
Although manufacturing is probably their most well-known field of application, control charts are also used in many other areas such health surveillance or service industries. 18 As already indicated in Section 1, HMMs are frequently applied in these fields. A health-related example is described by Sebastian et al, 19 where a time series of monthly Vibrio cholerae (VC) counts is modeled by a three-state Poisson HMM. There, the three states express "mild," "moderate," and "severe" VC epidemic, and they are related to climate conditions. In what follows, we consider a data example that can be attributed to the field of inventory management. There has been a lot of research activity on the use of control charts for inventory management systems during the last 25 years, see, for example, other works. [20][21][22] There, control charts are applied, among others, to demand time series to detect changes in the operating environment, which, in turn, would affect the inventory system's performance. Related to such applications, we analyze a count time series regarding the demand for a certain soap product in a supermarket. 14 It contains 242 data points, x 1 , … , x 242 , each observation representing the weekly number of sales of the considered soap product. For this time series, it has been established that a stationary three-state Poisson HMM (i. e., d Q = 2) provides an excellent fit, see MacDonald and Zucchini. 14 Therefore, we shall use this data set to further illustrate the interpretation and model fitting of HMMs and to design and apply appropriate control charts for process monitoring, where the latter is done in Section 5.
The nine model parameters (transition matrix A and state-dependent Poisson parameters λ q , q = 0,1,2) are again estimated by ML estimation, that is, by a direct numerical maximization of the likelihood function with θ and P(x) denoting the parameter vector and the diagonal matrix PðxÞ : = diag pðxj0Þ, …,pðxjd Q Þ ð Þ 2 ½0; 1� ðd Q + 1Þ × ðd Q + 1Þ , respectively (also see Appendix A). The corresponding implementation in R can be found in the supplemental materials (together with all other codes used in this article). The state-dependent means of the fitted Poisson HMM (see (2.6)) areλ = ðλ 0 ,λ 1 ,λ 2 Þ≈ð3:74,8:44,14:93Þ, whichfrom a practical point of viewcan be interpreted as a low (state 0), medium (state 1) and high (state 2) level of demand for the particular soap product. Together with the estimates for the states' PMF,π = ðπ 0 ,π 1 ,π 2 Þ > ≈ð0:722, 0:220, 0:058Þ > , they yield a model mean of 5.42 and a model variance of 14.72, both being very close to the corresponding sample values 5.44 and 15.40, respectively. As before, the fitted model was also used for global decoding of the complete series of hidden states, q 1 , … , q 242 , by using the Viterbi algorithm (see Appendix A). Figure 2 shows the data and the sequence of states (left-hand side) as well as the sample and model ACF against increasing time lags (right-hand side). Clearly, the sample ACF indicates a significant degree of positive serial dependence, which is almost exactly captured for the first four lags by the fitted HMM. As the time series in Figure 2 nicely shows, state 0 is predominant with an overall probability of 72.2% for being in this state, whereas the occurrence of state 2 seems to be most unlikely (5.8%). Also when looking at the estimatesâ 3j1 ,â 3j2 in one can conclude that a period of high demand seems to start very rarely. However, once the demand for soap reaches state 2, it tends to remain at this level instead of changing directly. Furthermore, it is interesting to note that the fitted model nearly excludes the transition from state 2 to state 0 (the displayed values are roundings), which means that we would not expect a decline in demand to be that radical.

| MONITORING HMMS WITH C CHART, CUSUM CHART, AND LOG-LR CUSUM CHART
To detect changes in, for example, the sales activity as early as possible, we shall from now on concentrate on control chart approaches for count processes following an HMM. Because in most applications, the detection of increases in the counts' mean is particularly relevant (e. g., to request supply of sales products, or to detect an outbreak of an epidemic disease), we focus on upper-sided control charts; so the investigated charts require only one (that is, an upper) control limit. But it would easily be possible to adapt our monitoring concepts to the lower-sided or two-sided case. As the monitoring approaches, we consider the c chart as a common benchmark chart on the one hand, and different types of CUSUM charts on the other hand, providing an improved sensitivity towards small process changes thanks to an inherent memory.
The c chart is a basic Shewhart-type control chart, 3 which simply plots the counts X t as they arrive in time. It triggers an alarm if its specified control limit u > 0 is violated, that is, if X t > u. Because the monitored data are autocorrelated, we cannot use simple 3-σ or quantile limits but have to adapt the control limits according to the actual serial dependence structure (so a "modified Shewhart chart" in the sense of Schmid 23 ). Like for the remaining charts, this is done based on average run length (ARL) considerations. The c chart is often well suited for detecting large shifts, and it also serves as a benchmark for performance analysis in this work, analogously to, for example, the literature. 24,25 The statistics of the (standard) upper-sided CUSUM chart are calculated as with k > 0 and c 0 ≥ 0 as the reference and starting value, respectively. The parameter k is commonly chosen larger than but close to the in-control mean μ 0 . A value c 0 > 0 is referred to as a Fast Initial Response (FIR) feature. (Generally, employing a FIR feature may help to detect an initial out-of-control situation more quickly; at the same time, it may also cause more false alarms when the process is initially in-control. Fixing this leads to slightly delayed detection of late changes.) If C t exceeds the CUSUM chart's control limit h > 0, an alarm indicates that the process might be out of control. In view of an exact performance evaluation, we choose rational or even integer design parameters k, h, c 0 . The ordinary CUSUM recursion (3.1) can be motivated by the log-likelihood ratio test applied to i. i. d. Poisson counts. If explicitly considering the log-likelihood ratio for HMMs, a more complex type of CUSUM chart is obtained, the log-LR CUSUM control chart. To prevent a possible numerical underflow in later computations, the following recursive scheme (similar to the one recommended by Zucchini et al 1 ) is utilized for calculating the likelihood function given the count time series x 1 , … , x T (also see Appendix A): The log-likelihood function then equals The log-likelihood ratio requires another parameter vector which targets a certain out-of-control scenario, say θ 1 , whereas θ 0 denotes the in-control parameter vector. The log-likelihood ratio reads as ð3:3Þ where "0" and "1" in the index of w t,� refer to the in-control and anticipated out-of-control parameters, respectively, used for the computations in (3.2). Note that the HMM allows for different θ yielding the same process mean μ. Therefore, in Section 4, we investigate two log-LR CUSUM charts with different θ 1 : one, where only the state-dependent means λ q change compared to θ 0 (thus P(x) in (3.2)); and another one, where the PMF π of the hidden states Q t changes (thus the transition probabilities A in (3.2)). We refer to these specific charts as the log-LR λ and log-LR π CUSUM chart, respectively. However, both parameter vectors θ 1 result in the same out-of-control mean μ 1 ; further explanations are provided in Section 4. The log-LR CUSUM chart is now defined as with ℓR t = ln(w t,1 ) − ln(w t,0 ). Again, if ℓC t > ℓh > 0 with ℓh as the control limit, this indicates that the process might be out of control. It is worth mentioning that Fuh and Mei 7 also consider a log-LR CUSUM statistic but applied to a twostate HMM with normal state-dependent distributions (thus variables data).

| PERFORMANCE EVALUATION OF C CHART, CUSUM CHART, AND LOG-LR CUSUM CHART
We assess the effectiveness of all investigated charts by their ARL. In the field of SPC, the ARL is a common performance measure, which indicates the expected number of points plotted on a control chart before a possible out-ofcontrol condition is signaled. Thus, a control chart's ARL should be as large as possible while the monitored process is in control, whereas the out-of-control ARL should be preferably small. In the following, we restrict our investigations to the most popular ARL concept, the zero-state ARL, for evaluating both the in-control and out-of-control performance.
(Although there exist a few more ARL concepts, see Knoth, 26 the zero-state ARL is particularly meaningful in the CUSUM case as it expresses some kind of worst-case behavior.) In order to compute the ARL values of the c and the CUSUM chart, the MC approach of Brook and Evans 27 is utilized. This is possible, because ðX t , Q t Þ N 0 and ðX t , Q t , C t Þ N 0 constitute a bivariate and trivariate MC, respectively, see Appendix B for further details. Because the processes ðX t Þ N 0 and ðC t Þ N only take discrete values, applying the MC approach yields exact numerical results. ARL values of the log-LR charts are obtained by simulations (with 10 6 replications per ARL value); so the results for the log-LR charts are only approximations, nevertheless sufficiently accurate ones. The charts' design parameters are chosen such that the in-control ARL (hereafter signified by adding the subscript "0") is between 200 and 300. In order to be able to compare the investigated charts, the differences between their ARL 0 values should be as small as possible. Because the c chart has only one parameter affecting its ARL, that is, its upper control limit u 2 N, the chart is most inflexible in terms of adjusting the ARL 0 to a given target value. More flexibility is shown by the CUSUM chart with two design parameters (neglecting a FIR feature, i. e., setting c 0 := 0), whichin contrast to the c chart's control limit ucan even take noninteger values for a further fine-tuning of the ARL performance. For practical reasons concerning software implementations (see the corresponding R-code in the supplemental materials) and without loss of generality, the CUSUM parameters are chosen to be rational, that is, k,h 2 Q + , and to have the same denominator. So as the first step of our performance study, an appropriate value for u is chosen. Then the CUSUM chart's parameters are determined with respect to the c chart's ARL 0 , and finally, the control limits of both log-LR charts are selected such that their ARL 0 values locate between the ones of the c and the CUSUM chart. Usually, one can find more than just one relevant tuple (k, h) as CUSUM chart design having nearly the same ARL 0 ; the same holds for both log-LR charts. For designing the CUSUM chart, we first have to set the reference value k, which we decided to be close to the in-control mean μ 0 (analogous to the recommendation of Weiß, 24 who investigated CUSUM charts for an alternative integer-valued model with autocorrelation, namely, first-order integer-valued autoregressive (INAR(1)) processes). Subsequently, an appropriate value for h was determined. For the log-LR charts, a certain out-ofcontrol mean μ 1 must be predefined, then a corresponding control limit ℓh is identified. Now, for comparison purposes, the log-LR charts should aim for the detection of smaller (positive) shifts as well, that is, μ 1 should be reasonably close to μ 0 (one may use the formula k = ðμ 1 − μ 0 Þ=ðlnμ 1 −lnμ 0 Þ in Lucas, 28 which was derived for i. i. d. Poisson counts, as a rough guide for adjusting the different CUSUM concepts). In the following, sustained (positive) shifts in the marginal mean μ of the observable count process are considered. These shifts are obtained in three different ways: • all state-dependent means λ q are multiplied by the same factor simultaneously (see the two upper graphs in Figure 3); • only one of the state-dependent means λ q is multiplied by a certain factor (see Figure 4); • the probability mass within the states' PMF π is shifted towards the state with the largest mean (i. e., towards state d Q ), whereas the remaining probabilities are reduced by the same relative amount (see the lower graphs in Figure 3).
All scenarios are adjusted to give the same increase δ = μ − μ 0 in the marginal mean μ (in the last scenario, μ is bounded by the value λ d Q of the largest state-dependent mean).
Our presented results specifically refer to an in-control three-state Poisson HMM with λ=(λ 0 , λ 1 , λ 2 ) = (1, 2, 5) and π=(π 0 , π 1 , π 2 ) > = (0.5, 0.35, 0.15) > , leading to μ 0 = 1.95. To easily control the dependence level while keeping π fixed, we use the DAR(1) model from Section 2, with the parameter ϕ taking the three different values 0.2, 0.5 and 0.8. But neither the charts nor the schemes for ARL computation are limited to such DAR(1) hidden states, which is illustrated later in Section 5. The respective parameter values of the investigated charts, corresponding to the ARL performances in Figures 3-5, are displayed in Table 1. Here, the target value μ 1 = 3.0225 for the log-LR λ chart results from multiplying λ with 1.55, and the one for the log-LR π chart from the "shifted" marginal distribution π 1 = (0.324, 0.227, 0.449) > . Note that the above shift scenarios used for evaluating the charts' performance include both of these anticipated out-ofcontrol situations. For the c chart design, we always ended up with the same value u = 9 for the control limit because of discreteness. As a consequence, the ARL 0 values in Table 1 slightly increase with increasing dependence parameter ϕ, a F I G U R E 3 ARL performances against mean shifts δ = μ − μ 0 , μ 0 = 1.95, for two different degrees of dependence: ϕ = 0.2 (left side) and ϕ = 0.8 (right side). Mean shift δ caused by uniform shift in λ in upper panel, and by change in π in lower panel well-known phenomenon for Shewhart-type charts. 23,25 To ensure a fair comparison, also, the different CUSUM charts have been designed to (roughly) meet these c charts' ARL 0 values.
The two upper graphs in Figure 3 show the ARL curves where all state-dependent means are shifted uniformly. For the low dependence level 0.2 (left part), all charts show almost the same performance. At least for smaller shifts, say δ ≤ 2, there are some minor differences in the charts' performances, such as the log-LR λ chart tends to exhibit the best performance. This becomes clearer when ϕ rises up to 0.8 (right part), now the log-LR λ chart performs clearly best concerning smaller shifts. Actually, this superior performance is not surprising as the log-LR λ chart is specifically designed for such a uniform shift in λ. Moreover, it can be concluded that the log-LR π chart is least suited for this scenario as its ARL curve locates above the other ones throughout all shift sizes. This can be explained by the fact that the true out-ofcontrol situation strongly differs from the anticipated one.
The opposite situation is considered in the two lower graphs in Figure 3, where the mean shifts are caused by changes within the states' PMF π. By an analogous reasoning as before, it is not surprising that the log-LR π chart performs very well for both dependence levels. However, it is interesting to observe that the ordinary CUSUM chart has F I G U R E 4 ARL performances against mean shifts δ = μ − μ 0 , μ 0 = 1.95, for two different degrees of dependence: ϕ = 0.2 (left side) and ϕ = 0.8 (right side). Mean shift δ caused by shift in λ 0 (upper panel), by shift in λ 1 (middle panel) and by shift in λ 2 (lower panel) nearly the same sensitivity towards small shifts, and even a slightly better sensitivity towards medium-sized shifts for ϕ = 0.8 (note that the largest possible shift size in this scenario is δ = λ 2 − μ 0 = 3.05). In contrast, the c chart exhibits a rather bad ARL performance, irrespective of ϕ. The log-LR λ chart is even worse for ϕ = 0.8 (even a nonmonotonic behavior in δ). So together with our analyses of the upper panel in Figure 3, it becomes clear already at this point that the log-LR charts are very sensitive regarding deviations from the anticipated out-of-control scenario. Figure 4 shows the remaining shift scenarios of our investigations (again left column ϕ = 0.2 and right column ϕ = 0.8), that is, only one of the three state-dependent means λ q is shifted such that the increase in the (overall) mean μ is still the same as in the two preceding scenarios. Starting with shifts of the lowest state-dependent mean λ 0 , we can observe that the c chart performs very poorly concerning smaller shifts for both ϕ values. Keeping in mind the c chart's simple structure (i. e., the observed counts X t are plotted against the upper limit u), it seems plausible that its ARL does not significantly decline in this case: If the hidden state equals 1 or 2, counts are emitted with the same state-dependent mean as before, and if the hidden state equals 0, the corresponding mean λ 0 is somewhat larger but still far away from the control limit u. For large shifts, however, λ 0 exceeds u and the c chart consequently shows its well-known advantage over other control charts. A similar pattern appears in case of λ 1 -shifts, though with an improved performance for small shifts, since already the in-control value of λ 1 is closer to u than λ 0 . Finally, for the λ 2 -shifts (with the in-control value of λ 2 being closest to u), the c chart has even one of the best performances throughout all shift sizes. This time, its simple structure has a positive effect on the ARL, even for smaller shifts, because λ 2 moves even closer to u, and exceeds the control limit much earlier (in terms of δ-shifts) than λ 0 and λ 1 do in the previous shift scenarios.
The most striking performances in Figure 4 are caused by the log-LR λ chart. While the chart performs pretty well for ϕ = 0.2, its ARL curves become nonmonotonic for ϕ = 0.8, note the (local) maxima in both curves for λ 0 -and for λ 1shifts. To understand this behavior, it is important to know that these local maxima are attained if λ 0 and λ 1 , respectively, are approximately equal to λ 2 = 5. In the boundary cases λ 0 = λ 2 and λ 1 = λ 2 , respectively, the three-state HMM actually reduces to a two-state HMM with increased probability mass in the uppermost stateso these λ q -shifts are equivalent to the pure π-changes towards (0, 0.35, 0.65) and (0.5, 0, 0.5), respectively. (Although we omitted to print these numbers in the manuscript, it is worth mentioning that for ϕ = 0.9, a further local maximum appears in the ARL curve of the λ 0 -shifts, located where λ 0 approximately reaches the value of λ 1 .) Generally, the out-of-control states are not uniquely characterized in the sense that, for example, the sole shift λ 0 ↦ 3 cannot be distinguished from the joint shift λ 0 ↦ 2, λ 1 ↦ 3, π 0 ↦ 0.35, π 1 ↦ 0.5. So corroborating the results in the lower panel of Figure 3, the log-LR λ chart seems to be severely affected by deviations from its anticipated out-of-control scenario (though becoming apparent only if the dependence level is sufficiently large). For the λ 2 -shifts, in contrast, we have a strictly monotone ARL curve as λ 2 is the largest among the state-dependent means, that is, the ordering among the means does not change by such shifts. To sum up, the log-LR λ chart exhibits serious weaknesses when a shifted λ q approaches another state-dependent mean λ r , r ≠ q, which is more pronounced for larger ϕ.
Also the log-LR π chart shows a clear pattern: the ARL converges to some limiting level for increasing shifts, which is lowest in the upper panel, somewhat larger in the middle panel and largest in the lower panel of Figure 4. Again, this behavior is considerably more pronounced for ϕ = 0.8. Our explanation is as follows: In the upper panel, once λ 0 has exceeded λ 2 , the states' PMF π (corresponding to the ordered state-dependent means λ q ) is (0.35, 0.15, 0.5), which is pretty close to the anticipated out-of-control π given by (0.324, 0.227, 0.449). Analogously, in the middle panel, we end up with an ordered π of (0.5, 0.15, 0.35) again close to (0.324, 0.227, 0.449). For the λ 2 -shifts in the lower panel, in contrast, the order within π remains (0.5, 0.35, 0.15) throughout all shift sizes and thus exhibits the least fit with (0.324, 0.227, 0.449).
Let us complete the analysis of the single λ q -shifts by briefly discussing the ordinary CUSUM chart's performance. Obviously, the CUSUM chart yields an almost equally good performance in all panels of Figure 4 (this impression is also confirmed for ϕ = 0.5 in Figure 5 below), that is, the CUSUM chart's ARL performance exhibits the best overall robustness against both the actual autocorrelation level and the specific shift scenario. Figure 5 displays the resulting performances for the moderate dependence level ϕ = 0.5, now sorted by the control charts instead of the different shift scenarios. As already indicated above, the CUSUM chart shows the most homogeneous performance regarding the different types of shifts across all investigated charts. Any other chart turns out to be very sensitive to the actual shift type: The c chart has strongly varying ARL values for small shifts, the log-LR π chart performs very badly for λ 2 -shifts, and the log-LR λ chart has nonmonotonic ARL curves for the λ 0 -, λ 1 -and π-shifts. So if we cannot be sure about the expected type of out-of-control scenario, only the ordinary CUSUM chart appears to be a reliable choice.

| APPLICATION TO SALES COUNTS DATA
To illustrate the application of the c and CUSUM chart as well as the log-LR approach, we continue the example about sales counts as introduced in Section 2. There, it was argued that a stationary three-state Poisson HMM with transition matrix This model is considered as the in-control model. For Phase II monitoring, we simulate a time series ðx 1 ,q 1 Þ, …, ðx 100 ,q 100 Þ of length 100, which is initialized by the last hidden state of our Phase I data according to the global decoding shown in Figure 2, that is, byq 0 = 0. The first 20 observations are simulated according to the in-control model. But then, there is a change point (marked by the vertical dotted lines in Figure 6), at which the two lower F I G U R E 5 ARL performances against mean shifts δ = μ − μ 0 , μ 0 = 1.95, of c, CUSUM, log-LR λ and log-LR π chart for ϕ = 0.5 state-dependent means increase from λ 0 = (3.74, 8.44, 14.93) to λ 1 = (6, 12, 14.93), whereas all other parameters remain the same. Consequently, the model mean (variance) rises from 5.42 (14.72) to 7.84 (17.01). The remaining 80 observations are simulated according to this out-of-control model. The resulting Phase II data, that is, 100 weekly sales counts together with the true hidden states (or their means, respectively), are plotted in the first graph of Figure 6. Note that there are three possible hidden states before and after the dotted line, but their actual mean levels (the grey squares) change from λ 0 to λ 1 . In terms of a practical meaning of this shift, the out-of-control model assumes the same mean number of sales in the high-demand state, whereas in weeks with a low or medium demand, the supermarket is confronted with an increased sales activity.
As the next step, we have to select and design the control charts to be used for Phase II monitoring. As before, we consider the c, the CUSUM and the log-LR λ chart for process monitoring, but the log-LR π chart shall not be applied in this section. The reason for this is as follows: For our performance analyses in Section 4, we assumed a DAR(1) structure for the hidden states' MC, because this model allowed us to separate the marginal PMF π from the extent of serial dependence as expressed by the parameter ϕ. The hidden MC of our sales counts data, however, turned out to be more complex than just DAR(1), and the relation between π and A is not obvious. Thus, it would be difficult to change (A, π) for receiving an anticipated μ 1 . So we concentrate on c, CUSUM and log-LR λ chart this time, and the charts' design parameters are chosen as • u = 20 for the c chart, yielding (the in-control) ARL 0 245.35, • (k, h) = (7, 47) for the CUSUM chart, yielding ARL 0 244.37, and • (μ 1 , ℓh) = (8.4057, 3.57) for the log-LR λ chart, yielding ARL 0 245.15.
These charts are now applied to the Phase II data, see Figure 6 (the displayed results can be replicated by applying the corresponding R code in the supplemental materials). After the shift, the ARLs reduce to 152.38 (c chart), 53.37 (CUSUM chart) and 23.84 (log-LR λ chart). In line with these values, the log-LR λ chart is the first chart triggering an alarm for the considered simulation, that is at t = 36 (delay 15). It is followed by the CUSUM chart, whose control limit is first exceeded at t = 55 (delay 34). Before the change point though, the statistics of both charts are always almost 0 such that no false alarms are triggered, see the corresponding graphs in Figure 6. In contrast, the c chart does not signal any alarm within the considered time interval.
The bad performance of the c chart is mainly due to the shift-type that does not affect the value of λ 2 . We thus have a situation similar to what has been previously discussed in Section 4 regarding the two upper panels of Figure 4, where the value of λ 2 remains unchanged likewise. In the present case, the chosen shift-type only concerns those statedependent means, which are farthest below the control limit u and which additionally even remain below the value of λ 2 after the shift. The deficiency of the c chart is also confirmed by the large out-of-control ARL of 152.38. At the same time, it seems reasonable that the log-LR λ chart outperforms both other charts, as its anticipated out-of-control scenario, 1.55 � λ 0 = (5.797, 13.082, 23.1415), is quite close to λ 1 . But also the simple CUSUM chart does quite well for these data, which confirms our conclusion in Section 4 that the CUSUM appears to be a good compromise: It is well suited for detecting mean shifts independent of the actual shift structure.

| ON STATE-DEPENDENT CONTROL CHARTS
The control charts discussed so far consider the HMM for chart design, but the HMM's structure does not become visible to the operator during process monitoring. However, one may also think about a state-dependent control chart, where the action taken at time t tries to make use of the current state q t . For example, one might define a statedependent c chart, where still the observed counts x t are plotted on the chart, but the control limit u q t varies according to the underlying state q t . A related approach was proposed by Alshraideh and Runger, 5 who define a residuals-based Shewhart chart with 3-σ limits. Here, the residual r t at time t considers the previous state q t − 1 in the following way: In both cases, however, the true hidden states q t cannot be observed, so we have to rely on a (locally) decoded stateq t . This is possible by using (3.2): the most probable state at time t, given the observed counts x 1 , … , x t , equalsq t = argmax q v t,q , see Appendix A. So the state-dependent c chart actually plots x t against uq t , and the residuals chart plotsr t := x t − P d Q j = 0 λ j � a jjq t − 1 . Because we do not know the true data-generating HMM, the local decoding has to be based on the supposed incontrol model, whereas the generated data might be out-of-control. This, however, may lead to systematic misidentifications of the hidden states, that is, we permanently haveq t ≠ q t , which, in turn, deteriorates the charts' out-of-control performance. This can be seen in Figure 7, where we applied both types of investigated state-dependent charts to the Phase II data from Section 5. Comparing the decoded states with the true ones, see the table at the bottom of Figure 7, we observe that during the out-of-control period, the decoded states are often to high. This, however, causes the state-dependent c chart to choose the "wrong" control limit, and the residuals are not sufficiently increased in value. In simulations, we also observed that the ARL performance often shows strong breaks of monotonicity similar to the performance of the above log-LR λ chart, but much more pronounced at the corresponding shift sizes.
In view of the previous results, it appears to be more promising to develop such control charts, which make use of a possible misidentification of the hidden states (in order to detect the out-of-control scenario), but do not suffer from it. A possible approach could be to run multiple (two-sided) Bernoulli CUSUM control charts as introduced by Reynolds and Stoumbos 29 in parallel, to directly monitor the hidden states. For example, we may use one chart for each type of hidden state, see the analogous application by Ryan et al 30 for the monitoring of a categorical process (but Bernoulli CUSUM charts for certain combinations of states might also turn out to be useful). In case of a three-state HMM, this means that three two-sided Bernoulli CUSUM charts with individual reference values and control limits constitute one large scheme (thus, 12 design parameters-in contrast to just 2 parameters for the ordinary CUSUM chart from Section 4). For each q 2 Q , the corresponding qth Bernoulli CUSUM accumulates the binary indicators 1ðq t ≠qÞ, where 1ð�Þ denotes the indicator function. Thus, this chart has the potential to detect any type of out-of-control scenario that effects the identification of the hidden states. Figure 8 presents some first results from a simulation experiment (again with 10 6 replications). These results, however, do not show a clear advantage of the joint Bernoulli CUSUM. Its ARL graphs are close to the ones of the log-LR π chart (which is the natural counterpart among our HMM CUSUM schemes), but with the latter having less design parameters. In the case of λ 2 -shifts, the joint Bernoulli CUSUM chart performs considerably worse than the log-LR π chart, whereas it seems to be more sensitive towards λ 0 -shifts.
Altogether, it appears worth to be analyzed in a future research if the design of the joint Bernoulli CUSUM chart can be optimized to perform superiorly in further out-of-control scenarios.

| CONCLUSIONS
In this work, we considered the task of monitoring counts being generated by a Poisson HMM. Besides the wellestablished c and CUSUM chart, we also proposed types of log-LR CUSUM charts for this purpose. In comparing their ARL performance with regard to diverse out-of-control scenarios, it turned out that both the c chart and the log-LR charts may show a rather bad ARL performance, depending on the actual type of mean shift. The log-LR charts can only be recommended if one can be sure that the type of out-of-control scenario to be expected agrees well with the anticipated one. The ordinary CUSUM chart, in contrast, showed nearly the same ARL performance for any scenario and any dependence level. Therefore, it can generally be recommended as a reliable chart for monitoring mean shifts in a Poisson HMM.
Although the ordinary CUSUM chart turned out to be a good "all-purpose chart", it is still a relevant question to ask if the log-LR approach could be modified to make it an attractive alternative for practice. A solution could be to consider a GLR chart instead, also see previous studies, 7,31,32 where the fixed target value θ 1 in (3.3) is omitted by considering the statistic G t = max τ, θ ln PðX t , …, X τ + 1 jX τ , …, X 1 ; θÞ PðX t , …, X τ + 1 jX τ , …, X 1 ; θ 0 Þ at time t. This statistic, however, does not only require to maximize in τ (i. e., to find the most probable position of the change point), but for each value of τ 2 {1, … , t}, a maximization also has to be done in θ (i. e., altogether t maximum likelihood estimations). Because ML estimation is computationally very demanding for HMMs, this does not seem to be feasible in practice, even if using a moving-window technique like in Wang and Reynolds. 31 Therefore, it would be an interesting task for future research to find out if the computational burden could be substantially reduced by employing appropriate constraints on the out-of-control parameter vector θ, for example, by keeping the parameters related to ðQ t Þ N fixed but by jointly maximizing about all state-dependent means λ q . Certainly, it then has to be checked if such a constrained GLR chart does show a more stable ARL behavior. Another relevant topic for future research was already presented in Section 6, namely, trying to find a successful way of state-based process monitoring. Although the state-dependent Shewhart charts suffer from the problem of decoding the hidden states, the idea of a joint Bernoulli CUSUM chart for state monitoring appears to be quite promising. Finally, the presented performance analyses rely on specified model parameters; future research should also investigate the effect of parameter estimation on the charts' performance. F I G U R E 8 ARL performances against mean shifts δ = μ − μ 0 , μ 0 = 1.95, of log-LR π and joint trivariate Bernoulli CUSUM chart for ϕ = 0.8 respectively, for all q 2 Q, and they can be computed recursively from α 1 = Pðx 1 Þπ, α t = Pðx t ÞAα t − 1 ; β T = 1, β > t = β > t + 1 Pðx t + 1 ÞA: This, in turn, enables a recursive computation of the likelihood function, because LðθÞ = 1 > α T holds. It was also used to derive the recursive scheme (3.2), which is computationally more stable because of the scaling v t = α t /(w 1 … w t − 1 ) with w t = 1 > v t .
For an h-step-ahead forecasting of the HMM's observations, we compute the forecast distribution and take the mode or median thereof as the point forecast value. Analogously, we can predict future hidden states from where e q is the qth unit vector. Concerning a decoding of the already realized (but invisible) hidden states, we distinguish between the "local decoding" of a single hidden state, and the "global decoding" of all hidden states. To decode q t with 1 ≤ t ≤ T, given the observations x 1 , … , x T , we computê q t : = arg max q PðQ t = qjx T , …, x 1 Þ = arg max q α t,q � β t,q 1 > α T : For an online local decoding, we set T = t in this scheme (or h = 0 in the above forecast distribution for Q t + h ), leading toq t = argmax q v t,q . A global decoding can be done by using the "Viterbi algorithm" to maximize P(Q T = q T , … , Q 1 = q 1 j x T , … , x 1 ). For each q 2 Q, we define the probabilities m 1,q := PðQ 1 = q, X 1 = x 1 Þ = pðx 1 jqÞπ q , m t + 1,q := max q 1 ,…,q t PðQ t + 1 = q, Q t = q t , …, Q 1 = q 1 , X t + 1 = x t + 1 , …, X 1 = x 1 Þ: These are computed recursively as m t + 1,q = pðx t + 1 jqÞ � max r fm t,r � a qjr g for t≥1: Then the globally decoded states arê q T := arg max q fm T,q g,q t := arg max q fm t,q � aq t + 1 jq g for t = T −1, …,1:

A P P END I X B: IMPLEMENTATION OF MC APPROACH
Let us describe the implementation of the MC approach by Brook and Evans 27 for the c chart and the CUSUM chart applied to a HMM; also see Weiß 33 for general details concerning the MC approach. The approach assumes that the monitored process can be represented as a finite homogeneous MC ðZ t Þ N with state space S = N [ fag, where the set N contains the "no-alarm states" and the set {a} consists of a single "alarm state" "a." We define Z > to be the transition matrix for the states in N , that is, Z > := ðp ijj Þ i,j2N . Furthermore, let I be the identity matrix and 1 the vector of ones. To compute the zero-state ARL, we first have to solve the linear equation (I − Z)μ = 1 in μ. If z denotes the vector of initial probabilities P(Z 1 = j) for j 2 N , then the zero-state ARL is computed as 1 + μ > z.
Note that the choice of z is not unique if computing out-of-control ARLs, see the discussion of Scheme 2.3 in Weiß. 33 For computational simplicity, we followed the approach in Weiß 33 and defined z as the stationary marginal distribution of the actual out-of-control model. We also experimented with different initializations z, but the effect was very small.

ARL computation for c chart
The bivariate process ðX t , Q t Þ N 0 is a MC with transition probabilities P ðX t , Q t Þ = ðx, qÞjðX t − 1 ,Q t − 1 Þ = ðy, rÞ ð Þ = PðX t = x jQ t = qÞPðQ t = qjQ t − 1 = rÞ = pðxjqÞa qjr : The set of "no-alarm states" is given by which is a finite set of size ðu + 1Þðd Q + 1Þ such that the MC approach is applicable.

ARL computation for CUSUM chart
The trivariate process ðX t , Q t ,C t Þ N 0 is an MC with transition probabilities P ðX t , Q t , C t Þ = ðx, q, cÞjðX t − 1 ,Q t − 1 , C t − 1 Þ = ðy, r, dÞ ð Þ = PðC t = cjX t = x,C t − 1 = dÞPðX t = x jQ t = qÞPðQ t = qjQ t − 1 = rÞ where 1ð�Þ denotes the indicator function. For ease of implementation, we choose k, h 2 Q + (with same denominator n), so C t 2 C := f0, 1=n,2=n, …, hg. It can then be concluded with same arguments as in Weiß 24 that the set of "no-alarm states" equals Note that not all potential "no-alarm states" can be attained by (X t , Q t , C t ) given the triple ðX t − 1 , Q t − 1 , C t −1 Þ 2 N . Therefore, using sparse matrix techniques for the transition matrix is an efficient way of implementing the MC approach in order to compute the ARL.