Multivariate control charts to monitor the monthly frequency of vehicle robberies in São Paulo city

Accepted Manuscript Multivariate control charts to monitor the monthly frequency of vehicle robberies in São Paulo city. Adriano B. Moala, Linda L. Ho, Roberto C. Quinino

PII: DOI: Reference:

S2211-6753(16)30173-7 https://doi.org/10.1016/j.spasta.2018.09.002 SPASTA 326

To appear in:

Spatial Statistics

Received date : 2 December 2016 Accepted date : 28 September 2018 Please cite this article as: Moala A.B., et al., Multivariate control charts to monitor the monthly frequency of vehicle robberies in São Paulo city.. Spatial Statistics (2018), https://doi.org/10.1016/j.spasta.2018.09.002 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Multivariate control charts to monitor the monthly frequency of vehicle robberies in Sa˜o Paulo city. Adriano B. Moala, Linda L. Ho University of S˜ ao Paulo, Brazil

Roberto C. Quinino Federal University of Minas Gerais, Brazil

Abstract In this paper, vectors of the deviance residuals (after fitting a STARMA model) are used to build MCUSUM and MEWMA control charts to monitor multivariatespace-time count series. Chart parameters are estimated by simulation to meet a desired in-control average run length and to minimize out-of-control average run length. A complementary simulation study is performed to measure the impact of the omission of the spatial dependencies on the performance of the control charts. Results highlight that false alarms will be signaled much earlier on average. For illustrative purposes, consider the data set of the monthly rates of vehicle robberies registered in 93 police districts located in S˜ao Paulo City (Brazil). Data from January 2001 to December 2013 are used to fit the STARMA model. Chart parameters are searched by simulation and used to draw the control charts to monitor observed vehicle robberies from January 2014 to April 2016. The control charts signal as out-of-control for almost all months of 2014 and the beginning of 2015. An exploratory analysis is used to Preprint submitted to Spatial Statistics

July 18, 2018

identify which districts are responsible for these signals. In the case of omission of spatial dependencies, in the current application, these control charts will give false alarms on average 2.5 (MEWMA) and 7 months (MCUSUM) earlier due to this fault. Keywords: MCUSUM; MEWMA; STARMA; Negative binomial distribution; ARL0 ; Car robbery 02.50.Sk 02.50.Tt 02.60.Pn 02.70.Uu 05.45.Tp 62H11 62M10 62P99 37M10

2

Multivariate control charts to monitor the monthly frequency of vehicle robberies in Sa˜o Paulo city. Adriano B. Moala, Linda L. Ho University of S˜ ao Paulo, Brazil

Roberto C. Quinino Federal University of Minas Gerais, Brazil

1. Introduction Traditionally, control charts have been used to monitor manufacturing processes. Usually at a fixed sampling interval, items are sampled, value(s) of quality characteristic(s) is (are) collected at each sampled unit, a statistic is calculated and if it meets some decision criterion, it is said that the process is in-control and thus production continues (that is, the process is stable and running only under its natural variability). Otherwise it is said that the process is out-of-control; thus, the production is stopped to look for the assignable causes. As the control charts are relatively easy to implement and with the advent of total quality management, their application was adopted in other fields of knowledge rather than only in industrial processes, so as to monitor quality of services, public health, syndrome surveillance, epidemic surveillance, and so on. However, these processes cannot be truly stopped when the control chart signals an out-of-control situation like in manufacturing processes. For Preprint submitted to Spatial Statistics

July 18, 2018

example, when a control chart that monitors the number of hospitalizations due to a disease signals an out-of-control situation, this may be viewed as an outbreak and the public health authorities have to take immediate action to avoid an epidemic level. Other examples in public health can be cited, such as epidemic surveillance for diseases like Zika, yellow fever, dengue, or for deforestation control in the Amazon region. All these problems are considered priorities in Brazil. The usual control charts assume independence among the observations. However, they are not if we consider that the responses are count time series collected in several locations. Further, there is interest in monitoring them to detect changes in the expected value, but the usual control charts cannot be employed without some adjustment. In this paper, let us consider the monthly frequency of vehicle robberies registered in police districts (a total of 93 districts) in S˜ao Paulo City. Vectors of the deviance residuals (after fitting a Spatial-temporal Autoregressive Moving Average (STARMA) model) are used to build Multivariate Cumulative Sum (MCUSUM) and Multivariate Exponentially Weighted Moving Average (MEWMA) control charts to monitor multivariate space-time count series. Chart parameters are searched by simulation to meet a desired incontrol average run length (ARL0 ) and to minimize the out-of-control average run length (ARL1 ). A complementary simulation study is performed to measure the impact of the omission of the spatial dependencies on the performance of the control charts. Results point out that false alarms will be signaled averagely much time before. 4

A STARMA model is fitted for the events occurring from January 2002 to December 2013 using as a response variable the monthly rate, defined as the ratio of the monthly frequency of vehicle robberies by the total of registered vehicles in the same month (in thousands). The residual analysis confirms the goodness of fit of the proposed model. Count data are simulated (using as input the predicted values provided by the proposed STARMA model) to build the run length distributions and thus the control limits of MCUSUM and MEWMA control charts are searched to meet a desired value of ARL0 but also to yield a minimum value of ARL1 (as the performance of these control charts is compared in terms of ARL1 ). With these control limits, control charts are applied to monitor the events observed from January 2014 to April 2016. In the case of omission of spatial dependencies, in the current application, these control charts will give false alarms on average 2.5 (MEWMA) and 7 months (MCUSUM) earlier due to this fault. We emphasize that more flexible models are also considered, for example the STARMAG model (Di Giacinto, 2006), which is capable of capturing the heterogeneity of the characteristics of each location, and the STARIMA model, which is a generalization of a family of space–time model STARMA (with an inclusion of a parameter to indicate the number of differences). In this paper, we evaluate the estimates provided by a non-parametric bootstrap procedure and by the package STARMA in R (which employs a Kalman filter and assumes a plain covariance-variance matrix). The impact of these methods is described in Section 4. 5

This paper is organized as follows: the Introduction section is first; in Section 2, we point out the relevance of monitoring the vehicle robberies and the availability of data set for this aim; in Section 3, the main models and control charts found in the literature to accommodate spatio-temporal dependency are briefly reviewed. The Spatial-Temporal Autoregressive Moving Average (STARMA) model fitted to the rates of vehicle robberies occurring from January 2002 to December 2013 is the subject of Section 4; the determination of parameters of the MCUSUM and MEWMA control charts and the comparison of the performance of these control charts in terms of ARL1 are in Section 5. In Section 6, a simulation study is conducted to measure the impact of the performance on the control charts if the spatial components are neglected. In Section 7, the control charts for the observed rates of car robberies from January 2014 to April 2016 are drawn and discussed. In addition, an exploratory analysis is developed to identify locations where the average robbery rate may have increased. The conclusions and discussion are outlined in Section 8. 2. The problem of vehicle robberies in S˜ ao Paulo city. S˜ao Paulo is a city of twelve million inhabitants (and eight million vehicles) located in the Southeastern region of Brazil and it is the capital of the S˜ao Paulo State. The spatial-temporal dependency of car robberies and also their negative impact on economic activities (a loss of 30 billion dollars in 2004, which is approximately 5% of the Brazilian Gross Natural Product, according to The Institute for Applied Economic Research IPEA-Brazil)is well 6

