Sequential monitoring of high‐dimensional time series

In the paper we derive new types of multivariate exponentially weighted moving average (EWMA) control charts which are based on the Euclidean distance and on the distance defined by using the inverse of the diagonal matrix consisting of the variances. The design of the proposed control schemes does not involve the computation of the inverse covariance matrix and, thus, it can be used in the high‐dimensional setting. The distributional properties of the control statistics are obtained and are used in the determination of the new control procedures. Within an extensive simulation study, the new approaches are compared with the multivariate EWMA control charts which are based on the Mahalanobis distance.

The first multivariate control chart was introduced by Hotelling (1947). It is based on the Mahalanobis distance between the vector of observations and the target vector of the characteristics. A control chart based on an MEWMA (multivariate exponentially weighted moving average) recursion was suggested by Lowry et al. (1992). Crosier (1988), Pignatiello and Runger (1990), and Ngai and Zhang (2001) proposed several multivariate CUSUM (cumulative sum) control charts.
The above-mentioned multivariate control procedures were designed for independent observations, which appears to be a rather restrictive assumption in practice (see Alwan & Roberts, 1988). Theodossiou (1993) proposed an extension of the multivariate CUSUM control chart designed to detect changes in the mean of vector autoregressive moving average (VARMA) processes. The extension of the MEWMA control chart to dependent observations was given in Kramer and Schmid (1997), while Schmid (2007, 2011) derived several CUSUM schemes for detecting changes in the mean vector of multivariate time series.
The problem of sequential monitoring of process parameters becomes extremely challenging, when the process dimension is very large. The classical approaches usually cannot be used in such a situation since the control statistics depend on the inverse covariance matrix. The computation of the inverse covariance matrix can be extremely time consuming in practice and also numerical issues can take place. As a result, modifications of the existing approaches are necessary when the parameters of a high-dimensional time series should be monitored.
Several statistical control procedures based on the high-dimensional data have been proposed recently. Wang and Jiang (2009) considered a variable-selection-based multivariate SPC procedure for monitoring the process parameters and fault diagnosis in high dimension. A forward-selection algorithm was utilized to screen out potential out-of-control variables followed by control charts to monitor the suspicious variables. A high-dimensional control chart approach for profile monitoring was suggested by Chen and Nembhard (2011) which is based on the adaptive Neyman test statistic for the coefficients of the discrete Fourier transform of profiles. Li et al. (2014a) proposed a chart that starts monitoring with the second observation regardless of the dimensionality and reduces the average run length (ARL) in detecting early shifts in high-dimensionality measurements. More recently, Wang et al. (2017b) constructed a hybrid control chart in the case of independent multivariate Poisson data which is based on a goodness-of-fit test. A number of challenges appear when these approaches are adopted to the case of dependent data. Although several high-dimensional time series models exist in the literature (see, e.g., Aue et al., 2017aAue et al., , 2017bChudik & Pesaran, 2013;Dias & Kapetanios, 2018;Gupta & Robinson, 2015Hall et al., 2016;Han & Liu, 2013;Han et al., 2015;Kock & Callot, 2015;Lam et al., 2012;Liu et al., 2015;Wang et al., 2017a), it seems that the problem of sequential monitoring of their parameters have not been treated up to now.
Statistical sequential procedures, especially control charts, are widely spread to detect changes in several fields of science. Applications in economics and finance can be found in Theodossiou (1993), Golosnoy and Schmid (2007), Bodnar (2009), Bodnar et al. (2014), Bersimis et al. (2018), while applications in medicine can be found in Shmueli and Burkom (2010), Schiöler and Frisén (2012), Meyer et al. (2017), Barone and Chakhunashvili (2022). Anderson and Thompson (2004) applied multivariate control charts to detect changes in environmental processes, while Ye and Chen (2001) dealt with applications in network monitoring. Despite of the wide spread of sequential control procedures in different areas of research, most of applications are still present in engineering (see, e.g., Ding et al., 2006;Montgomery, 2020).
The rest of the paper is organized as follows. The change-point model is presented in Section 2, while Section 3 provides the results for the multivariate EWMA control chart based on the Mahalanobis distance. The new control schemes designed for high-dimensional time series are proposed in Section 4. Section 5 shows the results of the comparison study. Section 6 summarizes the obtained results, while technical proofs are moved to the appendix (see, Section A).