known. Therefore, monitoring the process by a control chart would better guide public policy actions and decisions, the logistic organization, and the allocation of resources, for example. The occurrence of vehicle robbery in S˜ao Paulo City (Brazil) is registered in 93 police districts and the data set can be requested after filling a form at the site http://www.sic.sp.gov.br/. Each police district acts in a delimited region (see Figure 1). Figure 1 shows the aggregated rate of the robberies registered in the 93 police districts for the months of 2002. Each panel is related to one month of 2002, that is, 200201 is for January, 200202 for February, and so on. Similar figures for other years are also drawn but they will not be shown here. Another important issue for analyzing the problem is that the data processing center can only provide this information monthly. In operational and logistic terms, monthly decision-making is considered viable by the public security department. Therefore, even if there were weekly information, which would enable a faster detection of the increase in vehicle robbery rate, it would not be operationally feasible. A plot of the aggregated rate of vehicle robbery by month/year is shown in Figure 2. Note that the rate decreased in the period 2002-2008 and then remained almost stable until 2011. After that, the rate started to increase until 2014 and again declined. The methodology of extra allocation of police forces in the 93 districts is subjective and based on the managers’ experience and information from the responsible official of each district. Often this raises doubts and inefficiency 7

Rate of Vehicle Robbery in 2002 200201

200202

200203

200204

200205

200206

Rate of Vehicle Robbery 30 200207

200208

200209

200210

200211

200212

20 10 0

Figure 1: Rates of vehicle robbery in S˜ ao Paulo city (Brazil) in 2002.

(observed a posterior), consequently implying a waste of the scarce resources destined for extra police forces. The major problem is that the aggregate data in Figure 2 need to be compared with an expected estimate (under a confidence interval perspective) taking into account that the components of the aggregate data (the 93 districts) are spatio-temporal correlated. In addition, when an aggregate value is greater than the expected, the manager still needs to decide in which districts the problem exists and how serious it is in each one of them. This is not an easy task without the use of an appropriate methodology. Thus, the main objective of the study was to provide more accurate information to the public security managers for better decision-making. The process control allows us to evaluate whether the robbery rate is out of control, considering the spatio-temporal structure, and if 8

Rate of Vehicle Robbery

1200

● ●● ● ●

1000

● ● ● ● ● ● ● ● ●● ● ● ●● ● ● ●

●

● ● ●

800

● ●● ● ●● ● ●

●

● ●●● ● ● ●● ●● ●●

● ●

● ● ● ●● ● ●●● ● ●

600

●● ●

● ● ●● ● ● ● ●●● ●● ●●● ● ● ● ● ● ●● ● ● ●●● ● ●● ● ● ● ●●● ●● ●●● ● ● ●● ● ● ● ● ●

●● ● ● ● ● ● ● ●

●●

● ● ●●●●

●

● ● ●●● ● ● ● ● ● ● ●● ●● ●

●

● ● ●● ●● ● ●

400 2004

2006

2008

2010

2012

●●

● ● ●● ●

●

2002

●

● ●●

2014

● ●● ●

2016

Date

Figure 2: Vehicle robbery rate per one thousand vehicles.

it is, to provide subsidies to the districts that may be contributing most to the problem. 3. Reviewing models and control charts for spatio-temporal processes In this Section, the main models and control charts for dealing with space and/or time processes are briefly reviewed. To accommodate the time autocorrelations, the well-known ARIMA models can be used to fit the data. In the literature, models dealing with only spatial correlations can be found. We can mention, for example, the Spatial Error Model (SEM) proposed by Ord (1975) and the Spatial Lag or Spatial Autoregressive Model (SAR) studied by Ord (1975) and Anselin (1988). An extension of the last model is the High-Order Mixed-Regressive Spatial Autoregressive Model (MRSAR), described in Lee and Liu (2010). 9

The most known models considering the two-type correlation components are Spatio-temporal Autoregressive Moving Average (STARMA) and Space-time autoregressive integrated moving average (STARIMA), studied by Cliff and Ord (1975) and Pfeifer and Deutsch (1980a,b); Deutsch and Pfeifer (1981a). The first is a broader model that is an extension of ARMA models, and the second generalizes the STARMA model including one more parameter related to the number of differences. Panel Spatio-Temporal models are proposed by Anselin (2003) but consider only lag of order one for the temporal component. Progress has been made in introducing spatial heterogeneity and/or temporal dynamics to space–time models. Di Giacinto (2006) propose a generalized STARMA model that generalizes the STARMA model to include spatially varying autoregressive and moving average coefficients (GSTARMA). The generalized STARIMA model of Min et al. (2010) also allows the autoregressive and moving average parameters to vary by location and outperforms the standard STARIMA model in terms of forecasting accuracy. Min and Wynter (2011) discuss the Multivariate Spatial-Temporal Auto-Regressive Moving Average (MSTARMA) model presenting real-time traffic prediction of speed and volume over 5-min intervals for up to one hour in advance. Cheng et al. (2014) proposed the LSTARMA model that captures the autocorrelation locally (heterogeneity) and dynamically (nonstationarity). The control chart for monitoring an autocorrelated production process is the subject of many studies. Alwan (1992) studied the effects of autocorrelation on control chart performance considering a general ARMA (p, q) 10

process. Wardell et al. (1994) derived the run-length distribution of the special cause control charts for correlated processes. Zhang (1998) proposed a control chart for stationary process data instead of the residual chart, as this suggestion demands fitting a model before the design of the control chart. Chao-Wen and Reynolds Jr (1999) suggested the EWMA control chart to monitor the mean of autocorrelated processes. Jiang et al. (2000) presented a control chart named “ARMA” control chart to monitor the ARMA statistics of the original observations. Lin et al. (2012) used genetic algorithms to find the parameters for an economical ARMA control chart. Although many contributions related to control charts for autocorrelated data can be found in the literature, most of them dealt with low lag order time series models, such as AR(1), AR(2), or ARMA(1,1), as these are the most common models found in production processes. Additionally, few applications of control charts to analyze crime rates are listed in the literature, such as Anderson and Diaz (1996); Chamlin (1988), but both used ARIMA models. To monitor a vector of autocorrelated data, some contributions are found in the literature. It is not an exhaustive list but most of them are concerned with the temporal component. For example, Cheng and Thaga (2005) proposed a multivariate max-CUSUM chart that is able to detect simultaneously shifts in both the mean vector and covariance matrix of a multivariate process. Noorossana and Vaghefi (2006) studied the effect of autocorrelation on the performance of the MCUSUM control chart. Jarrett and Pan (2007) also presented vector autoregressive (VAR) control charts for multivariate autocorrelated processes. Bodnar and Schmid (2007) described several CUSUM 11

charts to monitor the mean of a multivariate time series. Some are modified charts and based on residual schemes. Fricker Jr et al. (2008) compared two directionally sensitive MCUSUM and MEWMA procedures with application in biosurveillance. Bodnar and Schmid (2011) studied a CUSUM chart for monitoring the mean of a multivariate Gaussian process focusing on a VARMA(1,1) time series. However, few contributions dealing with a mean vector that simultaneously takes into account spatial and temporal correlations are found. Recently Jiang et al. (2011) proposed a set of multivariate surveillance schemes in a multivariate statistical process control based on likelihood ratio considering also spatial correlations, but the response variables follow a multivariate Gaussian distribution. 4. Fitting a STARMA model In this Section, a STARMA model is fitted for the robbery car rates occurring from January 2002 to December 2013. The stages of the STARMA modeling are identification, estimation, and diagnostic checking. The procedure was applied repeatedly until we obtained an appropriate model. First, in Subsection 4.1 the STARMA model is presented. Following the notation, assumptions, and inputs used to fit the model is presented in Subsection 4.2 and the choice of the final model is in Subsection 4.3. 4.1. STARMA Model An extension of the well-known ARMA model to deal with space-time dependency is the STARMA(p, λp , q, δq ) model proposed by Cliff and Ord

12

(1975) and Pfeifer and Deutsch (1980a,b), which is expressed as Z∗t =

p λk X X k=1 j=0

φk,j Wj Z∗t−k −

q δk X X

θk,j Wj εt−k + εt ,

(1)

k=1 j=0

where • p and q are the lags of autoregressive and moving average components, respectively; • λk is the degree of spatial dependency within the k-th autoregressive lag component; • δk is the degree of spatial dependency within the k-th moving average lag component; • φk,j are the parameters of the autoregressive components; and • θk,j are the parameters of the moving average component. Let Xi,t be the monthly number of vehicle robberies at the i-th location at time t. It is assumed that Xi,t follows a Negative Binomial distribution (as overdispersion is also observed, as in Alencar et al. (2015) in their surveillance Xi,t ∗ monitoring) and Zi,t = × 106 is the respective rate of vehicle robberies, f lt where f lt is the fleet of registered vehicles at time t and i = 1, . . . , N . For ∗ ∗ ∗ N (here N = 93) locations, at instant t, Z∗t = (Z1,t , Z2,t , · · · , ZN,t ). In an

absence of time information this vector is replaced by Z∗ = (Z1∗ , Z2∗ , · · · , ZN∗ ). In this context, Zi∗ indicates the random response at the i-th location.

13

The model assumes that εt = {ε1,t , ..., εN,t } in (1) are normally distributed with µ = 0 and variance-covariance matrix Σ = σ 2 IN , where IN is the N × N identity matrix. The standardized matrix Wj , with dimension N × N , is used to describe the spatial neighborhood relationship of order j among N locations. For an order j, the elements wimj > 0 indicate the neighborhood relationship P strength between i-th and m-th locations with wiij = 0 and N m=1 wimj = 1,

where i = 1, ..., N . For j = 0, the matrix W0 = I. So φk,0 and θk,0 are, respectively, the “pure” temporal components of the autoregressive and the moving average in the STARMA model.

Two particular models are derived from STARMA models: STAR(p, λp )—models with only spatio-temporal autoregressive components, and STMA(q, δq )—models with only spatio-temporal moving average components. 4.2. Definitions, Assumptions, Inputs In the current application, some definitions, assumptions, and inputs are needed. While anisotropy is important in spatio-temporal modelling, the exact locations (that is, the latitude and longitude) of all occurrences of robbery were not available, but only the police district responsible for the locality in which the robbery occurred. Accordingly, as our raw data only show the number of robberies in each of the 93 police districts over the years, we did not use anisotropic aspects in the model, but the concept of neighboring police districts. This means that the distance calculated corresponds to the borders of the police districts, not to the place of robbery. 14

So let i and j be two sites, located at (Pi,xi , Pi,yi ) and (Pj,xj , Pj,yj ), respectively (expressed in degrees). The Haversine distance takes into account the Earth curvature radius (R) in its calculation, thus the Haversine distance d(i, j), in km, between the site i and j is expressed as s         P π P π d π d π d(i, j) j,y i,y x y j i −1 = sin + cos cos sin2 , sin2 2R 180 180 180 180 (2) where R ≈ 6, 372.79548 km, dy = (Pj,yj − Pi,yi )/2, and dx = (Pj,xj − Pi,xi )/2. Consider two police districts A and B, delimited by polygons defined, respectively, by the points {A1 , A2 , . . . , AI } and {B1 , B2 , . . . , BJ }. Let d(Ai , Bj ) be the Haversine distance of points Ai and Bj as expressed in (2). In this study, the distance between districts A and B is defined as δAB = min{d(A1 , B1 ), . . . , d(A1 , BJ ), d(A2 , B1 ), . . . , d(A2 , BJ ), . . . , d(AI , BJ )}. Let (a1 , a2 , . . . , ak ) be a k-tuple of distances with a1 < a2 < . . . < ak . The police districts are neighbors of order 1 if the minimum distance δ is ≤ a1 km, neighbors of order 2 if their minimum distance δ ∈ ]a1 , a2 ], . . . , and finally are neighbors of order k if their minimum distance δ ∈ ]ak−1 , ak ]. Moran’s index has been used to estimate which boundary distances would enhance the spatial autocorrelation as it uses the district grouping described by the data and is also easily interpreted. So, let It (ai ), be the Moran Index (Moran, 1950) related to the neighbors of order i at time t. An overall index related to the configuration (a1 , a2 , . . . ,

15

ak ) at time t is defined as It (a1 , a2 , . . . , ak ) =

k X

It (ai )