CHANGE-POINT MODEL
We denote the observed process by {X t }. It is assumed to be a p-dimensional time series. The target process is denoted by {Y t }. The relationship between both processes is given by the change point model where a ≠ 0 and ∈ N ∪ {∞}. Thus both processes may differ by a change in the mean behavior. We say that the observed process is in control if = ∞, else it is said to be out of control. In the following the symbols E ∞ (.), Var ∞ (.), and Cov ∞ (.) will denote the mean, the variance, and the covariance matrix, respectively, computed under the assumption of no change, that is in the in-control state. Throughout the paper it is supposed that {Y t } is a weakly stationary time series with

THE MULTIVARIATE EWMA CHART FOR TIME SERIES
We consider the multivariate EWMA control chart for time series introduced by Kramer and Schmid (1997). It is based on the recursion with Z 0 = where R = diag(r 1 ,…, r p ) with r 1 ,…, r p ∈ (0, 1]. The recursion can also be written as and thus for t = 1, 2,…. Then it holds for p fixed that (cf, Kramer & Schmid, 1997) provided that { (v)} is absolutely convergent. Note that the mean change does not influence the covariance matrix of Z t , that is, Cov The MEWMA control chart is constructed by determining the control statistics at each time point t as the Mahalanobis distance of the vector Z t given by depending on whether the exact or asymptotic covariance matrix of Z t is used (cf, Kramer & Schmid, 1997). It is concluded that there is a change at time t if the control statistics are sufficiently large. Next we will analyze the distributional behavior of T Mah,t and T MahInf,t , in particular in the high-dimensional setting. The following notations will be used throughout the paper.
First, we focus on the control statistic based on the exact covariance matrix.
(c) Suppose that lim p→∞ ,t,p p < ∞. Let t be fixed, then Proof. It holds that Z t − ∼  p (a t− , t,p ). Thus the proofs of parts (a) and (b) follow immediately for (Z t − ) ′ −1 t,p (Z t − ). In order to prove part (c) for (Z t − ) ′ −1 t,p (Z t − ) we apply lemma 3 of Bodnar and Reiß (2016). ▪ In the next lemma the behavior of the control statistic based on the asymptotic covariance matrix for t tending to infinity is analyzed. Since this quantity is easier to determine it is frequently used in applications. However, in the high-dimensional situation it is questionable whether such an approach makes sense. If t is small and p is large the approximation to the exact covariance matrix may be weak.
(c) Suppose that lim p→∞ ,t,p p < ∞. Let t be fixed, then and thus the first part of (a) follows. Now Z t − ∼  p (a t− , t,p ) and the second part of (a) follows with (3.2b.10) of Mathai and Provost (1992). Applying that is a quadratic form its exact distribution can be written as a series expression (cf, Mathai & Provost, 1992). Furthermore, from part (b) it holds that for p fixed The practical calculation of t,p and l;p turns out to be difficult. Explicit formulas can only be obtained for special cases.
Lemma 3. Let {Y t } be a stationary VAR(1) process given by Moreover, let R = r I with r ∈ (0, 1]. Then, it holds that with c t (r, ) = r 2 − r The proof of Lemma 3 is given in Appendix. In the case 1 − r = we interpret the value of c t (r, ) as the limit if converges to 1 − r and the limit coincides with the quantity given in (9).
The MEWMA chart was mainly applied in low-dimensional spaces. However, recently we are faced with situations where the amount of collected data is huge and thus it is of interest to monitor high-dimensional processes. Thus the question arises how good this chart behaves in a high-dimensional context. Following this approach, we might face a number of difficulties in a high-dimensional context. One of them is the computation of the inverse matrix −1 t,p , which in general depends on t and has to be inverted at each time point when it is checked whether the process is still in-control. Even the application of the asymptotic covariance matrix −1 l;p might lead to some difficulties. For example, due to the dimensionality issue, the resulting covariance matrix might be ill-conditioned or its computation might be time consuming.