i=1



  and let F I(a1 , a2 , . . . , ak ) = median It (a1 , a2 , . . . , ak ), t = 1, . . . , T . The selected configuration (a01 , a02 , . . . , a0k ) is searched such that

  (a01 , a02 , . . . , a0k ) = arg(a1 ,a2 ,...,ak ) M ax F [I(a1 , a2 , . . . , ak )] . In this paper, (a01 , a02 , a03 ) = (0.5, 3, 6); as for the higher order level matrix Wj , j ≥ 4 does not aggregate information. Thus, the following order neighbor matrices Wj , j = 1, 2, 3 are considered: • Order 1: Matrix W1 —for the police districts whose minimum distance δ is ≤ 0.5 km; • Order 2: Matrix W2 —for the police districts whose minimum distance δ ∈ ]0.5km; 3km]; • Order 3: Matrix W3 —for the police districts whose minimum distance δ ∈ ]3km; 6km] The model (1) assumes that εt = {ε1,t , ..., εN,t } are normally distributed with µ = 0 and variance-covariance matrix Σ = σ 2 IN , where IN is the N ×N identity matrix. In the initial checking stage of the STARMA model, such assumptions were not considered reasonable, and we opted to evaluate possible transformations in the response variable that would allow considering the normality assumption more reasonable (as final diagnostic checking pro16

cedures after fitted the model are in the next steps). After evaluating several possibilities, we observed that those that yield a normal variable presented the best results in the sense of getting closer observed and estimated values. One transformation is the square root of the random response (Figure 3 DR panel (b) ) and the another is based on the deviance residuals Zi,t (Figure

3 - panel (c) ) expressed as DR Zi,t

=

∗ sign(Zi,t

q 2 , − µi ) gi,t