Introduction of the control schemes
Recently, several authors considered two-sample tests for high dimensions (e.g., Bai & Saranadasa, 1996;Chen & Qin, 2010). They proposed in principle to make use of the Euclidean distance or approximations to this distance in a high-dimensional setting. Motivated by the high-dimensional discriminant analysis (see, e.g., , we consider three types of control statistics which are based on where d;t,p and d;l;p are obtained from t,p and l;p by setting all their nondiagonal elements equal to zero. The process is concluded to be out of control if the value of a control statistic is larger than a control limit c > 0.
The choice of the control limit c plays a central role in the analysis. If c is chosen small the chart will trigger many signals while for a large value of c signals will rarely occur. Frequently, c is chosen such that the in-control ARL is equal to a predetermined value. In engineering, the in-control ARL is often taken to be 500, in finance equal to 60. Essentially, this depends on the frequency of the obtained observations. For a single-step test on the mean vector of the multivariate normal distribution the control limit c is chosen by fixing the probability of the type I error. In this case, the application of the Mahalanobis distance in the construction of the test statistics is advantageous, since the control limit c can be determined as a quantile of a 2 -distribution and thus it is independent of the covariance matrix. However, in the time series context and when a sequential procedure is employed, the advantage of the usage of the Mahalanobis distance in the computation of the control limit becomes redundant, in particular, due to the fact that the joint distribution of the sequentially constructed control statistics involves the autocorrelation function, which in fact depends on the covariance matrix of the error process. In such a situation, the application of other norms to construct control statistics may be considered with the Euclidean distance being a natural alternative to the Mahalanobis distance.
The calculation of the ARL turns out to be quite challenging in the present case. Li et al. (2014b) provide an exhausting review of existent approaches for the computation of the ARL and of the average time to signal (ATS). On the other side, we note that for univariate time series an explicit formula for the ARL is only known for special type of processes as, for example, exchangeable variables (cf, Schmid, 1995). Thus mostly simulations are used to determine the control limit. This approach is of course time consuming since for each parameter constellation and dimension the value must be calculated. Here we are confronted again with the problem of dimensionality.
Practitioners often choose a more simpler procedure. They are working with 3 control limits. In order to apply such a procedure it is first necessary to determine characteristic quantities of the control statistics.

Determination of the control design
Following the proof of Lemma 2 we get the following two results.
Eu;l;i 2 1, Eu;l;i , where Eu;l;1 ,…, Eu;l;p are the eigenvalues of t,p and Eu;l;1 ,…, Eu;l;p are the components of the vector (a) Then t,p and dEu;l;1 ,…, dEu;l;p are the components of the vector dEu; The results of Theorems 1 and 2 are used to determined the designs of several control charts, that are based on the multivariate EWMA recursion of Section 3 with the control statistics computed by using the Euclidean norm and the norm employing the diagonal matrix. Depending on the usage of the exact mean and variance of the control statistics or their asymptotic counterparts, several control schemes are obtained which are the following , , , , , , , , , .
The control statistics T 1,t , T 2,t , T 3,t , T 4,t , and T 5,t are more computationally efficient, since no inversion is used in their computations. On the other side, the control statistics T 6,t , T 7,t , T 8,t , T 9,t , and T 10,t might be more sensitive to detect changes, when the components of X t have different variability. To this end, we note that the control statistics T 5,t and T 10,t are applicable in special situations when the limits can analytically be presented in terms of a finite number of model parameters. This is possible, for example, in the case of VARMA models where all parameter matrices are proportional to the identity matrix. In the simulation study of Section 5 we will consider a more general model for the target process {Y t }, and for this reason these two statistics will not be included. Finally, we note that where d;l,p is obtained from l,p by setting all nondiagonal elements of l,p to zero.

Invariance property
Note that contrary to independent samples the MEWMA chart for stationary processes is in general not directionally invariant. Kramer and Schmid (1997) proved an invariance property of the MEWMA chart for a VAR(1) process in the in-control state. They gave a sufficient condition such that the in-control ARL does not depend on the covariance matrix of the white noise process. If this holds, then the analysis of the MEWMA chart turns out to be much simpler. Moreover, Schmid (2007, 2011) presented conditions under which multivariate CUSUM control charts defined for time series possess the invariance property in the out-of-control state.
In Theorem 3 we present conditions required for the validity of the invariance property in the case of the MEWMA control chart based on the Mahalanobis distance. The proof of the theorem is given in Appendix. (1) the distributions of the run lengths of the MEWMA control charts based on the Mahalanobis depend on a only through a ′ (0)a, that is, if RL denotes the run length of such a chart, then for all a 1 , a 2 with a ′ 1 (0)a 1 = a ′ 2 (0)a 2 it holds that where the symbol P a, (RL > i) denotes the probability computed under model (1) when a shift a takes place at time .
It is important to note that the results of Theorem 3 are derived only for the MEWMA control charts based on the Mahalanobis distance. On the other hand, the MEWMA schemes determined by using the Euclidean distance are not directionally invariant due to the definition of the control statistics used in their construction. Although the invariance property appears to be very helpful in the study of properties of the control charts in the out-of-control situation and the comparison of their performance, in practice, however, it has only a minor impact on the determination of the design of the control charts. Furthermore, the design of the control charts based on the Euclidean distance does not depend in the inverse of the covariance matrix of Z t , which makes the procedures stable in the high-dimensional setting since no inverse of a large-dimensional matrix is required in their construction.

SIMULATION STUDY
The aim of this section is to investigate the performance of the control charts proposed in Section 4 and to compare them with the approaches based on the Mahalanobis distance as described in Section 3.

Results in the in-control state
To investigate the properties of the considered control schemes via simulations, one has to define both the target and the observed process. In this section, we use a VAR(1) process as a target process {Y t }, whose coefficient matrix is assumed to be proportional to the identity matrix, that is = I with = 0.5. For this type of time series model, one can use the results of Lemma 3 to obtain the analytical expressions of t,p and l;p . The covariance matrix of the error process is a correlation matrix with = 0.5. To make the results of the simulation study more flexible, the values of SDs d 1 ,…, d p are drawn randomly from the uniform distribution on the interval [0.5, 2]. Furthermore, it is remarkable that if the matrix D is proportional to the identity matrix, then the control chart T 1,t coincides with T 6,t , T 2,t with T 7,t , T 3,t with T 8,t , T 4,t with T 9,t , and T 5,t with T 10,t . The mean of the VAR(1) process {Y t } is set to be zero vector, that is, = 0. To this end, we put R = r I with r ∈ {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1} and consider two cases of the dimension of the observed process and the target process p to be equal to 50 and 100, respectively.
In the definition of several control statistics in Section 4, the asymptotic mean and the asymptotic variance was used. For the target process as defined in the beginning of the section, we get that the matrices t,p and l;p are proportional following the results derived in Lemma 3. As such, the control charts T 6,t and T 7,t have the same behaviour as T 8,t and T 9,t , respectively. Thus, we will omit presenting the results for T 8,t and T 9,t in this simulation study.
To study the impact of the asymptotic approximation on the performance of the control charts T 1,t , T 2,t , T 3,t , and T 4,t , we compare the exact mean and the exact variance of Z ′ t Z t , t = 1,…, 30, as derived in Theorem 1 with the corresponding asymptotic values. The results are obtained in the case of p = 50 for r ∈ {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1} and are depicted in Table 1. Large differences between the exact and the asymptotic quantities occur for r ≤ 0.7, especially when t = 1. For small values of r considerable differences are also present for larger values of t, up to t = 10. Such a poor performance of the asymptotic mean and variance to provide an accurate approximation of the exact ones might have a strong impact on the control schemes whose definitions involve the asymptotic mean and variance of Z ′ t Z t . This is especially the case with the T 2,t chart where the asymptotic mean and the exact variance is used. This control scheme performs very poorly both in the in-control and out-of-control state and, for that reason it is excluded from the comparison study.
The in-control ARL is used to calibrate the control charts T 1,t , T 3,t , T 4,t , T 6,t , T 7,t , T Mah,t , and T MahInf,t , which are used in the comparison study in Section 5.2. The control limits computed for these control schemes in the case of the considered VAR(1) process are given in Table 2 for r ∈ {0.1, 0.2,…, 1.0} and for both p = 50 and p = 100. The in-control ARL is taken equal to 200. The computation of the control limits is performed by numerical computation where the ARL is found in each iteration based on Monte Carlo simulations with 10 4 independent repetitions. For all considered control charts the control limits increase as r becomes larger. The only exception is the control chart T 4,t , where the control limits first become smaller and start to increase only when r = 0.3. Such behavior of the control limits for the control chart T 4,t is explained by the results presented in Table 1 and the fact that both the asymptotic mean and the asymptotic variance are used in the construction of the control statistics. Since the true values of the mean and of the variance deviate considerably from their asymptotic counterparts when r is small, it leads to the unusual behavior of the computed control limits.
In case r = 1 the MEWMA control schemes become Shewhart control charts. In this case, only the observations at a given time point are used to make a decision about the existence of a change in the model structure, previous observations are not taken into account. The control limits are the largest for r = 1, while they drop monotonically as r decreases. Moreover, we observe that the control limits of T 1,t , T 3,t , T 4,t , T 6,t , and T 7,t charts are considerably smaller than the ones computed for the control schemes based on the Mahalanobis distance. This result follows directly from the observations, that T Mah,t and T MahInf,t are not normalized, while the other five control schemes are centered and normalized by the exact (asymptotic) standard deviations.    TA B L E 2 Control limits computed for T 1,t , T 3,t , T 4,t , T 6,t , T 7,t , T Mah,t , and T MahInf,t control charts in the case of VAR (1)

Results in the out-of-control state
In this section we study the ability of the new MEWMA control charts to detect changes in the location behaviour of high-dimensional VAR(1) processes as defined in the previous section and compare it with the two benchmark approaches that are based on the Mahalanobis distance. The changes in the mean vector are generated according to the change point model (1). Namely, the vector a = (a,…, a, 0,…, 0) ′ with a ∈ {−2, −1.5, −1, −0.5, 0.5, 1, 1.5, 2} is added to the target process {Y t } in the out-of-control state where the number of nonzero elements of a is l ∈ {12, 25, 50} for p = 50 and l ∈ {25, 50,100} for p = 100, respectively. The selection of l corresponds to changes in 25%, 50%, and 100% of the first components of the observed process {X t }.
The maximum expected delay (MED) is used as a performance measure in the out-of-control state. The expected delay is defined as the average delay between the time of a change and the time when the control chart detects a change under the condition that there is no false alarm before the change takes place. The MED takes the maximum value of the expected delays with TA B L E 3 Maximum expected delays for T 1,t , T 3,t , T 4,t , T 6,t , T 7,t , T Mah,t , and T MahInf,t control charts in the case of a VAR(1) process, p = 50, r ∈ {0.1, 0.2,…, 1.0} and a = 0.5  Tables 3-10. For each control chart and out-of-control situation the minimum value of the MED with respect to r is highlighted bold, while the smallest value across all the control schemes and the values of r are presented bold cursive. As expected, changes of small magnitude are detected faster when the small values of r are employed, for changes of moderate size one should prefer to set r ∈ {0.1, 0.2, 0.3}, while large deviations in the mean behaviour of a VAR(1) process are best monitored when r = 1 or when r is close to one. These findings are in agreement with the previous results of Kramer and Schmid (1997) obtained for the multivariate EWMA control charts based on the Mahalanobis distance.
The control charts based on the distance computed with respect to the inverse of the diagonal matrix of variances show the best performance for all considered values of a and the number of elements l whose mean values are shifted. The control schemes T 6,t and T 7,t are followed by the control procedures where the Mahalanobis distance is used in the construction of the control statistics. The worst performance is documented for the control charts based on the Euclidean distance, which are better than the control schemes based on the Mahalanobis distance only when p = 50 and l = 12.

CONCLUSION
Monitoring changes in the parameters of a multivariate time series is a complicated task due to large number of parameters that need to be controlled simultaneously. Although several control schemes exist in statistical literature and are successfully implemented in real-life applications (cf, Bodnar & Schmid, 2005, 2007Kramer & Schmid, 1997), the problem has not been investigated in detail in the high-dimensional setting, that is, when the dimension of the target process can be very large. In a such situation, sequential surveillance of the model parameter becomes extremely challenging, especially in the case when data consist of dependent observations. We contribute to the literature by deriving new types of control procedures based on the MEWMA recursion where in the definition of the control statistics other types of norms are used. In particular, instead of using the Mahalanobis distance, control charts based on the Euclidean distance are suggested as well as control schemes where the norm is computed by employing the inverse of the diagonal matrix consisting of the variances. The two latter approaches possess several advantages with respect to the procedures that are based on the Mahalanobis distance. First, since no inverse of the covariance matrix must be computed, the new approaches can be easily implemented also in the high-dimensional setting. Second, avoiding the computation of the inverse covariance matrix speeds up the computation time considerably. Moreover, the new approaches, especially the ones that uses the inverse of the diagonal matrix in the computation of the control statistics, possess the best performance in the conducted simulation study independently of the model used to generate changes in the mean behavior of a multivariate autoregressive process. When the magnitude of the change is small, we also obtain that the control schemes based on the Euclidean norm outperforms the ones based on the Mahalanobis distance.
One possible direction to extend the results obtained in the paper is to study the impact of the model misspecification on the performance of both types of the control charts based on the Mahalanobis and Euclidean distances. The issues related to the model misspecification appear to be very challenging especially in the high-dimensional time series context. They might be handled by imposing some regularization assumptions on the model parameters. This point is left for future research. then Proof. From z ∼  p ( , ), we get that and hence, where Now we can apply theorem 5 in Dette and Dörnemann (2020) because Moreover, we get with , that (iv) max 1≤i≤p g p (i) 2 Var(u (1 − r) i I |j−i| I (0)(1 − r) j I ⋅ rI (1 − r) i+j |j−i| .
If = 0 then 1 − r ≠ 0 = and we get t,p = r 2 1 − (1 − r) 2t 1 − (1 − r) 2 (0) = (1 − (1 − r) 2t ) which is also obtained by setting = 0 in (7). Results for l;p follow from the respective formulas for t,p by letting t tends to ∞. Now let 1 − r = . Then Note that this result is also obtained from (7) by letting converge to 1 − r. ▪ Proof of Theorem 3. Following Crosier (1988), first we note that the MEWMA control statistics will remain unchanged when a nonsingular transformation is applied to the observed process {X t }. Namely, let M be a nonsingular p × p matrix. Then, where a ≠ 0 and ∈ N ∪ {∞} and̃= E(X t ) = ME ∞ (X t ) = M and Cov(X t+h ,X t ) = MCov ∞ (X t+h , X t )M ′ = M (h)M ′ .
As a result, the multivariate EWMA recursion computed for the process {X t } becomes withZ 0 =̃= M = MZ 0 where we use that R = rI, r ∈ (0, 1]. Moreover, it holds that E(Z t ) = ME(Z t ) and Cov(Z t ) = MCov(Z t )M ′ . Using that that is, they remain the same when a nonsingular transformation is applied to the observation process {X t }.
The rest of the proof follows directly from the proof of theorem 1 in Bodnar and Schmid (2007). ▪