2 where gi,t is calculated as follows:    ∗ 2γi ln (1 + µi /γi ) =0 if Zi,t 2 gi,t =    1+Z ∗ /γi  ∗ ∗ ∗ ∗ 2Zi,t /γi ) ln 1+µi,ti /γi /µi ) − 2γi (1 + Zi,t ln (Zi,t > 0. if Zi,t

(3)

(4)

∗

µi in (3) and (4) is replaced by the sample mean Z i for any t, and parameter ∗ ) = γi πi /(1−πi )2 , with πi = µi /(µi +γi ). γi is estimated by satisfying V ar(Zi,t

These equalities are due to the assumption that the original random count variable Xi,t follows a Negative Binomial Distribution (Alencar et al., 2015). Figure 3 shows three histograms. The left one (Panel a) is of the nontransformed random response (showing asymmetry). The middle one (Panel b) is built with the standardized square root values and the right one (Panel c) DR is the histogram of the deviance residuals Zi,t . It is important to emphasize DR that the transformation (3) yields standardized values Zi,t , with the sample DR

mean Z i , i = 1, . . . , 93 closer to zero, their standard deviations closer to one, and most of values in the range [-3;+3] for every district.

17

(a)

(b)

2500

(c)

0.4

0.4

0.3

0.3

1000

Density

1500

Density

Quantity

2000

0.2

0.1

500

0

0.1

0.0 0

20

40

Rate of Vehicle Robbery

60

0.2

0.0 −6

−3

0

3

Square Root of Rate

6

−6

−3

0

3

6

Deviance Residual of Rate

Figure 3: Histogram of the rate of vehicle robbery and its transformations.

Due to asymmetry of the original random response, the STARMA model is fitted using as response the deviance residual expressed in (3) as it provides a higher p-value when Anderson-Darling normality test is applied on aggregate data. The square root transformation is also considered; however, the observed and the fitted values revealed larger discrepancies, so this option is discarded. Hereafter the notations Zi,t and Zt = (Z1,t , Z2,t , · · · , ZN,t ) are kept in the text for the STARMA models but they are all related to the transformation DR Zi,t expressed in (3).

The choice of the lags p and q, respectively, for the autoregressive and the moving average components in ARM A(p, q) models for time series is made by examining the behavior of the autocorrelation function (ACF) and the partial autocorrelation function (PACF). Similarly, the usage of the space time autocorrelation function(STACF) and the space time partial autocorrelation 18

function (STPACF) is primordial in determining the lags in the STARMA model. Additionally, STACF and STPACF will be repeatedly calculated and analyzed in each modeling process improvement until obtaining independent random residuals (free of the correlations in spatio-time components) analogous to the residual diagnostics in ARMA models. DR STACF and STPACF (of the Zi,t series) for the spatial orders 0-3 are

obtained (but not shown here). The initial autocorrelations’ spatio-temporal decay is slow in all four orders but the STPACF shows a strong decay in temporal lags. The partial autocorrelation between the temporal lag one and the spatial order zero is high as is also the partial autocorrelation between the temporal lag one and the spatial order one. 4.3. The choice of the final model When the STACF presents a slow decay, Pfeifer and Deutsch (1980a,b) recommend to start using a simple model like a STAR (1,1), expressed as Zt = φ1,0 W0 Zt−1 + φ1,1 W1 Zt−1 = φ1,0 Zt−1 + φ1,1 W1 Zt−1 .

(5)

STACF and STPACF of the residuals of the model (5) are analyzed. The initial model (5) is improved by introducing new lag components according to the behavior of STACF and STPACF. This procedure is repeated until the observed autocorrelations and the partial autocorrelations of the residuals of the proposed model are closer to null. After some iterations, the final fitted 19

Spatial Order 0

Spatial Order 1

Spatial Order 2

Spatial Order 3

1.0

0.5

0.0

Autocorrelation

−0.5

−1.0 1.0

0.5

0.0

−0.5

−1.0 0

5

10

15

20

0

5

10

15

20

Temporal Lag

Figure 4: STACF of the residuals of final model STARMA(12;2;1;0).

model is STARMA (12; 2; 1; 0), expressed in (6). Note that this model is an incomplete STARMA (12; 2; 1; 0), as the coefficients associated to some lags are statistically not significant (for example the lags (t − 3), (t − 4), etc.), so they are absent. The final model expressed in (6) is chosen as the STACF and STPACF of their residuals are very close to zero (See Figures 4 and 5). Zt = φ1,0 Zt−1 + φ12,0 Zt−12 + φ1,1 W1 Zt−1 + φ2,1 W1 Zt−2 − θ1,0 εt−1 + εt (6) The coefficient estimates of model (6), their respective standard errors (S. E.) and p-values are in Table 1. In this model, we have two types of elements: one related only to time components (as the “pure” autoregressive terms of lag one and twelve, “pure” moving average terms of lag one) and another one related to spatio-temporal components (spatial of order one in the temporal lags of order one and two).

20

Spatial Order 0

Spatial Order 1

Spatial Order 2

Spatial Order 3

1.0

0.5

Partial Autocorrelation

0.0

−0.5

−1.0 1.0

0.5

0.0

−0.5

−1.0 0

5

10

15

20

0

5

10

15

20

Temporal Lag

Figure 5: STPACF of the residuals of final model STARMA(12;2;1;0).

Table 1: Estimates of the parameters of STARMA(12,2,1,0): Final Model.

Elements Parameter Estimate S.E AR Temporal Lag 1 φ1,0 0.80409 0.02172 AR Temporal Lag 12 φ12,0 0.04687 0.00605 MA Temporal Lag 1 θ1,0 0.51406 0.02115 AR Spatio-Temporal Lag 1 φ1,1 0.22849 0.01645 AR Spatio-Temporal Lag 2 φ2,1 -0.10124 0.01848

21

P-value < 0.0001 < 0.0001 < 0.0001 < 0.0001 < 0.0001

ˆ t can be written as a linear combination By expression (7), the vector Z of the past vector rates: ≈ 80% of the (t-1) rate and ≈ 5% of the (t-12) rate. Related to the spatial elements, only neighbor locations of order 1 are included: ≈ 23% of the neighbor location rates at (t-1) and ≈ −10% of the

ˆ t−1 , the prediction neighbor location rates at (t-2). Further, as εˆt−1 = Zt−1 − Z ˆ t is corrected by a percentage ≈ 51% of the residual at (t-1). for Z ˆ t = + 0.80409 Zt−1 + 0.04687 Zt−12 Z + 0.22849W1 Zt−1 − 0.10124W1 Zt−2 − 0.51406 εˆt−1

(7)

In summary, a five-parameter model whose coefficients are easily interpreted explains the variability in most police districts. This final model incorporates the annual seasonality and highlights the strong dependency on the previous month’s rate and also on the previous month neighbors’ rates. Diagnostic checking of this model included an examination of the sample space-time autocorrelation functions of the model residuals (Figures 4 and 5). These autocorrelations were judged to be representative of random residuals. An initial exploratory analysis of the estimated residuals b εt , t = 1, 2, . . . , 132,

confirms that the hypothesis of normality is true for most cases if individual Anderson-Darling normality tests are applied at the 5% level. However, in 10% of cases this hypothesis was rejected and we concluded that this is due to the presence of possible outliers in 13 of the 132 total cases. Moreover, the assumption of independency of the residuals b εt and b εs , s 6= t, arose as the 22

hypothesis of a null correlation was rejected in 9% (8826 pairs) of the possible pairs of b εt and b εs . Furthermore, the variances of estimated residuals b εt , t = 1, 2, . . . , 132, are in the range [0.3; 0.81] and the hypothesis test of equality variances using Bartlett’s test was rejected at the 5% level. Generally, the initial exploratory analysis indicated that problems in the assumptions of the model’s hypotheses were due to a small number of districts. This led us to conduct another analysis. The final diagnostic check concerned the tentative hypotheses H0 : Σ = σ 2 I (errors in each of the N districts are independent with equal variance) and H0 : Σ = I (errors in each of the N districts are independent). The sample covariance matrix of the residuals wa nonspherical and nondiagonal (for details on these tests, see Deutsch and Pfeifer (1981b)), which could compromise the effectiveness of the minimum square estimates and their inferences. Accordingly, the parameters were reestimated using a nonparametric bootstrap procedure as nonassumptions relied on the covariance matrix and the normality distribution of the residuals. This bootstrap procedure consisted of resampling (with replacement) the elements of each vector of the residuals b εt , t = 1, 2, . . . , 132. New values of Zt by (6) were then generated and the parameters reestimated using the minimum square

method. This procedure was repeated 50 thousand times. Details of the bootstrap procedure can be found in Kreiss and Lahiri (2012). The standard error was also estimated using the 50 thousand values yielded by the bootstrap procedure. The bootstrap results confirmed that the analysis on the estimates obtained by the minimum square method was robust in terms of the assumptions of the STARMA model as they are very similar in terms of 23

decision-making. Table 2 presents the estimates obtained from the minimum square method and the non-parametric bootstrap. Additionally, estimates provided by the STARMA package in R are included in the same table. The latter were obtained using the Kalman filter, which assumes a plain covariance-variance matrix. No significant differences were observed in the point estimates and the estimated standard errors. We therefore concluded that the minimum square method was reasonable for decision-making in this case once the hypotheses H0 : Σ = σ 2 I, H0 : Σ = I were rejected and normalizing the component b εt did not significantly change the point estimates of the

parameters and their respective standard errors. Such a situation is not uncommon; for example, Pfeifer and Bodily (1990) employed a STARMA model in a study related to the locations of hotels. The hypothesis H0 : Σ = σ 2 I in their STARMA model was rejected and they employed a plain covariancevariance matrix; in their study, the reestimation similarly had a minor effect on point estimates and their respective standard errors. After defining the model expressed in (7), we show the model to managers with experience with the problem to learn their opinion. First, we do not show the estimates of the coefficients to avoid bias in their evaluation. They reported that Lag-1 is expected to be the most important variable. They pointed out that the police force in each district is similar throughout all the months and the extra force is triggered only in emergencies. Thus, Lag2, Lag-3, Lag-4, and so on (of months in the same year) tend not to add information. Lag-12 is relevant but is less important. It is possible to observe that in the months of July and December we have low robbery frequencies in 24

Table 2: Estimates of the parameters of STARMA(12,2,1,0): Min Square, Bootstrap and Starma Package.

Parameter Method φ1,0 Min. Square Bootstrap Starma AR Temp. Lag 12 φ12,0 Min. Square Bootstrap Starma MA Temp. Lag 1 θ1,0 Min. Square Bootstrap Starma AR Spatio-Temp. Lag 1 φ1,1 Min. Square Bootstrap Starma AR Spatio-Temp. Lag 2 φ2,1 Min. Square Bootstrap Starma AR Temp. Lag 1

Estimate 0.8041 0.8039 0.8067 0.04687 0.0535 0.0566 0.5141 0.5031 0.5471 0.22850 0.2175 0.2195 -0.1012 -0.0884 -0.0889

S.E. 0.0217 0.0189 0.0192 0.00605 0.0070 0.0078 0.0212 0.0235 0.0211 0.0165 0.02231 0.0168 0.0185 0.0263 0.0193

p-value < 10−4 < 10−4 < 10−4 < 10−4 < 10−4 < 10−4 < 10−4 < 10−4 < 10−4 < 10−4 < 10−4 < 10−4 < 10−4 < 10−4 < 10−4

S˜ao Paulo City, as these months are months of school vacations, and parents usually take vacations in the same periods. Thus, the flow of vehicles in the city falls, and so do the dynamics of the families. As a consequence, the opportunities for criminal action are reduced, and such actions occur mostly at the departure or arrival at home in the morning and at night. In addition, Lag-12 may be associated with commemorative dates (Mother’s Day, Valentine’s Day, Children’s Day, etc.). It is important to point out that the spatio-temporal lags W1 Zt are obtained as the average of neighbor locations’ rates as the matrix W1 is standardized and has zeroes for diagonal elements. In the current study, the value of Zt is positively influenced by the events occurring in spatial lag at (t − 1). However, police actions in neighboring locations will yield a negative effect

25

after two months. This indicates a possible migration of vehicle robberies to the neighbor locations. That means some locations may see the vehicle robbery rate increase if an ostensible police action has happened at (t − 2) in the neighbours. Thus, the model suggests that police actions have been locally planned for the short term. This allows the spread of vehicle robberies in the neighbor locations like a virus in an absence of surveillance. When fitting a STARMA model, we are interested in obtaining the pre∗ . They are determined by a search algorithm as the dicted robbery rates Zbi,t

inverse function of (3) does not have an explicit expression. Therefore, for

∗ DR ∗ − − sign(Zbi,t such that Q = |Zbi,t a predicted ZbitDR , we have to find Zbi,t q ∗ 2 | = 0.. This was done by evaluating a long list of values of Zbi,t such µi ) gbi,t

that Q < ϕ, ϕ being a stipulated precision.

Figures plotting observed versus fitted values are not shown here, but the fitted values closely follow the observed ones.

Considering all 93 districts,

the mean absolute error between the fitted and observed values for robbery rates was 1.57 with a standard deviation of 1.61 (both per one million vehicles). This error was considered small by the policymakers responsible for the prevention of vehicle robbery. In other words, greater precision would not have changed the decision-making process. The analysis in this section was conducted using the freeware R 3.2.4 for Windows, the plots through the GGPLOT2 package, and the STARMA model parameters are estimated using the least-squares method through the OPTIM package and STARMA package of R. The code of the programs developed in R and the data analyzed in this 26

study are available as supplementary material for download. The code is in Markdown format to facilitate reading. The data are available in txt format. 5. Setting MCUSUM and MEWMA charts for a STARMA process The literature contains many MCUSUM procedures (such as Woodall and Ncube (1985); Pignatiello and Runger (1990), for example) but, in this paper, we will use the one proposed by Healy (1987); Crosier (1988) and the directional approach presented by Fricker Jr et al. (2008) as we are concerned with detecting locations with increases in the expected criminal rates in the illustrative example. Therefore, recalling the notation Zt = {Z1,t , ..., ZN,t } previously introDR (see expression (3)), at each time t the cumulative duced in Section 4 for Zi,t

statistic St (Crosier, 1988) is obtained as St =

   (St−1 + Zt − µ) 1 −  0, otherwise

k dt



if dt > k;

,

(8)

0

where µ= E(Zt ), k is the solution for k 2 = k Σ−1 k, h i 21 0 dt = (St−1 + Zt − µ) Σ−1 (St−1 + Zt − µ) , and St = (S1,t , S2,t , . . . , SN,t ) where Sj,t



  k = max 0, (Sj,t−1 + Zj,t − µj ) 1 − dt

(9)

for j = 1, . . . , N , in order to include the directional approach presented by Fricker Jr et al. (2008). Starting with S0 = 0, the control chart signals 27



0

−1

whenever Ct = St Σ St

 21

> h.

For MEWMA procedures, in this paper we will consider the extension of Lowry et al. (1992) of EWMA for the multivariate case and the directional approach presented in Joner et al. (2008). Therefore at time t,

Yt = max [0, λ(Zt − µ) + (1 − λ)Yt−1 ]

(10)

can be obtained, choosing a weight λ ∈ ]0; 1[, Yt = (Y1,t , Y2,t , . . . , YN,t ) with Yj,t = max[0; λ(Zj,t − µi ) + (1 − λ)Yj,t−1 ]

(11)

for j = 1, . . . , N . Starting at t = 0 with Y0 = 0, the MEWMA chart signals 0

whenever Et = Yt Σ−1 Yt Yt > b, where

Σ Yt =

λ[1 − (1 − λ)2t ] ΣZ , 2−λ

(12)

ΣZ = σ 2 IN , and IN is the N × N identity matrix.

Assuming model (7) and εt = {ε1,t , ..., εN,t } ∼ N (µ0 = 0; Σ = σ 2 IN ), data are simulated (for the period from January, 2014 to April, 2016) and used to build the empirical run length distribution with the aim to search for the values of parameters k and h for MCUSUM and b for MEWMA in order to meet a desired value of in-control average run length (ARL0 ). In this paper, five thousand simulated run lengths are considered. In production processes, values of ARL0 such as 370 or 500 are the most commonly used if daily data are available, but in the current context, the rates are available monthly. Therefore, an integer number of ARL0 =20 28

months (a year and eight months or ≈ 600 days) is chosen as a reasonable value having an overall type I (αoverall ) error equal to 5%.

Assuming that

all locations have an equal type I error of α, and the type I error of location i is independent of the location j, then αoverall = 1 − (1 − α)93 = 0.05, giving α ≈ 0.055%, which yields an individual ARL0 ≈ 2000 for each location i, i = 1, . . . , N = 93. For determining the control limits, new data under STARMA models are again simulated considering shifts in the vector of the means from in-control µ0 = (0, 0, . . . , 0) to µ1 to obtain ARL1 . In this study, we opt to consider the case where the mean of all locations has shifted equally to some level δ > 0. That is, µ1 = µ0 + δ × σ1, 1 = (1, 1, . . . , 1) (keeping Σ unchanged). The values of δ = (0.25, 0.50, 0.75, 1.00, 1.25, 1.50, 1.75, 2.00) are chosen for illustrative purposes. For example if δ = 1.25, this means that the vector of means has shifted from µ0 = (0, 0, . . . , 0) to µ1 = (1.25, 1.25, . . . , 1.25). The values of ARL, their respective standard deviations (SD), the control limits (b or h), and the constants (λ or k) are given in Table 3. Of course, other levels of shift δi > 0 of interest for location i could be selected to compose the vector µ1 = (δ1 , δ2 , . . . , δN ). By simulation, some pairs of k and h meet the ARL0 for MCUSUM; the pair of values reported in Table 3 is the one that yields the lowest ARL1 . Similarly, the same procedure is adopted to estimate the parameter of MEWMA. Concerning MEWMA, the results denoted “Exact” in the 2nd column (in Table 3) are obtained using expression (12), whereas those denoted “Asymp-

29

Table 3: ARL values for MCUSUM and MEWMA - Model (7)

δ 0 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00 b or h λ or k

MEWMA Exact Asymptotic ARL SD ARL SD 20.10 17.18 20.10 17.18 4.11 0.64 4.67 0.59 2.76 0.44 3.24 0.43 2.01 0.11 2.94 0.23 1.92 0.28 2.06 0.24 1.22 0.42 2.00 < 10−2 1.00 0.04 2.00 < 10−2 −2 1.00 < 10 2.00 < 10−2 1.00 < 10−2 2.00 0.03 213.06 213.06 0.2 0.2

MCUSUM ARL SD 19.88 12.36 3.96 1.03 2.22 0.48 1.57 0.50 1.04 0.19 1.00 < 10−2 1.00 < 10−2 1.00 < 10−2 1.00 < 10−2 2.16 0.70

λ ΣZ , as (1 − λ)2t → 0 2−λ as t → ∞. As expected, “exact” control limits always result in a lower ARL1 . totic” in the 3rd column are obtained with ΣYt =

The ARL01 s are very similar for low levels of shifts (δ < 1.00). Comparing the ARL1 ’s, both control charts (“exact” MEWMA and MCUSUM) have a similar performance for small and large shifts. The largest difference is observed for δ = 1.00, which is almost one month favoring MCUSUM. 6. Impact of neglecting the spatial components on the performance of the control charts In practical terms, ARMA models are well-known and easily fitted for any data set. Thus, neglecting the spatial components may be the most common practice. To measure the impact of this omission (in terms of ARL), new

30

simulations are conducted considering the model Zt = 0.80409Zt−1 + 0.04687Zt−12 − 0.51406εt−1 + εt

(13)

Expression (13) is exactly model (7) omitting the spatial correlations and keeping only the time components. New control limits are determined to reach ARL0 =20 (5 thousand of run lengths are simulated) for MCUSUM and MEWMA for model (13) and then applied to the data simulated with the model expressed in (7). The ARLs are shown in Table 4. The omission of the spatial components will lead to an increase of false alarms, on average 2.5 and 7.0 months, respectively, for MEWMA and MCUSUM before the expected. However, the average time to detect increases of levels of shifts δ is not affected by the use of the equivocal choice of model. Table 4: ARL values for MCUSUM and MEWMA when the spatial terms are omitted

δ 0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00

MEWMA ARL SD 17.68 14.68 3.94 0.62 2.64 0.48 2.00 0.07 1.84 0.37 1.13 0.33 1.00 0.03 1.00 < 10−2 1.00 < 10−2

31

MCUSUM ARL SD 13.09 7.50 3.37 0.79 2.03 0.48 1.35 0.48 1.01 0.11 1.25 < 10−2 1.00 < 10−2 1.00 < 10−2 1.00 < 10−2

7. Application of the MCUSUM and MEWMA control chart in a real data set MCUSUM and MEWMA are applied for the monthly car robberies occurring from January 2014 to April 2016 (see Figure 6). Both charts consider as out-of-control almost all months of 2014 and the first months of 2015. In the current application, MCUSUM signals more cases than MEWMA. It is important to remember that in the multivariate control charts the null hypothesis is that the average robbery rate is stable in all districts, against the alternative hypothesis that at least one average robbery rate has shifted. Further, when the null hypothesis is rejected, there is interest in identifying which components have shifted. To deal with this problem, an exploratory analysis is performed. High values of the statistics (9) and/or (11) for some location j indicate that the observed robbery rate of district j may deviate from its stable mean value and contribute strongly to result in high values of Ct for MCUSUM or Et for MEWMA, respectively, leading to the rejection of the null hypothesis. Figures 7-8 are built using the statistics (9) (for MCUSUM) and (11) (for MEWMA) to identify locations where the average robbery rate may have increased. By Figures 7-8, the districts (in darker colors) located in the periphery areas (such as far East, far South, and far North) have shifted the average robbery rate. Probably, police action is concentrated more in the police districts located downtown causing a migration of robberies to periphery zones.

32

500

500 ●

●

400

●

MCUSUM

●

●

5

● ●

●

●

●

● ●

●

●

●

●

●

● ● ● ●

●

● ● ● ●

●

●

200

● ● ● ● ●

● ●

300

● ●

● ●

0

●

●

●

100

●

● ● ●

●

0

2014−01 2014−07 2015−01 2015−07 2016−01

●

400

●

●

● ●

●

● ●

●

300 ●

● ● ●

200

● ● ●

●

●

●

●

100

● ●

● ● ●

●

●

0

2014−01 2014−07 2015−01 2015−07 2016−01

Date

●

●

● ●

●

MEWMA

10

●

●

MEWMA Asymptotic

●

● ●

2014−01 2014−07 2015−01 2015−07 2016−01

Date

Date

Figure 6: MCUSUM and MEWMA control chart.

8. Conclusions and discussion This paper contributes to applying MCUSUM and MEWMA control charts to deal with a STARMA process in a high dimension case. For illustrative purposes, the monthly rate of robberies of vehicles registered in the 93 police districts located in S˜ao Paulo City (Brazil) from 2002 to 2013 are considered. A five-parameter parsimonious STARMA model is fitted to the data and used for simulation to obtain control limits in order to meet a desired level of ARL0 . In this context, a value of ARL0 = 20 months (one year and 8 months) was chosen as the monthly observations are available. Comparing their performance in terms of ARL1 , both charts have a similar performance; MCUSUM is slightly better for δ ≤ 1.00. For δ = 0.75, 1.00, the differences are two weeks and one month, respectively. For practitioners, it is always recommendable to employ an “exact‘’ control chart rather than “asymptotic” if MEWMA is a choice. Both EWMA and MCUSUM present a similar performance, although the latter has two parameters h and k to estimate and is slightly faster if we think that actions could be taken two weeks or one month in advance (for δ = 0.75, 1.00).

33

201401

201402

201403

201404

201405

201406

201407

201408

201409

201410

201411

201412

201501

201502

201503

201504

201505

201506 Sj,t 2.5 2.0 1.5 1.0 0.5 0.0

201507

201508

201509

201510

201601

201602

201603

201604

201511

201512

Figure 7: Values of Sj,t in (9) versus time and locations

34

201401

201402

201403

201404

201405

201406

201407

201408

201409

201410

201411

201412

201501

201502

201503

201504

201505

201506 Yj,t 1.6 1.2 0.8 0.4 0.0

201507

201508

201509

201510

201601

201602

201603

201604

201511

201512

Figure 8: Values of Yj,t in (11) versus time and locations

35

In practical terms, ARMA models are well-known and supposed to be more widely applied. The impact of spatial components wrongly omitted in the STARMA model will be false alarms on average 2.5 and 7 months earlier than the expected, respectively, for MEWMA and MCUSUM control charts in the numerical example. However, the average time to detect increases of levels of δ is not affected by the equivocal choice of model. When the control chart signals a shift, it is very important to be able to detect which component variables have shifted. In this study, the identification of out-of-control locations was made by a descriptive analysis as the directional approach is included with this aim. Due to the high number of dimensions (almost one hundred locations), an effort to solve this problem is a direction for future studies. The current application of MCUSUM and MEWMA is restricted to the occurrences in S˜ao Paulo City. However, S˜ao Paulo City is neighbored by other cities, so an immediate extension is to consider fitting a similar model to the Metropolitan Region of S˜ao Paulo as it aggregates 49 municipalities and twenty million inhabitants, and then build similar control charts for the monitoring. However, this will cause an increase in the number of locations. For the future directions, one possibility is to consider the MSTARMA model proposed by Min et al. (2010) and build the control charts for monitoring several types of crime simultaneously. Another is the development of a new spatial-temporal model to extend STARMA models for distributions other than the Gaussian distribution, like the recent Generalized autoregressive moving average (GARMA) models developed by Benjamin et al. (2003) 36

(which extended the univariate Gaussian ARMA time series model to a flexible observation-driven model for non-Gaussian time series data). In addition, in this case, statistics other than the vector of deviance residuals can be used to build the control charts. Acknowledgement The authors would like to acknowledge The Brazilian National Council for Scientific and Technological Development (CNPq - 304670/2014-6) for partial financial support. References Alencar, A. P., Lee Ho, L., Albarracin, O. Y. E., 2015. CUSUM control charts to monitor series of Negative Binomial count data. Statistical Methods in Medical Research, http://dx.doi.org/10.1177/0962280215592427. Alwan, L. C., 1992. Effects of autocorrelation on control chart performance. Communications in Statistics-Theory and Methods 21 (4), 1025–1049. Anderson, E. A., Diaz, J., 1996. Using process control chart techniques to analyse crime rates in Houston, Texas. Journal of the Operational research Society, 871–881. Anselin, L., 1988. Spatial Econometrics: Methods and Models. Springer Netherlands. Anselin, L., 2003. Spatial econometrics. A Companion to Theoretical Econometrics, 310–330. 37

Benjamin, M. A., Rigby, R. A., Stasinopoulos, D. M., 2003. Generalized autoregressive moving average models. Journal of the American Statistical Association 98 (461), 214–223. Bodnar, O., Schmid, W., 2007. Surveillance of the mean behavior of multivariate time series. Statistica Neerlandica 61 (4), 383–406. Bodnar, O., Schmid, W., 2011. CUSUM charts for monitoring the mean of a multivariate Gaussian process. Journal of Statistical Planning and Inference 141 (6), 2055–2070. Chamlin, M. B., 1988. Crime and arrests: An autoregressive integrated moving average (ARIMA) approach. Journal of Quantitative Criminology 4 (3), 247–258. Chao-Wen, L., Reynolds Jr, M. R., 1999. EWMA control charts for monitoring the mean of autocorrelated processes. Journal of Quality Technology 31 (2), 166. Cheng, S. W., Thaga, K., 2005. Multivariate max-CUSUM chart. Quality Technology & Quantitative Management 2 (2), 221–235. Cheng, T., Wang, J., Haworth, J., Heydecker, B., Chow, A., 2014. A dynamic spatial weight matrix and localized space–time autoregressive integrated moving average for network modeling. Geographical Analysis 46 (1), 75–97. Cliff, A. D., Ord, J. K., 1975. Space-Time Modelling with an Application to

38

Regional Forecasting. Transactions of the Institute of British Geographers 64, pp. 119–128. Crosier, R. B., 1988. Multivariate Generalizations of Cumulative Sum Quality-Control Schemes. Technometrics 30 (3), 291–303. Deutsch, S. J., Pfeifer, P. E., 1981a. Space-Time ARMA Modeling with Contemporaneously Correlated Innovations. Technometrics 23 (4), pp. 401– 409. Deutsch, S. J., Pfeifer, P. E., 1981b. Space-time arma modeling with contemporaneously correlated innovations. Technometrics 23 (4), 401–409. Di Giacinto, V., 2006. A generalized space-time ARMA model with an application to regional unemployment analysis in Italy. International Regional Science Review 29 (2), 159–198. Fricker Jr, R. D., Knitt, M. C., Hu, C. X., 2008. Comparing directionally sensitive MCUSUM and MEWMA procedures with application to biosurveillance. Quality Engineering 20 (4), 478–494. Healy, J. D., 1987. A Note on Multivariate CUSUM Procedures. Technometrics 29 (4), 409–412. Jarrett, J. E., Pan, X., 2007. The quality control chart for monitoring multivariate autocorrelated processes. Computational Statistics & Data Analysis 51 (8), 3862–3870.

39

Jiang, W., Han, S. W., Tsui, K.-L., Woodall, W. H., 2011. Spatiotemporal surveillance methods in the presence of spatial correlation. Statistics in Medicine 30 (5), 569–583. Jiang, W., Tsui, K.-L., Woodall, W. H., 2000. A new SPC monitoring method: The ARMA chart. Technometrics 42 (4), 399–410. Joner, M. D., Woodall, W. H., Reynolds, M. R., Fricker, R. D., 2008. A one-sided MEWMA chart for health surveillance. Quality, Reliability Engineering International 24 (5), 503–518. Kreiss, J.-P., Lahiri, S. N., 2012. Bootstrap Methods for Time Series. In: Handbook of statistics. Vol. 30. Elsevier, pp. 3–26. Lee, L.-F., Liu, X., 2010. Efficient GMM estimation of high order spatial autoregressive models with autoregressive disturbances. Econometric Theory 26 (01), 187–230. Lin, S.-N., Chou, C.-Y., Wang, S.-L., Liu, H.-R., 2012. Economic design of autoregressive moving average control chart using genetic algorithms. Expert systems with applications 39 (2), 1793–1798. Lowry, C. A., Woodall, W. H., Champ, C. W., Rigdon, S. E., 1992. A Multivariate Exponentially Weighted Moving Average Control Chart. Technometrics 34 (1), 46–53. Min, W., Wynter, L., 2011. Real-time road traffic prediction with spatio-

40

temporal correlations. Transportation Research Part C: Emerging Technologies 19 (4), 606–616. Min, X., Hu, J., Zhang, Z., 2010. Urban traffic network modeling and shortterm traffic flow forecasting based on gstarima model. In: Intelligent Transportation Systems (ITSC), 2010 13th International IEEE Conference on. IEEE, pp. 1535–1540. Moran, P. A. P., 1950. Notes on continuous stochastic phenomena. Biometrika, 17–23. Noorossana, R., Vaghefi, S., 2006. Effect of autocorrelation on performance of the MCUSUM control chart. Quality and Reliability Engineering International 22 (2), 191–197. Ord, K., 1975. Estimation Methods for Models of Spatial Interaction. Journal of the American Statistical Association 70 (349), pp. 120–126. Pfeifer, P. E., Bodily, S. E., 1990. A test of space-time Arma modelling and forecasting of hotel data. Journal of Forecasting 9 (3), 255–272. Pfeifer, P. E., Deutsch, S. J., 1980a. A STARIMA Model-Building Procedure with Application to Description and Regional Forecasting. Transactions of the Institute of British Geographers 5 (3), pp. 330–349. Pfeifer, P. E., Deutsch, S. J., 1980b. A three-stage iterative procedure for space-time modeling. Technometrics 22 (1), 35–47.

41

Pignatiello, J., Runger, G., 1990. Comparison of multivariate Cusum charts. Journal of Quality Technology 22 (3), 173–186. Wardell, D. G., Moskowitz, H., Plante, R. D., 1994. Run-length distributions of special-cause control charts for correlated processes. Technometrics 36 (1), 3–17. Woodall, W. H., Ncube, M. M., 1985. Multivariate CUSUM quality control procedures. Technometrics 27, 285–292. Zhang, N. F., 1998. A statistical control chart for stationary process data. Technometrics 40 (1), 24–38.

42