Control charts for high-dimensional time series with estimated in-control parameters

Abstract

In this article, we study the effect of misspecification caused by fitting the target process in the Phase I analysis of the monitoring procedure on the behavior of several types of multivariate exponentially weighted moving average (MEWMA) control charts in the high-dimensional setting. In particular, the classical MEWMA control charts, whose control statistics are based on the exact and asymptotic Mahalanobis distance, are considered together with the novel approaches where the Euclidean distance and the diagonalized Euclidean distance are employed in the construction of control statistics. The high-dimensional distributions of the control statistics are deduced at each time. These results are later used to assess the performance of the considered control charts under misspecification. Both theoretical and empirical findings lead to the conclusion that the control charts based on the Euclidean distance and the diagonalized Euclidean distance are robust to misspecification for moderate dimensions of the data-generating model, whereas they tend to overestimate the in-control average run lengths (ARLs) in the case of larger dimensions. On the other hand, the control schemes based on the Mahalanobis distance are considerably affected by the estimation of the parameters of the target process, and their application results in drastically smaller values of the ARLs, especially when the dimension of the data-generating model is large.

Keywords:

• High-dimensional time series
• MEWMA control chart
• sequential surveillance
• vector autoregressive process

1. INTRODUCTION

Statistical process control (SPC) plays a special role in the monitoring of production processes. The methods of SPC are also widely used in other fields of science, like in engineering, economics, medicine, chemistry, biology, and finance (see, e.g., Frisén 1992; Schipper and Schmid 2001; Sonesson and Bock 2003; Andersson, Bock, and Frisén 2004; Schmid and Tzotchev 2004; Lawson and Kleinman 2005; O. Bodnar 2007, 2009b; Messaoud, Weihs, and Hering 2008; Golosnoy et al. 2011; O. Bodnar and Schmid 2017).

In the setup of the monitoring procedure, the relationship between the observed and target process should be specified. Let $\left\{{\mathbf{Y}}_t\right\}$ denote the p-dimensional target process and let $\left\{{\mathbf{X}}_t\right\}$ be the p-dimensional observed process. Under the target process, we consider a process that fulfills quality requirements, whereas the observed process is the actual process. In the following, the relationship between the target and observed processes is described by the change point model expressed as (1.1) $${\mathbf{X}}_t=\{\begin{array}{lll}{\mathbf{Y}}_t& \text{for}& t<\tau \\ {\mathbf{Y}}_t+\mathbf{a}& \text{for}& t\ge \tau \end{array},\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}t\in \mathbb{Z},$$(1.1) where $\mathbf{a}\ne \mathbf{0}$ and $\tau \in \mathbb{N}\cup \left\{\infty \right\}.$ If $\tau =\infty ,$ then the observed process is called an in-control process. Otherwise, it is called out-of-control process. The symbols ${\mathbb{E}}_{\infty }(.),\mathbb{V}a{r}_{\infty }(.),$ and $\mathbb{C}o{v}_{\infty }(.)$ will denote the mean, variance, and covariance matrix, respectively, computed under the assumption of the in-control state. We assume that the target process $\left\{{\mathbf{Y}}_t\right\}$ is a weakly stationary process with mean vector $\mathit{\mu }$ and autocovariance matrix at lag h, denoted by $\mathbf{\Gamma }(h).$

Control charts present the mostly spread tool of SPC (see Montgomery 2020). While the first control charts were designed for detecting changes in the location behavior based on univariate independent observations (cf. Shewhart 1926; Page 1954; Roberts 1959), they were extended to time series by Alwan and Roberts (1988), Schmid (1995, 1997a), Schmid and Schöne (1997), and Knoth and Schmid (2002), among others. Another line of research led to the surveillance of the parameters in multivariate models.

The first multivariate control chart was proposed in Hotelling (1947), who introduced a control scheme based on the Mahalanobis distance to monitor the mean vector of the independent observations coming from the multivariate normal distribution. This approach was later extended by Crosier (1988), Pignatiello and Runger (1990), Lowry et al. (1992), and Ngai and Zhang (2001), who proposed several multivariate control charts based on the multivariate exponentially weighted moving average (MEWMA) recursion and the cumulative sum approach. Multivariate control charts for monitoring the parameters of multivariate time series have recently become a hot topic of research. Control charts of the mean behavior were discussed in Theodossiou (1993), Kramer and Schmid (1997), O. Bodnar and Schmid (2007, 2011), and O. Bodnar (2009a), and control charts for monitoring the covariance matrix were introduced in Śliwa and Schmid (2005) and O. Bodnar and Schmid (2017).

Due to the rapid development of computer technology, monitoring the parameters of complex high-dimensional processes has become possible and has attracted many researchers to this challenging field of science. In the high-dimensional setting, it is assumed that the dimensions of the data-generating model grow at the same rate as the sample size when the latter tends to infinity (see, e.g., Bai and Silverstein 2010; T. Bodnar, Dette, and Parolya 2019). K. Wang and Jiang (2009) considered a variable selection–based multivariate SPC procedure under the high-dimensional setting, and a high-dimensional control chart for profile monitoring was suggested by Chen and Nembhard (2011). The control scheme is based on the adaptive Neyman test statistic for the coefficients of the discrete Fourier transform of profiles. Li et al. (2014) suggested a new control 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. Z. Wang, Li, and Zhou (2017) constructed a hybrid control chart in the case of independent multivariate Poisson data, and R. Bodnar, Bodnar, and Schmid (2023) introduced several MEWMA-type control charts for high-dimensional time series where the Euclidean distance and the diagonalized Euclidean distance are employed in the construction of the control statistics instead of the Mahalanobis distance.

All of the abovementioned control charts were designed under the assumption that the parameters of the target process $\left\{{\mathbf{Y}}_t\right\}$ are known. However, this assumption appears to be very restrictive in many practical situations (see, e.g., Kramer and Schmid 2000; Albers and Kallenberg 2004; Jensen et al. 2006; Saleh et al. 2015; Jardim, Chakraborti, and Epprecht 2020; Sarmiento et al. 2022). In practice, the target process should be fitted in the Phase I analysis of the monitoring procedure, whereas the control statistics for the Phase II are constructed by replacing the unknown true parameters with their corresponding estimators. This approach leads to the misspecified control charts, and the effect of the misspecification has to be studied before the monitoring scheme is applied in practice. In this article, we contribute to the literature by developing new theoretical results that allow assessment of the effect of the misspecification in the high-dimensional setting. In particular, we show that the MEWMA control charts based on the Euclidean distance and the diagonalized Euclidean distance are quite robust to the misspecification caused by the estimation of the parameters of the target process in Phase I, whereas the control schemes based on the Mahalanobis distance can be strongly affected by the misspecification effect in high dimensions.

The derived results can be applied in several fields of science. One direction of possible applications lead to economics and finance with a special emphasis on optimal portfolio theory. Several control charts for monitoring the structure of optimal portfolios were suggested in O. Bodnar and Schmid (2007), Golosnoy and Schmid (2007), O. Bodnar (2009b), and Golosnoy et al. (2011), among others. Whereas these procedures were developed in the case of small dimensions of data-generating model, the introduced approaches extend the existent methods to the high-dimensional case. Moreover, because the model parameters of the data-generating model in Phase I are not known in financial applications, the obtained findings potentially introduce new methods for monitoring the structure of optimal portfolios, where the Euclidean distance and the diagonalized Euclidean distance are employed in the construction of the test statistics. Another line of possible applications is in environmental science, where problems of monitoring the parameters of the high-dimensional spatial processes may be present (see, e.g., Otto and Schmid 2023).

The rest of the article is organized as follows. In Section 2, the MEWMA recursion is discussed and its basic properties are presented. Section 3 introduces the MEWMA control charts for high-dimensional time series under model misspecification, and the distributional properties of the considered statistics are derived in Section 4. The results of the simulation study are given in Section 5, and final remarks are summarized in Section 6.

2. CONTROL CHARTS BASED ON MEWMA RECURSIONS

The MEWMA recursion is defined by (2.1) $${\mathbf{Z}}_t=(\mathbf{I}-\mathbf{R}){\mathbf{Z}}_{t-1}+\mathbf{R}{\mathbf{X}}_t,\hspace{1em}t\ge 1$$(2.1) with ${\mathbf{Z}}_0=\mathit{\mu }.$ In the following, we set $\mathbf{R}=\text{diag}(r_1,\dots ,r_p)$ with $r_i\in (0,1]$ for $i=1,\dots ,p$ being known and deterministic (see, e.g., Qiu 2013; Montgomery 2020).

Following Kramer and Schmid (1997), it holds that for $t\in N$ and p fixed, $$\mathbb{E}({\mathbf{Z}}_t)=\mathit{\mu }+{\mathbf{a}}_{t-\tau }\underset{t\to \infty }{\overset{}{\to }}\mathit{\mu }+\mathbf{a}{I}_N(\tau )\hspace{1em}\text{with}\hspace{1em}{\mathbf{a}}_{t-\tau }=(\mathbf{I}-{(\mathbf{I}-\mathbf{R})}^{t-\tau +1})\mathbf{a}{I}_{\{0,1,\dots \}}(t-\tau )$$ and $$\begin{array}{l}{\mathbf{\Sigma }}_t=\mathbb{C}\mathit{\text{ov}}({\mathbf{Z}}_t)=\mathbf{R}\sum _{i,j=0}^{t-1}{(\mathbf{I}-\mathbf{R})}^i\mathbf{\Gamma }(j-i){(\mathbf{I}-\mathbf{R})}^j\mathbf{R}\\ \underset{t\to \infty }{\overset{}{\to }}{\mathbf{\Sigma }}_l=\mathbf{R}\sum _{i,j=0}^{\infty }{(\mathbf{I}-\mathbf{R})}^i\mathbf{\Gamma }(j-i){(\mathbf{I}-\mathbf{R})}^j\mathbf{R},\end{array}$$ provided that $\{\mathbf{\Gamma }(v)\}$ is absolutely summable; that is, that ${\sum }_{v=0}^{\infty }\left|\right|\mathbf{\Gamma }(v)\left|\right|<\infty ,$ where $\left|\right|.\left|\right|$ denotes the Euclidean norm (see, e.g., Kramer and Schmid 1997). Furthermore, to ensure that ${\mathbf{\Sigma }}_t$ is positive definite, it is assumed that the covariance matrix of $({\mathbf{Y}}_1^{\prime },\dots {\mathbf{Y}}_t^{\prime })\prime ,$ defined as a block matrix with the (j, i)th block equal to $\mathbf{\Gamma }(j-i),$ is positive definite for any t. Alternatively, the positive definiteness of ${\mathbf{\Sigma }}_t$ can be ensured by assuming that the process $\left\{{\mathbf{Y}}_t\right\}$ has a positive definite spectral density matrix denoted by $\mathbf{f}(\lambda ),$ because for any vector u it holds that $$\begin{array}{l}\mathbf{u}\prime {\mathbf{\Sigma }}_t\mathbf{u}=\mathbf{u}\prime \mathbf{R}\sum _{v,j=0}^{t-1}{(\mathbf{I}-\mathbf{R})}^v\left({\int }_{-\pi }^{\pi }{e}^{i(j-v)\lambda }\mathbf{f}(\lambda )\mathrm{d}\lambda \right){(\mathbf{I}-\mathbf{R})}^j\mathbf{R}\mathbf{u}\\ ={\int }_{-\pi }^{\pi }\left(\sum _{v=0}^{t-1}{(\mathbf{I}-\mathbf{R})}^v{e}^{-\mathit{\text{iv}}\lambda }\mathbf{R}\mathbf{u}\right)\prime \mathbf{f}(\lambda )\left(\sum _{j=0}^{t-1}{(\mathbf{I}-\mathbf{R})}^j{e}^{-\mathit{\text{ij}}\lambda }\mathbf{R}\mathbf{u}\right)\mathrm{d}\lambda \ge 0\end{array}$$ with the equality if and only if $$\sum _{j=0}^{t-1}{(\mathbf{I}-\mathbf{R})}^j{e}^{-\mathit{\text{ij}}\lambda }\mathbf{R}\mathbf{u}={\left(\mathbf{I}-(\mathbf{I}-\mathbf{R})e^{-i\lambda }\right)}^{-1}\left(\mathbf{I}-{(\mathbf{I}-\mathbf{R})}^t{e}^{-\mathit{\text{it}}\lambda }\right)\mathbf{u}=\mathbf{0},$$ which is equivalent to $\mathbf{u}=\mathbf{0}$ due to the definition of R at the beginning of Section 2. It is noted that the spectral density matrix $\mathbf{f}(\lambda )$ is positive definite for many stationary multivariate time series, such as, for example, for a stationary vector autoregressive moving average (VARMA) process. Similar results in the univariate case can be found in Schmid (1997b). Whereas the mean vector of ${\mathbf{Z}}_t$ depends whether the observed process is in the in-control state or in the out-of-control state, its covariance matrix is same in both states. In the following, we will always assume that $\text{rk}({\mathbf{\Sigma }}_t)=p=\text{rk}({\mathbf{\Sigma }}_l).$

The first control chart based on the MEWMA recursion was suggested by Lowry et al. (1992) for independent multivariate observations, and Kramer and Schmid (1997) extended this chart to monitor changes in the mean vector of a multivariate time series. This control chart is based on the Mahalanobis distance $({\mathbf{Z}}_t-\mathit{\mu })\prime {\mathbf{\Sigma }}_t^{-1}({\mathbf{Z}}_t-\mathit{\mu }).$ Because for this scheme the covariance matrix has to be determined at each time point, practitioners prefer to work with the statistic $({\mathbf{Z}}_t-\mathit{\mu })\prime {\mathbf{\Sigma }}_l^{-1}({\mathbf{Z}}_t-\mathit{\mu }),$ where the exact covariance matrix of ${\mathbf{Z}}_t$ is replaced by its limit as t tends to infinity.

Recently, R. Bodnar, Bodnar, and Schmid (2023) proposed two further versions of control statistics based on the MEWMA recursion. The first one uses the Euclidean distance to compute the control statistic at each time t and is given by $({\mathbf{Z}}_t-\mathit{\mu })\prime ({\mathbf{Z}}_t-\mathit{\mu }).$ Let ${\mathbf{\Sigma }}_{d;t}$ be a diagonal matrix that consists of the diagonal elements of ${\mathbf{\Sigma }}_t.$ The second control statistic proposed in R. Bodnar, Bodnar, and Schmid (2023) is based on the diagonalized Euclidean distance expressed as $({\mathbf{Z}}_t-\mathit{\mu })\prime {\mathbf{\Sigma }}_{d;t}^{-1}({\mathbf{Z}}_t-\mathit{\mu }).$

Using these approaches. R. Bodnar, Bodnar, and Schmid (2023) developed suitable control statistics by normalizing these quantities. They distinguished between several possibilities. The control statistics are centered by subtracting the exact in-control mean, the asymptotic in-control mean as t goes to infinity, and the asymptotic in-control mean as p goes to infinity, respectively. Further they normalized these expressions by dividing by the square root of the exact variance, the asymptotic variance as t goes to infinity, and the asymptotic variance as p goes to infinity, respectively. Altogether, 12 control schemes were considered and compared with each other.

3. MEWMA CONTROL CHARTS WITH ESTIMATED PARAMETERS

The control charts in the previous section depend on certain parameters as—for example, $\mathit{\mu },\mathbf{\Gamma }(h),h\ge 0$—which were assumed to be known in the previous studies. In practice, however, these parameters have to be estimated. This is usually done by a prerun in engineering or by using a historical sample in, for example, finance. In the Phase I analysis, the unknown parameters of the target process are estimated by a previous sample, and these parameter estimators are used within the Phase II analysis, the monitoring phase.

We will assume in the following that a sample of the underlying target process—that is, for $\tau =\infty$—is available to estimate the parameters in Phase I. This sample is assumed to be independent of the observation vectors in Phase II, which are used to construct the MEWMA recursion and the corresponding control charts to monitor the mean behavior of the underlying high-dimensional time series.

To run a control chart, one should estimate the expectation and the autocovariance matrices of the underlying stationary process using the sample from Phase I by a certain estimation procedure. As such, the estimators of $\mathit{\mu }$ and $\mathbf{\Gamma }(h),h\ge 0$ are obtained and are used instead of the true population counterparts in the definition of the control statistics. In particular, this means that the estimators of $\mathit{\mu }$ and $\mathbf{\Gamma }(h),$ denoted by ${\mathit{\mu }}^{*}$ and $\mathbf{\Gamma }{(h)}^{*},$ are assumed to be deterministic in practice. Of course, we cannot estimate $\mathbf{\Gamma }(h)$ for all values of h for an arbitrary stationary process. However, it is usually assumed that the underlying target process $\left\{{\mathbf{Y}}_t\right\}$ follows a vector autoregressive (VAR) or VARMA process and then the autocovariance matrices can be estimated for all h using the estimators of the coefficient matrices and the covariance matrix of the white noise process.

If $\left\{{\mathbf{Y}}_t\right\}$ follows a VAR(1) process given by (3.1) $${\mathbf{Y}}_t=\mathit{\mu }+\mathbf{\Phi }({\mathbf{Y}}_{t-1}-\mathit{\mu })+{\mathit{\epsilon }}_t,$$(3.1) where $\left\{{\mathit{\epsilon }}_t\right\}$ are independent and normally distributed with $\mathbb{E}({\mathit{\epsilon }}_t)=\mathbf{0}$ and $\mathbb{C}\mathit{\text{ov}}({\mathit{\epsilon }}_t,{\mathit{\epsilon }}_t)=\mathbf{\Sigma },$ then it holds that $\mathbb{E}({\mathbf{Y}}_t)=\mathit{\mu }={\mathbb{E}}_{\infty }({\mathbf{X}}_t)$ and (see, e.g., Brockwell and Davis 1991; Reinsel 1993; Lütkepohl 2005) (3.2) $$\mathbf{\Gamma }(h)={\mathbf{\Phi }}^h\mathbf{\Gamma }(0)\hspace{1em}\text{and}\hspace{1em}\mathbf{\Gamma }(-h)=\mathbf{\Gamma }(h)\prime \hspace{1em}\text{for}\hspace{1em}h=1,2,\dots ,$$(3.2) where $\mathbf{\Gamma }(0)$ is the solution of the following matrix equation: (3.3) $$\mathbf{\Gamma }(0)=\mathbf{\Phi }\mathbf{\Gamma }(0)\mathbf{\Phi }\prime +\mathbf{\Sigma }.$$(3.3)

For estimation of the autocovariance matrices $\mathbf{\Gamma }(h)$ for $h\ge 0,$ the coefficient matrix $\mathbf{\Phi }$ and the covariance matrix $\mathbf{\Sigma }$ of the white noise process’s various procedures were proposed in the statistical literature (e.g., Brockwell and Davis 1991; Hamilton 1994). In the following, three estimation methods are considered.

Let ${\mathbf{X}}_{1:N_0}=({\mathbf{X}}_1,\dots ,{\mathbf{X}}_{N_0})$ denote the $p\times N_0$ observation matrix consisting of N₀ observation vectors taken in Phase I. Later, it is assumed that $N_0>p.$ This assumption ensures that the sample covariance matrix is positive definite with probability 1 when an independent sample is taken from the normal distribution (see, e.g., theorem 3.1.4 in Muirhead 1982). In the first and second approaches, the mean vector $\mathit{\mu }$ is estimated by its sample counterpart given by (Brockwell and Davis 1991) $$\overline{\mathbf{X}}=\frac{1}{N_0}\sum _{t=1}^{N_0}{\mathbf{X}}_t.$$

The first method is based on the nonparametric estimation of the autocovariance matrices $\mathbf{\Gamma }(h)$ expressed as $${\widehat{\mathbf{\Gamma }}}_{\mathit{\text{non}}}(h)=\{\begin{array}{ll}\frac{1}{N_0}\sum _{t=1}^{N_0-h}({\mathbf{X}}_{t+h}-\overline{\mathbf{X}})({\mathbf{X}}_t-\overline{\mathbf{X}})\prime & \text{for} 0\le h\le N_0-1,\\ \frac{1}{N_0}\sum _{t=-h+1}^{N_0}({\mathbf{X}}_{t+h}-\overline{\mathbf{X}})({\mathbf{X}}_t-\overline{\mathbf{X}})\prime & \text{for}-N_0+1\le h<0.\end{array}$$

The second approach uses the assumption of the VAR(1) model as given in (3.1)(3.1) $${\mathbf{Y}}_t=\mathit{\mu }+\mathbf{\Phi }({\mathbf{Y}}_{t-1}-\mathit{\mu })+{\mathit{\epsilon }}_t,$$(3.1) . In this case, $\mathbf{\Gamma }(0)$ and $\mathbf{\Gamma }(1)$ are estimated nonparametrically by ${\widehat{\mathbf{\Gamma }}}_{\mathit{\text{non}}}(0)$ and ${\widehat{\mathbf{\Gamma }}}_{\mathit{\text{non}}}(1),$ respectively, which are then used to estimate $\mathbf{\Phi }$ and $\mathbf{\Sigma }$ by $${\widehat{\mathbf{\Phi }}}_V={\widehat{\mathbf{\Gamma }}}_{\mathit{\text{non}}}(1){\widehat{\mathbf{\Gamma }}}_{\mathit{\text{non}}}{(0)}^{-1}$$ and $${\widehat{\mathbf{\Sigma }}}_V={\widehat{\mathbf{\Gamma }}}_{\mathit{\text{non}}}(0)-{\widehat{\mathbf{\Phi }}}_V{\widehat{\mathbf{\Gamma }}}_{\mathit{\text{non}}}(0){\widehat{\mathbf{\Phi }}}_V^{\prime }.$$

Finally, $\mathbf{\Gamma }(h)$ and $\mathbf{\Gamma }(-h)$ for $h\ge 1$ are estimated by $${\widehat{\mathbf{\Gamma }}}_V(h)={\widehat{\mathbf{\Phi }}}_V^h{\widehat{\mathbf{\Gamma }}}_{\mathit{\text{non}}}(0)\hspace{1em}\text{and}\hspace{1em}{\widehat{\mathbf{\Gamma }}}_V(-h)={\widehat{\mathbf{\Gamma }}}_V(h)\prime .$$

The third approach is based on the maximum likelihood estimation of $\mathit{\nu }=(\mathbf{I}-\mathbf{\Phi })\mathit{\mu },\mathbf{\Phi }$ and $\mathbf{\Sigma },$ where I denotes the identity matrix. Following section 11.1 of Hamilton (1994), the maximum likelihood estimator of $\mathbf{\Pi }\prime =[\mathit{\nu },\mathbf{\Phi }]$ is given by $$\widehat{\mathbf{\Pi }}{\prime }_{\mathit{\text{ML}}}=[{\widehat{\mathit{\nu }}}_{\mathit{\text{ML}}},{\widehat{\mathbf{\Phi }}}_{\mathit{\text{ML}}}]=\left[\sum _{t=1}^{N_0}{\mathbf{X}}_t{\mathbf{V}}_t^{\prime }\right]{\left[\sum _{t=1}^{N_0}{\mathbf{V}}_t{\mathbf{V}}_t^{\prime }\right]}^{-1}\hspace{1em}\text{with}\hspace{1em}{\mathbf{V}}_t=\left[\begin{array}{c}1\\ {\mathbf{X}}_{t-1}\end{array}\right].$$

Furthermore, the maximum likelihood estimator for $\mathbf{\Sigma }$ is expressed as $${\widehat{\mathbf{\Sigma }}}_{\mathit{\text{ML}}}=\frac{1}{T}\sum _{t=1}^T\widehat{{\mathit{\epsilon }}_t}\widehat{{\mathit{\epsilon }}_t}\prime \hspace{1em}\text{with}\hspace{1em}\widehat{{\mathit{\epsilon }}_t}={\mathbf{Y}}_t-\widehat{\mathbf{\Pi }}\prime {\mathbf{V}}_t.$$

Finally, using (3.2)(3.2) $$\mathbf{\Gamma }(h)={\mathbf{\Phi }}^h\mathbf{\Gamma }(0)\hspace{1em}\text{and}\hspace{1em}\mathbf{\Gamma }(-h)=\mathbf{\Gamma }(h)\prime \hspace{1em}\text{for}\hspace{1em}h=1,2,\dots ,$$(3.2) and (3.3)(3.3) $$\mathbf{\Gamma }(0)=\mathbf{\Phi }\mathbf{\Gamma }(0)\mathbf{\Phi }\prime +\mathbf{\Sigma }.$$(3.3) , we get the maximum likelihood estimators of $\mathbf{\Gamma }(h)$ given by $$\text{vec}({\widehat{\mathbf{\Gamma }}}_{\mathit{\text{ML}}}(0))=(\mathbf{I}-{\widehat{\mathbf{\Phi }}}_{\mathit{\text{ML}}}\otimes {\widehat{\mathbf{\Phi }}}_{\mathit{\text{ML}}})\text{vec}({\widehat{\mathbf{\Sigma }}}_{\mathit{\text{ML}}}),{\widehat{\mathbf{\Gamma }}}_{\mathit{\text{ML}}}(h)={\widehat{\mathbf{\Phi }}}_{\mathit{\text{ML}}}^h{\widehat{\mathbf{\Gamma }}}_{\mathit{\text{ML}}}(0),{\widehat{\mathbf{\Gamma }}}_{\mathit{\text{ML}}}(-h)={\widehat{\mathbf{\Gamma }}}_{\mathit{\text{ML}}}(h)\prime ,$$ where $\text{vec}(.)$ denotes the vec operator and the symbol ⊗ stands for the Kronecker product (see, e.g., Harville 1997).

Later, we use the symbol * to denote the estimators of the population quantities, which are used in the construction of the control statistics based on MEWMA recursion. This leads to ${\mathit{\mu }}^{*},$ $${\mathbf{\Sigma }}_t^{*}=\mathbf{R}\sum _{i,j=0}^{t-1}{(\mathbf{I}-\mathbf{R})}^i\mathbf{\Gamma }{(j-i)}^{*}{(\mathbf{I}-\mathbf{R})}^j\mathbf{R},$$ and ${\mathbf{\Sigma }}_{d;t}^{*}$ for the matrix consisting of the diagonal elements of ${\mathbf{\Sigma }}_t^{*}.$ By analogy, ${\mathbf{\Sigma }}_l^{*}$ is used, which is well defined if ${\sum }_{v=0}^{\infty }\left|\right|\mathbf{\Gamma }{(v)}^{*}\left|\right|\le \infty .$ Further, we use ${\mathbf{\Sigma }}_{d;l}^{*}$ for the limit of ${\mathbf{\Sigma }}_{d;t}^{*}$ as t tends to infinity.

Replacing $\mathit{\mu },{\mathbf{\Sigma }}_t,{\mathbf{\Sigma }}_l,$ and ${\mathbf{\Sigma }}_{d;t}$ by the misspecified values ${\mathit{\mu }}^{*},{\mathbf{\Sigma }}_t^{*},{\mathbf{\Sigma }}_l^{*},$ and ${\mathbf{\Sigma }}_{d;t}^{*}$ in the quadratic forms discussed in Section 2, we get the misspecified Mahalanobis quantities $$({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_t^{*\mathrm{ }-1}({\mathbf{Z}}_t-{\mathit{\mu }}^{*}),({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_l^{*\mathrm{ }-1}({\mathbf{Z}}_t-{\mathit{\mu }}^{*})$$ and the misspecified quantities based on the Euclidean distance $$({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime ({\mathbf{Z}}_t-{\mathit{\mu }}^{*}),({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}({\mathbf{Z}}_t-{\mathit{\mu }}^{*}).$$

R. Bodnar, Bodnar, and Schmid (2023) introduced several control charts based on the Euclidean distance. However, the effect of misspecification was not taken into account. In this article, we want to analyze the behavior of these control charts under misspecification. Following R. Bodnar, Bodnar, and Schmid (2023) and their notation, we consider the MEWMA control charts given by $$\begin{array}{l}{T}_{1,t}^{*}=\frac{({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime ({\mathbf{Z}}_t-{\mathit{\mu }}^{*})-\text{tr}({\mathbf{\Sigma }}_t^{*})}{\sqrt{2\text{tr}({\mathbf{\Sigma }}_t^{*\mathrm{ }2})}},\\ T_{2,t}^{*}=\frac{({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime ({\mathbf{Z}}_t-{\mathit{\mu }}^{*})-\text{tr}({\mathbf{\Sigma }}_l^{*})}{\sqrt{2\text{tr}({\mathbf{\Sigma }}_t^{*\mathrm{ }2})}},\\ T_{3,t}^{*}=\frac{({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime ({\mathbf{Z}}_t-{\mathit{\mu }}^{*})-\text{tr}({\mathbf{\Sigma }}_t^{*})}{\sqrt{2\text{tr}({\mathbf{\Sigma }}_l^{*\mathrm{ }2})}},\\ T_{4,t}^{*}=\frac{({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime ({\mathbf{Z}}_t-{\mathit{\mu }}^{*})-\text{tr}({\mathbf{\Sigma }}_l^{*})}{\sqrt{2\text{tr}({\mathbf{\Sigma }}_l^{*\mathrm{ }2})}},\\ T_{6,t}^{*}=\frac{({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}({\mathbf{Z}}_t-{\mathit{\mu }}^{*})-\text{tr}({\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t^{*})}{\sqrt{2\text{tr}({({\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t^{*})}^2)}},\\ T_{7,t}^{*}=\frac{({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}({\mathbf{Z}}_t-{\mathit{\mu }}^{*})-\text{tr}({\mathbf{\Sigma }}_{d;l}^{*\mathrm{ }-1}{\mathbf{\Sigma }}_l^{*})}{\sqrt{2\text{tr}({({\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t^{*})}^2)}},\\ T_{8,t}^{*}=\frac{({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}({\mathbf{Z}}_t-{\mathit{\mu }}^{*})-\text{tr}({\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t^{*})}{\sqrt{2\text{tr}({({\mathbf{\Sigma }}_{d;l}^{*\mathrm{ }-1}{\mathbf{\Sigma }}_l^{*})}^2)}},\\ T_{9,t}^{*}=\frac{({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}({\mathbf{Z}}_t-{\mathit{\mu }}^{*})-\text{tr}({\mathbf{\Sigma }}_{d;l}^{*\mathrm{ }-1}{\mathbf{\Sigma }}_l^{*})}{\sqrt{2\text{tr}({({\mathbf{\Sigma }}_{d;l}^{*\mathrm{ }-1}{\mathbf{\Sigma }}_l^{*})}^2)}},\\ T_{\mathit{\text{Mah}},t}^{*}=\frac{({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_t^{*\mathrm{ }-1}({\mathbf{Z}}_t-{\mathit{\mu }}^{*})-p}{\sqrt{2p}},\\ T_{\mathit{\text{MahInf}},t}^{*}=\frac{({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_l^{*\mathrm{ }-1}({\mathbf{Z}}_t-{\mathit{\mu }}^{*})-\text{tr}({\mathbf{\Sigma }}_l^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t^{*})}{\sqrt{2\text{tr}({({\mathbf{\Sigma }}_l^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t^{*})}^2)}}.\end{array}$$

These statistics were studied in R. Bodnar, Bodnar, and Schmid (2023); however, here we additionally take misspecification into account. The statistics $T_{5,t}$ and $T_{10,t}$ of R. Bodnar, Bodnar, and Schmid (2023) are not considered in this article, because they make use of the limit of the first and second moments as p tends to infinity, which in most cases is difficult to determine.

4. BEHAVIOR OF THE CONTROL STATISTICS UNDER MISSPECIFICATION

In this section, we analyze the distributional properties of the control statistics defined under misspecification in detail. In particular, the exact distributions of the quadratic forms present in the definition of the control statistics are derived and their asymptotic approximations are provided in two cases, when t tends to infinity and when p tends to infinity. These results shed light on the effect of misspecification on the performance of the considered MEWMA control charts, especially in the high-dimensional case.

The following notation will be used throughout the article: (4.1) $${\mathbf{a}}_{t-\tau }=(\mathbf{I}-{(\mathbf{I}-\mathbf{R})}^{t-\tau +1})\mathbf{a}{I}_{\{0,1,2,\dots \}}(t-\tau ).$$(4.1)

In Theorems 4.1 and 4.2, the results are presented in the case of the control schemes based on the Mahalanobis distance, and Theorems 4.3 and 4.4 provide the results of the control chart based on the Euclidean distance and the diagonalized Euclidean distance, respectively.

Theorem 4.1.

Let $\left\{{\mathbf{Y}}_t\right\}$ be a stationary Gaussian process with $\mathbb{E}({\mathbf{Y}}_t)=\mathit{\mu }$ and $\mathbb{C}\mathit{\text{ov}}({\mathbf{Y}}_{t+h},{\mathbf{Y}}_t)=\mathbf{\Gamma }(h).$ Let τ be fixed.

1. Suppose that $\text{rk}({\mathbf{\Sigma }}_t)=\text{rk}({\mathbf{\Sigma }}_t^{*})=p$. Let ${\mathbf{U}}_t$ be an orthogonal matrix such that$${\mathbf{U}}_t^{\prime }{\mathbf{\Sigma }}_t^{1/2}{\mathbf{\Sigma }}_t^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t^{1/2}{\mathbf{U}}_t=\text{diag}({\lambda }_{\mathit{\text{Mah}},1,t},\dots ,{\lambda }_{\mathit{\text{Mah}},p,t})$$and let ${\mathit{\delta }}_{\mathit{\text{Mah}},t}={\mathbf{U}}_t^{\prime }{\mathbf{\Sigma }}_t^{-1/2}(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})={({\delta }_{\mathit{\text{Mah}},i,t})}_{i=1,\dots ,p}$. Further, suppose that ζ₁, …, ζ_p are independent and standard normally distributed random variables.

   Then$$({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_t^{*\mathrm{ }-1}({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\stackrel{d}{=}\sum _{i=1}^p{\lambda }_{\mathit{\text{Mah}},i,t}{({\zeta }_i+{\delta }_{\mathit{\text{Mah}},i,t})}^2.$$Moreover,$$\mathbb{E}(({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_t^{*\mathrm{ }-1}({\mathbf{Z}}_t-{\mathit{\mu }}^{*}))=\text{tr}({\mathbf{\Sigma }}_t^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t)+(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_t^{*\mathrm{ }-1}(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})$$and$$\mathbb{V}\mathit{\text{ar}}(({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_t^{*\mathrm{ }-1}({\mathbf{Z}}_t-{\mathit{\mu }}^{*}))=2\text{tr}({({\mathbf{\Sigma }}_t^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t)}^2)+4(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_t^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t{\mathbf{\Sigma }}_t^{*-1}(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*}).$$

2. Suppose that $\text{rk}({\mathbf{\Sigma }}_t)=\text{rk}({\mathbf{\Sigma }}_l)=\text{rk}({\mathbf{\Sigma }}_t^{*})=\text{rk}({\mathbf{\Sigma }}_l^{*})=p$ and that $\{\mathbf{\Gamma }(v)\}$ and $\left\{\mathbf{\Gamma }{(v)}^{*}\right\}$ are absolutely summable. Let p be fixed and let ${\mathbf{U}}_l$ be an orthogonal matrix such that$${\mathbf{U}}_l^{\prime }{\mathbf{\Sigma }}_l^{1/2}{\mathbf{\Sigma }}_l^{*\mathrm{ }-1}{\mathbf{\Sigma }}_l^{1/2}{\mathbf{U}}_l=\text{diag}({\lambda }_{\mathit{\text{Mah}},1},\dots ,{\lambda }_{\mathit{\text{Mah}},p})$$and let ${\mathit{\delta }}_{\mathit{\text{Mah}}}={\mathbf{U}}_l^{\prime }{\mathbf{\Sigma }}_l^{-1/2}(\mathit{\mu }+\mathbf{a}{I}_N(\tau )-{\mathit{\mu }}^{*})={({\delta }_{\mathit{\text{Mah}},i})}_{i=1,\dots ,p}$.

   If, further,(4.2) $$\underset{t\to \infty }{\mathrm{\text{lim}}}{\mathbf{U}}_t={\mathbf{U}}_l\hspace{1em}\text{and}\hspace{1em}\underset{t\to \infty }{\mathrm{\text{lim}}}{\mathbf{\Sigma }}_t^{1/2}={\mathbf{\Sigma }}_l^{1/2},$$(4.2) then the asymptotic distribution of $({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_t^{*\mathrm{ }-1}({\mathbf{Z}}_t-{\mathit{\mu }}^{*})$ as t tends to infinity is equal to the distribution of ${\sum }_{i=1}^p{\lambda }_{\mathit{\text{Mah}},i}{({\zeta }_i+{\delta }_{\mathit{\text{Mah}},i})}^2$ with ζ_i as in part (a).

3. Let t be fixed. Suppose that $\text{rk}({\mathbf{\Sigma }}_t)=\text{rk}({\mathbf{\Sigma }}_t^{*})=p$ and that(4.3) $$\underset{p\to \infty }{\mathrm{\text{lim}}}\frac{{\mathrm{\text{max}}}_{1\le i\le p}{\lambda }_{\mathit{\text{Mah}},i,t}^2(1+2{\delta }_{\mathit{\text{Mah}},i,t}^2)}{\sum _{i=1}^p{\lambda }_{\mathit{\text{Mah}},i,t}^2(1+2{\delta }_{\mathit{\text{Mah}},i,t}^2)}=0.$$(4.3) Then$$\frac{({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_t^{*\mathrm{ }-1}({\mathbf{Z}}_t-{\mathit{\mu }}^{*})-\mathbb{E}(({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_t^{*\mathrm{ }-1}({\mathbf{Z}}_t-{\mathit{\mu }}^{*}))}{\sqrt{\mathbb{V}\mathit{\text{ar}}(({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_t^{*\mathrm{ }-1}({\mathbf{Z}}_t-{\mathit{\mu }}^{*}))}}\underset{p\to \infty }{\overset{d}{\to }}N(0,1).$$

Proof.

It holds that ${\mathbf{Z}}_t-\mathit{\mu }\sim {\mathcal{N}}_p({\mathbf{a}}_{t-\tau },{\mathbf{\Sigma }}_t)$ and thus ${\mathbf{\Sigma }}_t^{-1/2}({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\sim \mathcal{N}({\mathbf{\Sigma }}_t^{-1/2}(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*}),\mathbf{I}).$ Consequently, $$({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_t^{*\mathrm{ }-1}({\mathbf{Z}}_t-{\mathit{\mu }}^{*})=({\mathbf{\Sigma }}_t^{-1/2}({\mathbf{Z}}_t-{\mathit{\mu }}^{*}))\prime {\mathbf{\Sigma }}_t^{1/2}{\mathbf{\Sigma }}_t^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t^{1/2} ({\mathbf{\Sigma }}_t^{-1/2}({\mathbf{Z}}_t-{\mathit{\mu }}^{*})).$$

Thus, the proof of part (a) follows immediately using chapter 3.1a, corollary 3.2b.1, and theorem 3.2b.2 of Mathai and Provost (1992).

To prove part (b), we make use of (a) and the fact that ${\mathbf{\Sigma }}_t\underset{t\to \infty }{\to }{\mathbf{\Sigma }}_l.$ Further, we use that the eigenvalues of a matrix are continuous functions of the elements of the matrix (cf. theorem 9.6 in Lax 2007). Consequently, ${\lambda }_{\mathit{\text{Mah}},i,t}\underset{t\to \infty }{\to }{\lambda }_{\mathit{\text{Mah}},i}.$ Because of (4.2)(4.2) $$\underset{t\to \infty }{\mathrm{\text{lim}}}{\mathbf{U}}_t={\mathbf{U}}_l\hspace{1em}\text{and}\hspace{1em}\underset{t\to \infty }{\mathrm{\text{lim}}}{\mathbf{\Sigma }}_t^{1/2}={\mathbf{\Sigma }}_l^{1/2},$$(4.2) , part (b) follows.

To prove part (c), we apply lemma 7.1 of R. Bodnar, Bodnar, and Schmid (2023). □

Condition (4.2)(4.2) $$\underset{t\to \infty }{\mathrm{\text{lim}}}{\mathbf{U}}_t={\mathbf{U}}_l\hspace{1em}\text{and}\hspace{1em}\underset{t\to \infty }{\mathrm{\text{lim}}}{\mathbf{\Sigma }}_t^{1/2}={\mathbf{\Sigma }}_l^{1/2},$$(4.2) is needed to prove part (b). In principle, we need the eigenvectors of a symmetric matrix to be continuous functions of the elements of the matrix. This is in general not fulfilled and therefore we have to assume (4.2)(4.2) $$\underset{t\to \infty }{\mathrm{\text{lim}}}{\mathbf{U}}_t={\mathbf{U}}_l\hspace{1em}\text{and}\hspace{1em}\underset{t\to \infty }{\mathrm{\text{lim}}}{\mathbf{\Sigma }}_t^{1/2}={\mathbf{\Sigma }}_l^{1/2},$$(4.2) . A detailed discussion of this problem is given in, for example, chapter 9 of Lax (2007). A sufficient condition for (4.2)(4.2) $$\underset{t\to \infty }{\mathrm{\text{lim}}}{\mathbf{U}}_t={\mathbf{U}}_l\hspace{1em}\text{and}\hspace{1em}\underset{t\to \infty }{\mathrm{\text{lim}}}{\mathbf{\Sigma }}_t^{1/2}={\mathbf{\Sigma }}_l^{1/2},$$(4.2) to hold is that all eigenvalues of ${\mathbf{\Sigma }}_t^{1/2}{\mathbf{\Sigma }}_t^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t^{1/2}$ are simple. Condition (4.3)(4.3) $$\underset{p\to \infty }{\mathrm{\text{lim}}}\frac{{\mathrm{\text{max}}}_{1\le i\le p}{\lambda }_{\mathit{\text{Mah}},i,t}^2(1+2{\delta }_{\mathit{\text{Mah}},i,t}^2)}{\sum _{i=1}^p{\lambda }_{\mathit{\text{Mah}},i,t}^2(1+2{\delta }_{\mathit{\text{Mah}},i,t}^2)}=0.$$(4.3) is a technical one and is needed to apply a central limit theorem for a nonidentically distributed random sequence. In particular, this condition ensures that there is no dominating summand (with considerably larger variance) in the infinite sum of the random variables. Namely, the quadratic form $({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_t^{*\mathrm{ }-1}({\mathbf{Z}}_t-{\mathit{\mu }}^{*})$ can be presented as a weighted sum of independent random variables that are all non-central ${\chi }^2$-distributed with one degree of freedom. Furthermore, the denominator in condition (4.3)(4.3) $$\underset{p\to \infty }{\mathrm{\text{lim}}}\frac{{\mathrm{\text{max}}}_{1\le i\le p}{\lambda }_{\mathit{\text{Mah}},i,t}^2(1+2{\delta }_{\mathit{\text{Mah}},i,t}^2)}{\sum _{i=1}^p{\lambda }_{\mathit{\text{Mah}},i,t}^2(1+2{\delta }_{\mathit{\text{Mah}},i,t}^2)}=0.$$(4.3) is equal to the variance of the quadratic form with the summands corresponding the variance contribution of each random variable presented in the stochastic representation of $({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_t^{*\mathrm{ }-1}({\mathbf{Z}}_t-{\mathit{\mu }}^{*})$ mentioned above. As such, condition (4.3)(4.3) $$\underset{p\to \infty }{\mathrm{\text{lim}}}\frac{{\mathrm{\text{max}}}_{1\le i\le p}{\lambda }_{\mathit{\text{Mah}},i,t}^2(1+2{\delta }_{\mathit{\text{Mah}},i,t}^2)}{\sum _{i=1}^p{\lambda }_{\mathit{\text{Mah}},i,t}^2(1+2{\delta }_{\mathit{\text{Mah}},i,t}^2)}=0.$$(4.3) requires that no summand in the variance decomposition be dominating.

Now we study the statistic based on the limit covariance matrix.

Theorem 4.2.

Let $\left\{{\mathbf{Y}}_t\right\}$ be a stationary Gaussian process with $\mathbb{E}({\mathbf{Y}}_t)=\mathit{\mu }$ and $\mathbb{C}\mathit{\text{ov}}({\mathbf{Y}}_{t+h},{\mathbf{Y}}_t)=\mathbf{\Gamma }(h).$ Let τ be fixed.

1. Suppose that $\text{rk}({\mathbf{\Sigma }}_t)=\text{rk}({\mathbf{\Sigma }}_t^{*})=p$. Let ${\mathbf{U}}_{l;t}$ be an orthogonal matrix such that$${\mathbf{U}}_{l;t}^{\prime }{\mathbf{\Sigma }}_t^{1/2}{\mathbf{\Sigma }}_l^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t^{1/2}{\mathbf{U}}_{l;t}=\text{diag}({\lambda }_{\mathit{\text{MahInf}},1,t},\dots ,{\lambda }_{\mathit{\text{MahInf}},p,t})$$and let ${\mathit{\delta }}_{\mathit{\text{MahInf}},t}={\mathbf{U}}_{l;t}^{\prime }{\mathbf{\Sigma }}_t^{-1/2}(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})={({\delta }_{\mathit{\text{MahInf}},i,t})}_{i=1,\dots ,p}$. Further, suppose that ζ₁, …, ζ_p are independent and standard normally distributed random variables.

   Then$$({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_l^{*\mathrm{ }-1}({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\stackrel{d}{=}\sum _{i=1}^p{\lambda }_{\mathit{\text{MahInf}},i,t}{({\zeta }_i+{\delta }_{\mathit{\text{MahInf}},i,t})}^2.$$Moreover,$$\mathbb{E}(({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_l^{*\mathrm{ }-1}({\mathbf{Z}}_t-{\mathit{\mu }}^{*}))=\text{tr}({\mathbf{\Sigma }}_l^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t)+(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_l^{*\mathrm{ }-1}(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})$$and$$\mathbb{V}\mathit{\text{ar}}(({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_l^{*\mathrm{ }-1}({\mathbf{Z}}_t-{\mathit{\mu }}^{*}))=2\text{tr}({({\mathbf{\Sigma }}_l^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t)}^2)+4(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_l^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t{\mathbf{\Sigma }}_l^{*-1}(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*}).$$

2. Suppose that $\text{rk}({\mathbf{\Sigma }}_t)=\text{rk}({\mathbf{\Sigma }}_l)=\text{rk}({\mathbf{\Sigma }}_t^{*})=\text{rk}({\mathbf{\Sigma }}_l^{*})=p$ and that $\{\mathbf{\Gamma }(v)\}$ and $\left\{\mathbf{\Gamma }{(v)}^{*}\right\}$ are absolutely summable. Let p be fixed and let ${\mathbf{U}}_l$ be the orthogonal matrix defined in Theorem 4.1(b). If, further,(4.4) $$\underset{t\to \infty }{\mathrm{\text{lim}}}{\mathbf{U}}_{l;t}={\mathbf{U}}_l\hspace{1em}\text{and}\hspace{1em}\underset{t\to \infty }{\mathrm{\text{lim}}}{\mathbf{\Sigma }}_t^{1/2}={\mathbf{\Sigma }}_l^{1/2},$$(4.4) then the asymptotic distribution of $({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_l^{*\mathrm{ }-1}({\mathbf{Z}}_t-{\mathit{\mu }}^{*})$ as t tends to infinity is equal to the distribution of ${\sum }_{i=1}^p{\lambda }_{\mathit{\text{Mah}},i}{({\zeta }_i+{\delta }_{\mathit{\text{Mah}},i})}^2$ with ζ_i as in part (a).

3. Let t be fixed. Suppose that $\text{rk}({\mathbf{\Sigma }}_t)=\text{rk}({\mathbf{\Sigma }}_t^{*})=p$ and that$$\underset{p\to \infty }{\mathrm{\text{lim}}}\frac{{\mathrm{\text{max}}}_{1\le i\le p}{\lambda }_{\mathit{\text{MahInf}},i,t}^2(1+2{\delta }_{\mathit{\text{MahInf}},i,t}^2)}{\sum _{i=1}^p{\lambda }_{\mathit{\text{MahInf}},i,t}^2(1+2{\delta }_{\mathit{\text{MahInf}},i,t}^2)}=0.$$Then$$\frac{({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_l^{*\mathrm{ }-1}({\mathbf{Z}}_t-{\mathit{\mu }}^{*})-\mathbb{E}(({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_l^{*\mathrm{ }-1}({\mathbf{Z}}_t-{\mathit{\mu }}^{*}))}{\sqrt{\mathbb{V}\mathit{\text{ar}}(({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_l^{*\mathrm{ }-1}({\mathbf{Z}}_t-{\mathit{\mu }}^{*}))}}\underset{p\to \infty }{\overset{d}{\to }}N(0,1).$$ Of course, the limit distribution in Theorem 4.2(b) is the same as that in Theorem 4.1(b).

Next we analyze the statistics based on the Euclidean distance and the diagonalized Euclidean distance.

Theorem 4.3.

Let $\left\{{\mathbf{Y}}_t\right\}$ be a stationary Gaussian process with $\mathbb{E}({\mathbf{Y}}_t)=\mathit{\mu }$ and $\mathbb{C}\mathit{\text{ov}}({\mathbf{Y}}_{t+h},{\mathbf{Y}}_t)=\mathbf{\Gamma }(h).$ Let τ be fixed.

1. Suppose that $\text{rk}({\mathbf{\Sigma }}_t)=p$. Let ${\tilde{\mathbf{U}}}_t$ be an orthogonal matrix such that$${\tilde{\mathbf{U}}}_t^{\prime }{\mathbf{\Sigma }}_t{\tilde{\mathbf{U}}}_t=\text{diag}({\lambda }_{\mathit{\text{Eu}},1,t},\dots ,{\lambda }_{\mathit{\text{Eu}},p,t})$$and let ${\mathit{\delta }}_{\mathit{\text{Eu}},t}={\tilde{\mathbf{U}}}_t^{\prime }{\mathbf{\Sigma }}_t^{-1/2}(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})={({\delta }_{\mathit{\text{Eu}},i,t})}_{i=1,\dots ,p}$. Further, suppose that ζ₁, …, ζ_p are independent and standard normally distributed random variables.

   Then$$({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime ({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\stackrel{d}{=}\sum _{i=1}^p{\lambda }_{\mathit{\text{Eu}},i,t}{({\zeta }_i+{\delta }_{\mathit{\text{Eu}},i,t})}^2.$$Moreover,$$\mathbb{E}(({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime ({\mathbf{Z}}_t-{\mathit{\mu }}^{*}))=\text{tr}({\mathbf{\Sigma }}_t)+(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_t(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})$$and$$\mathbb{V}\mathit{\text{ar}}(({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime ({\mathbf{Z}}_t-{\mathit{\mu }}^{*}))=2\text{tr}({\mathbf{\Sigma }}_t^2)+4(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_t^2(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*}).$$

2. Suppose that $\text{rk}({\mathbf{\Sigma }}_t)=\text{rk}({\mathbf{\Sigma }}_l)=p$ and that $\{\mathbf{\Gamma }(v)\}$ is absolutely summable. Let p be fixed and let ${\tilde{\mathbf{U}}}_l$ be an orthogonal matrix such that$${\tilde{\mathbf{U}}}_l^{\prime }{\mathbf{\Sigma }}_l{\tilde{\mathbf{U}}}_l=\text{diag}({\lambda }_{\mathit{\text{Eu}},1},\dots ,{\lambda }_{\mathit{\text{Eu}},p})$$and let ${\mathit{\delta }}_{\mathit{\text{Eu}}}={\tilde{\mathbf{U}}}_l^{\prime }{\mathbf{\Sigma }}_l^{-1/2}(\mathit{\mu }+\mathbf{a}{I}_N(\tau )-{\mathit{\mu }}^{*})={({\delta }_{\mathit{\text{Eu}},i})}_{i=1\dots ,p.}$ If, further,(4.5) $$\underset{t\to \infty }{\mathrm{\text{lim}}}{\tilde{\mathbf{U}}}_t={\tilde{\mathbf{U}}}_l\hspace{1em}\text{and}\hspace{1em}\underset{t\to \infty }{\mathrm{\text{lim}}}{\mathbf{\Sigma }}_t^{1/2}={\mathbf{\Sigma }}_l^{1/2},$$(4.5) then the asymptotic distribution of $({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime ({\mathbf{Z}}_t-{\mathit{\mu }}^{*})$ as t tends to infinity is equal to the distribution of ${\sum }_{i=1}^p{\lambda }_{\mathit{\text{Eu}},i}{({\zeta }_i+{\delta }_{\mathit{\text{Eu}},i})}^2$ with ζ_i as in part (a).

3. Let t be fixed. Suppose that $\text{rk}({\mathbf{\Sigma }}_t)=p$ and that$$\underset{p\to \infty }{\mathrm{\text{lim}}}\frac{{\mathrm{\text{max}}}_{1\le i\le p}{\lambda }_{\mathit{\text{Eu}},i,t}^2(1+2{\delta }_{\mathit{\text{Eu}},i,t}^2)}{\sum _{i=1}^p{\lambda }_{\mathit{\text{Eu}},i,t}^2(1+2{\delta }_{\mathit{\text{Eu}},i,t}^2)}=0.$$Then$$\frac{({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime ({\mathbf{Z}}_t-{\mathit{\mu }}^{*})-\mathbb{E}(({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime ({\mathbf{Z}}_t-{\mathit{\mu }}^{*}))}{\sqrt{\mathbb{V}\mathit{\text{ar}}(({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime ({\mathbf{Z}}_t-{\mathit{\mu }}^{*}))}}\underset{p\to \infty }{\overset{d}{\to }}N(0,1).$$

Theorem 4.4.

Let $\left\{{\mathbf{Y}}_t\right\}$ be a stationary Gaussian process with $\mathbb{E}({\mathbf{Y}}_t)=\mathit{\mu }$ and $\mathbb{C}\mathit{\text{ov}}({\mathbf{Y}}_{t+h},{\mathbf{Y}}_t)=\mathbf{\Gamma }(h).$ Let τ be fixed.

1. Suppose that $\text{rk}({\mathbf{\Sigma }}_t)=\text{rk}({\mathbf{\Sigma }}_t^{*})=p$. Let ${\mathbf{U}}_{d;t}$ be an orthogonal matrix such that$${\mathbf{U}}_{d;t}^{\prime }{\mathbf{\Sigma }}_t^{1/2}{\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t^{1/2}{\mathbf{U}}_{d;t}=\text{diag}({\lambda }_{\mathit{\text{Eu}},d;1,t},\dots ,{\lambda }_{\mathit{\text{Eu}},d;p,t})$$and let ${\mathit{\delta }}_{\mathit{\text{Eu}},d;t}={\mathbf{U}}_{d;t}^{\prime }{\mathbf{\Sigma }}_t^{-1/2}(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})={({\delta }_{\mathit{\text{Eu}},d;i,t})}_{i=1,\dots ,p}$. Further, suppose that ζ₁, …, ζ_p are independent and standard normally distributed random variables.

   Then$$({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\stackrel{d}{=}\sum _{i=1}^p{\lambda }_{\mathit{\text{Eu}},d;i,t}{({\zeta }_i+{\delta }_{\mathit{\text{Eu}},d;i,t})}^2.$$Moreover,$$\mathbb{E}(({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}({\mathbf{Z}}_t-{\mathit{\mu }}^{*}))=\text{tr}({\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t)+(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})$$and$$\mathbb{V}\mathit{\text{ar}}(({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}({\mathbf{Z}}_t-{\mathit{\mu }}^{*}))=2\text{tr}({({\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t)}^2)+4(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t{\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*}).$$

2. Suppose that $\text{rk}({\mathbf{\Sigma }}_t)=\text{rk}({\mathbf{\Sigma }}_l)=\text{rk}({\mathbf{\Sigma }}_t^{*})=\text{rk}({\mathbf{\Sigma }}_l^{*})=p$ and that $\{\mathbf{\Gamma }(v)\}$ and $\left\{\mathbf{\Gamma }{(v)}^{*}\right\}$ are absolutely summable. Let p be fixed and let ${\mathbf{U}}_d$ be an orthogonal matrix such that$${\mathbf{U}}_d^{\prime }{\mathbf{\Sigma }}_l^{1/2}{\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}{\mathbf{\Sigma }}_l^{1/2}{\mathbf{U}}_d=\text{diag}({\lambda }_{\mathit{\text{Eu}},d;1},\dots ,{\lambda }_{\mathit{\text{Eu}},d;p})$$and let ${\mathit{\delta }}_{\mathit{\text{Eu}},d}={\mathbf{U}}_d^{\prime }{\mathbf{\Sigma }}_l^{-1/2}(\mathit{\mu }+\mathbf{a}{I}_N(\tau )-{\mathit{\mu }}^{*})={({\delta }_{\mathit{\text{Eu}},d;i})}_{i=1,\dots ,p}.$

   If, further,(4.6) $$\underset{t\to \infty }{\mathrm{\text{lim}}}{\mathbf{U}}_{d;t}={\mathbf{U}}_d\hspace{1em}\text{and}\hspace{1em}\underset{t\to \infty }{\mathrm{\text{lim}}}{\mathbf{\Sigma }}_t^{1/2}={\mathbf{\Sigma }}_l^{1/2},$$(4.6) then the asymptotic distribution of $({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}({\mathbf{Z}}_t-{\mathit{\mu }}^{*})$ as t tends to infinity is equal to the distribution of ${\sum }_{i=1}^p{\lambda }_{\mathit{\text{Eu}},d;i}{({\zeta }_i+{\delta }_{\mathit{\text{Eu}},d;i})}^2$ with ζ_i as in part (a).

3. Let t be fixed. Suppose that $\text{rk}({\mathbf{\Sigma }}_t)=\text{rk}({\mathbf{\Sigma }}_t^{*})=p$ and that$$\underset{p\to \infty }{\mathrm{\text{lim}}}\frac{{\mathrm{\text{max}}}_{1\le i\le p}{\lambda }_{\mathit{\text{Eu}},d;i,t}^2(1+2{\delta }_{\mathit{\text{Eu}},d;i,t}^2)}{\sum _{i=1}^p{\lambda }_{\mathit{\text{Eu}},d;i,t}^2(1+2{\delta }_{\mathit{\text{Eu}},d;i,t}^2)}=0.$$Then$$\frac{({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}({\mathbf{Z}}_t-{\mathit{\mu }}^{*})-\mathbb{E}(({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}({\mathbf{Z}}_t-{\mathit{\mu }}^{*}))}{\sqrt{\mathbb{V}\mathit{\text{ar}}(({\mathbf{Z}}_t-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}({\mathbf{Z}}_t-{\mathit{\mu }}^{*}))}}\underset{p\to \infty }{\overset{d}{\to }}N(0,1).$$ Note that in practice the in-control mean is frequently known and the mean must not be estimated. This case is obtained by setting ${\mathit{\mu }}^{*}=\mathit{\mu }$ in the above formulas.

Using the results of the theorems, we get under the conditions given that as $t\to \infty$ and with p fixed, $$\begin{array}{l}P(T_{1,t}^{*}\le x)\approx P(\sum _{i=1}^p{\lambda }_{\mathit{\text{Eu}},i}{({\zeta }_i+{\delta }_{\mathit{\text{Eu}},i})}^2\le \text{tr}({\mathbf{\Sigma }}_t^{*})+x\sqrt{2\text{tr}({\mathbf{\Sigma }}_t^{*\mathrm{ }2})}),\\ P(T_{2,t}^{*}\le x)\approx P(\sum _{i=1}^p{\lambda }_{\mathit{\text{Eu}},i}{({\zeta }_i+{\delta }_{\mathit{\text{Eu}},i})}^2\le \text{tr}({\mathbf{\Sigma }}_l^{*})+x\sqrt{2\text{tr}({\mathbf{\Sigma }}_t^{*\mathrm{ }2})}),\\ P(T_{3,t}^{*}\le x)\approx P(\sum _{i=1}^p{\lambda }_{\mathit{\text{Eu}},i}{({\zeta }_i+{\delta }_{\mathit{\text{Eu}},i})}^2\le \text{tr}({\mathbf{\Sigma }}_t^{*})+x\sqrt{2\text{tr}({\mathbf{\Sigma }}_l^{*\mathrm{ }2}}),\\ P(T_{4,t}^{*}\le x)\approx P(\sum _{i=1}^p{\lambda }_{\mathit{\text{Eu}},i}{({\zeta }_i+{\delta }_{\mathit{\text{Eu}},i})}^2\le \text{tr}({\mathbf{\Sigma }}_l^{*})+x\sqrt{2\text{tr}({\mathbf{\Sigma }}_l^{*\mathrm{ }2})}),\\ P(T_{6,t}^{*}\le x)\approx P(\sum _{i=1}^p{\lambda }_{\mathit{\text{Eu}},d;i}{({\zeta }_i+{\delta }_{\mathit{\text{Eu}},d;i})}^2\le \text{tr}({\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t^{*})+x\sqrt{2\text{tr}({({\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t^{*})}^2)}),\\ P(T_{7,t}^{*}\le x)\approx P(\sum _{i=1}^p{\lambda }_{\mathit{\text{Eu}},d;i}{({\zeta }_i+{\delta }_{\mathit{\text{Eu}},d;i})}^2\le \text{tr}({\mathbf{\Sigma }}_{d;l}^{*\mathrm{ }-1}{\mathbf{\Sigma }}_l^{*})+x\sqrt{2\text{tr}({({\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t^{*})}^2)}),\\ P(T_{8,t}^{*}\le x)\approx P(\sum _{i=1}^p{\lambda }_{\mathit{\text{Eu}},d;i}{({\zeta }_i+{\delta }_{\mathit{\text{Eu}},d;i})}^2\le \text{tr}({\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t^{*})+x\sqrt{2\text{tr}({({\mathbf{\Sigma }}_{d;l}^{*\mathrm{ }-1}{\mathbf{\Sigma }}_l^{*})}^2)}),\\ P(T_{9,t}^{*}\le x)\approx P(\sum _{i=1}^p{\lambda }_{\mathit{\text{Eu}},d;i}{({\zeta }_i+{\delta }_{\mathit{\text{Eu}},d;i})}^2\le \text{tr}({\mathbf{\Sigma }}_{d;l}^{*\mathrm{ }-1}{\mathbf{\Sigma }}_l^{*})+x\sqrt{2\text{tr}({({\mathbf{\Sigma }}_{d;l}^{*\mathrm{ }-1}{\mathbf{\Sigma }}_l^{*})}^2)}),\\ P(T_{\mathit{\text{Mah}},t}^{*}\le x)\approx P(\sum _{i=1}^p{\lambda }_{\mathit{\text{Mah}},i}{({\zeta }_i+{\delta }_{\mathit{\text{Mah}},i})}^2\le p+x\sqrt{2p}),\\ P(T_{\mathit{\text{MahInf}},t}^{*}\le x)\approx P(\sum _{i=1}^p{\lambda }_{\mathit{\text{Mah}},i}{({\zeta }_i+{\delta }_{\mathit{\text{Mah}},i})}^2\le \text{tr}({\mathbf{\Sigma }}_l^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t^{*})+x\sqrt{2\text{tr}({({\mathbf{\Sigma }}_l^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t^{*})}^2)}).\end{array}$$

As t is fixed and p tends to infinity, it holds under conditions given in Section 4 that control statistics are asymptotically normally distributed in the high-dimensional setting. Namely, we get that $$\begin{array}{l}P(T_{1,t}^{*}\le x)\approx \Phi (\frac{\text{tr}({\mathbf{\Sigma }}_t^{*})-\text{tr}({\mathbf{\Sigma }}_t)-(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_t(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})+x\sqrt{2\text{tr}({\mathbf{\Sigma }}_t^{*\mathrm{ }2})}}{\sqrt{2\text{tr}({({\mathbf{\Sigma }}_t)}^2)+4(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_t^2(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})}}),\\ P(T_{2,t}^{*}\le x)\approx \Phi (\frac{\text{tr}({\mathbf{\Sigma }}_l^{*})-\text{tr}({\mathbf{\Sigma }}_t)-(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_t(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})+x\sqrt{2\text{tr}({\mathbf{\Sigma }}_t^{*\mathrm{ }2})}}{\sqrt{2\text{tr}({({\mathbf{\Sigma }}_t)}^2)+4(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_t^2(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})}}),\\ P(T_{3,t}^{*}\le x)\approx \Phi (\frac{\text{tr}({\mathbf{\Sigma }}_t^{*})-\text{tr}({\mathbf{\Sigma }}_t)-(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_t(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})+x\sqrt{2\text{tr}({\mathbf{\Sigma }}_l^{*\mathrm{ }2})}}{\sqrt{2\text{tr}({({\mathbf{\Sigma }}_t)}^2)+4(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_t^2(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})}}),\\ P(T_{4,t}^{*}\le x)\approx \Phi (\frac{\text{tr}({\mathbf{\Sigma }}_l^{*})-\text{tr}({\mathbf{\Sigma }}_t)-(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_t(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})+x\sqrt{2\text{tr}({\mathbf{\Sigma }}_l^{*\mathrm{ }2})}}{\sqrt{2\text{tr}({({\mathbf{\Sigma }}_t)}^2)+4(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_t^2(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})}}),\\ P(T_{6,t}^{*}\le x)\approx \Phi (\frac{\text{tr}({\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}({\mathbf{\Sigma }}_t^{*}-{\mathbf{\Sigma }}_t))-(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})+x\sqrt{2\text{tr}({({\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t^{*})}^2)}}{\sqrt{2\text{tr}({({\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t)}^2)+4(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t{\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})}}),\\ P(T_{7,t}^{*}\le x)\approx \Phi (\frac{\text{tr}({\mathbf{\Sigma }}_{d;l}^{*\mathrm{ }-1}{\mathbf{\Sigma }}_l^{*}-{\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t)-(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})+x\sqrt{2\text{tr}({({\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t^{*})}^2)}}{\sqrt{2\text{tr}({({\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t)}^2)+4(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t{\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})}}),\\ P(T_{8,t}^{*}\le x)\approx \Phi (\frac{\text{tr}({\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}({\mathbf{\Sigma }}_t^{*}-{\mathbf{\Sigma }}_t))-(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})+x\sqrt{2\text{tr}({({\mathbf{\Sigma }}_{d;l}^{*\mathrm{ }-1}{\mathbf{\Sigma }}_l^{*})}^2)}}{\sqrt{2\text{tr}({({\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t)}^2)+4(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t{\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})}}),\\ P(T_{9,t}^{*}\le x)\approx \Phi (\frac{\text{tr}({\mathbf{\Sigma }}_{d;l}^{*\mathrm{ }-1}{\mathbf{\Sigma }}_l^{*}-{\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t)-(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})+x\sqrt{2\text{tr}({({\mathbf{\Sigma }}_{d;l}^{*\mathrm{ }-1}{\mathbf{\Sigma }}_l^{*})}^2)}}{\sqrt{2\text{tr}({({\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t)}^2)+4(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t{\mathbf{\Sigma }}_{d;t}^{*\mathrm{ }-1}(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})}}),\\ P(T_{\mathit{\text{Mah}},t}^{*}\le x)\approx \Phi (\frac{p-\text{tr}({\mathbf{\Sigma }}_t^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t)-(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_t^{*\mathrm{ }-1}(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})+x\sqrt{2p}}{\sqrt{2\text{tr}({({\mathbf{\Sigma }}_t^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t)}^2)+4(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_t^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t{\mathbf{\Sigma }}_t^{*\mathrm{ }-1}(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})}}),\\ P(T_{\mathit{\text{MahInf}},t}^{*}\le x)\approx \Phi (\frac{\text{tr}({\mathbf{\Sigma }}_l^{*\mathrm{ }-1}({\mathbf{\Sigma }}_t^{*}-{\mathbf{\Sigma }}_t))-(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_l^{*\mathrm{ }-1}(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})+x\sqrt{2\text{tr}({({\mathbf{\Sigma }}_l^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t^{*})}^2)}}{\sqrt{2\text{tr}({({\mathbf{\Sigma }}_l^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t)}^2)+4(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})\prime {\mathbf{\Sigma }}_l^{*\mathrm{ }-1}{\mathbf{\Sigma }}_t{\mathbf{\Sigma }}_l^{*\mathrm{ }-1}(\mathit{\mu }+{\mathbf{a}}_{t-\tau }-{\mathit{\mu }}^{*})}}).\end{array}$$

We conclude this section with the analysis of the above asymptotic distributions obtained under the misspecification, which are compared with the corresponding high-dimensional distribution obtained without the effect of misspecification; that is, when ${\mathit{\mu }}^{*}$ and $\mathbf{\Gamma }{(h)}^{*}$ coincide with the true population values $\mathit{\mu }$ and $\mathbf{\Gamma }(h).$ These results provide initial intuition about the effect of misspecification. More detailed analysis is obtained via a Monte Carlo study and is presented in Section 5.

In the setup of the simulation study, we use the same data-generating model as described in Section 5.1. Moreover, to investigate the effects of the sample size, we set $N_0\in \{100,250,500\}$ when $p\in \{20,50\}$ and $N_0\in \{2000,5000,10000\}$ when $p\in \{500,1000\}.$ The results for the four variants of the MEWMA control chart based on the Euclidean distance and diagonalized Euclidean distance are depicted in Figures 1 to 4, and the findings obtained for the MEWMA approaches based on the exact and asymptotic Mahalanobis distance are summarized in Figure 5.

In the figures, we observe that the misspecification has only a minor effect on the in-control performance of the MEWMA control charts that are based on the Euclidean distance and diagonalized Euclidean distance in the case of moderate dimensions of the data-generating model (see Figures 1 and 2). All four curves in each plot are located very close to the density of the standard normal distribution, which is the limiting distribution of all considered control statistics when no misspecification is present; that is, the parameters of the in-control process are known. Some minor departures are present only when the sample size in Phase I is equal to 100. Interestingly, when p = 20, minor shifts of the densities to the right take place, whereas the asymptotic distributions under misspecification appear to have smaller variances when p = 50.

Figure 1. Probabilities $P(T_{1,t}^{*}\le x),P(T_{2,t}^{*}\le x),P(T_{3,t}^{*}\le x),P(T_{4,t}^{*}\le x),P(T_{6,t}^{*}\le x),P(T_{7,t}^{*}\le x),P(T_{8,t}^{*}\le x),$ and $P(T_{9,t}^{*}\le x)$ as functions of x for t = 5, p = 20, and $N_0\in \{100,250,500\}.$ The red plot corresponds to the density of the distribution in the case without misspecification; that is, the standard normal distribution.

Figure 2. Probabilities $P(T_{1,t}^{*}\le x),P(T_{2,t}^{*}\le x),P(T_{3,t}^{*}\le x),P(T_{4,t}^{*}\le x),P(T_{6,t}^{*}\le x),P(T_{7,t}^{*}\le x),P(T_{8,t}^{*}\le x),$ and $P(T_{9,t}^{*}\le x)$ as functions of x for t = 5, p = 50, and $N_0\in \{100,250,500\}.$ The red plot corresponds to the density of the distribution in the case without misspecification; that is, the standard normal distribution.

Figures 3 and 4 present the results for large dimensions of the data-generating model. In this case, the empirical densities are shifted to the left, indicating a more conservative behavior of the MEWMA control chart based on the Euclidean distance and diagonalized Euclidean distance. As such, the effect of misspecification is expected to result in larger values of the in-control ARLs for these control schemes.

Figure 3. Probabilities $P(T_{1,t}^{*}\le x),P(T_{2,t}^{*}\le x),P(T_{3,t}^{*}\le x),P(T_{4,t}^{*}\le x),P(T_{6,t}^{*}\le x),P(T_{7,t}^{*}\le x),P(T_{8,t}^{*}\le x),$ and $P(T_{9,t}^{*}\le x)$ as functions of x for t = 5, p = 500, and $N_0\in \{2000,5000,10000\}.$ The red plot corresponds to the density of the distribution in the case without misspecification; that is, the standard normal distribution.

Figure 4. Probabilities $P(T_{1,t}^{*}\le x),P(T_{2,t}^{*}\le x),P(T_{3,t}^{*}\le x),P(T_{4,t}^{*}\le x),P(T_{6,t}^{*}\le x),P(T_{7,t}^{*}\le x),P(T_{8,t}^{*}\le x),$ and $P(T_{9,t}^{*}\le x)$ as functions of x for t = 5, p = 1,000, and $N_0\in \{2000,5000,10000\}.$ The red plot corresponds to the density of the distribution in the case without misspecification; that is, the standard normal distribution.

The situation is completely different in Figure 5, where the results for the control charts based on the Mahalanobis distance are provided. The impact of the misspecification is considerable for this type of the MEWMA control scheme, even for moderate dimensions of the data-generating model. The blue curves in the plots, which correspond to the case of $N_0=100,$ are drastically shifted to the right. Moreover, the increase in the sample size does not obviously lead to the desired behavior of the asymptotic distributions. Even in the case of $N_0=500$ and p = 50, the asymptotic densities deviate significantly from the distribution corresponding to the case without misspecification. For larger dimensions of the data-generating model, the empirical densities are completely shifted to the right. This behavior would result in high probabilities of false alarms when the MEWMA control schemes based on the Mahalanobis distance are used. Such an effect can be explained by the large amount of noise that is present in the nondiagonal elements of the inverse estimated covariance matrices ${\mathbf{\Sigma }}_t^{*}$ and ${\mathbf{\Sigma }}_l^{*}.$ In contrast to the control schemes based on the Mahalanobis distance, the statistics of MEWMA control charts based on the Euclidian distance and the diagonalized Euclidean distance do not use the nondiagonal elements of ${\mathbf{\Sigma }}_t^{*}$ and ${\mathbf{\Sigma }}_l^{*}$ in their definitions, and, as such, they appear to be robust against misspecification.

Figure 5. Probabilities $P(T_{\mathit{\text{Mah}},t}^{*}\le x)$ (left column) and $P(T_{\mathit{\text{MahInf}},t}^{*}\le x)$ (right column) as functions of x for t = 5. We set $p\in \{20,50\}$ with $N_0\in \{100,250,500\}$ and $p\in \{500,1000\}$ with $N_0\in \{2000,5000,10000\}.$ The red plot corresponds to the density of the distribution in the case without misspecification; that is, the standard normal distribution.

5. COMPARISON STUDY

The above results characterize the behavior of the considered MEWMA charts if the parameters of the process are misspecified. Thus, we can compare the charts and see how the charts react to deviations from the true parameters. This point will be illustrated in Section 5.2.

In practice, we estimate the process parameters within the Phase I analysis. Usually a preliminary run or historical data is used. Several approaches have been proposed in the literature, which we discussed in Section 3. If we apply the control statistics $T_{i,t}=T_{i,t}(\mathit{\theta })$ proposed above, then we have to estimate the parameter by the prerun and get ${\widehat{T}}_{i,t}=T_{i,t}(\widehat{\mathit{\theta }}).$ Here we will assume that the estimators obtained in the Phase I analysis are independent from the control statistics used in the Phase II analysis. With the law of total probability, we get (4.7) $$P(T_{i,t}(\widehat{\mathit{\theta }})\le x)={\int }_{}^{}P(T_{i,t}(\widehat{\mathit{\theta }})\le x|\widehat{\mathit{\theta }}={\mathit{\theta }}^{*})f_{\widehat{\mathit{\theta }}}({\mathit{\theta }}^{*})d{\mathit{\theta }}^{*}={\int }_{}^{}P(T_{i,t}({\mathit{\theta }}^{*})\le x)f_{\widehat{\mathit{\theta }}}({\mathit{\theta }}^{*})d{\mathit{\theta }}^{*}.$$(4.7)

We have exclusively analyzed $P(T_{i,t}({\mathit{\theta }}^{*})\le x)$ in the previous sections. To study (4.7)(4.7) $$P(T_{i,t}(\widehat{\mathit{\theta }})\le x)={\int }_{}^{}P(T_{i,t}(\widehat{\mathit{\theta }})\le x|\widehat{\mathit{\theta }}={\mathit{\theta }}^{*})f_{\widehat{\mathit{\theta }}}({\mathit{\theta }}^{*})d{\mathit{\theta }}^{*}={\int }_{}^{}P(T_{i,t}({\mathit{\theta }}^{*})\le x)f_{\widehat{\mathit{\theta }}}({\mathit{\theta }}^{*})d{\mathit{\theta }}^{*}.$$(4.7) , we have to know the distribution of the estimator $\widehat{\mathit{\theta }}.$ This distribution is unknown; we only know that these quantities are usually asymptotically normally distributed. This fact could be used to get an approximation to (4.7)(4.7) $$P(T_{i,t}(\widehat{\mathit{\theta }})\le x)={\int }_{}^{}P(T_{i,t}(\widehat{\mathit{\theta }})\le x|\widehat{\mathit{\theta }}={\mathit{\theta }}^{*})f_{\widehat{\mathit{\theta }}}({\mathit{\theta }}^{*})d{\mathit{\theta }}^{*}={\int }_{}^{}P(T_{i,t}({\mathit{\theta }}^{*})\le x)f_{\widehat{\mathit{\theta }}}({\mathit{\theta }}^{*})d{\mathit{\theta }}^{*}.$$(4.7) , but then the problem is to determine the resulting integral. Here we choose another procedure to evaluate (4.7)(4.7) $$P(T_{i,t}(\widehat{\mathit{\theta }})\le x)={\int }_{}^{}P(T_{i,t}(\widehat{\mathit{\theta }})\le x|\widehat{\mathit{\theta }}={\mathit{\theta }}^{*})f_{\widehat{\mathit{\theta }}}({\mathit{\theta }}^{*})d{\mathit{\theta }}^{*}={\int }_{}^{}P(T_{i,t}({\mathit{\theta }}^{*})\le x)f_{\widehat{\mathit{\theta }}}({\mathit{\theta }}^{*})d{\mathit{\theta }}^{*}.$$(4.7) and make use of simulations. Our procedure and the results are provided in Section 5.2.

5.1. Design of the Comparison Study

In Section 5.1 we described the design of our simulation study. We have chosen the in-control process to be a VAR(1) process; that is, ${\mathbf{Y}}_t=\mathbf{\Phi }{\mathbf{Y}}_{t-1}+{\mathit{\epsilon }}_t$ with independent white noise process $\left\{{\mathit{\epsilon }}_t\right\}$ and

1. $\mathbf{\Phi }=\phi \mathbf{I}$ with $\phi =0.5,$

2. ${\mathit{\epsilon }}_t\sim {\mathcal{N}}_p(\mathbf{0},\mathbf{\Sigma })$ with $\mathbf{\Sigma }=\mathbf{\text{DAD}},$

where $\mathbf{D}=\mathit{\text{diag}}(d_1,\dots ,d_p)$ is a diagonal matrix consisting of the standard deviations $d_1,\dots ,d_p$ and $$\mathbf{A}=\left(\begin{array}{ccccc}1& \alpha & {\alpha }^2& \dots & {\alpha }^{p-1}\\ \alpha & 1& \alpha & \dots & {\alpha }^{p-2}\\ {\alpha }^2& \alpha & 1& \dots & {\alpha }^{p-3}\\ & ⋮& ⋮& \ddots & ⋮\\ {\alpha }^{p-1}& {\alpha }^{p-2}& {\alpha }^{p-3}& \dots & 1\end{array}\right)\hspace{1em}$$

is a correlation matrix with $\alpha =0.5.$ The values of $d_1,\dots ,d_p$ are drawn randomly from the uniform distribution on the interval $[0.5,2].$

Additionally, we set $\mathbf{R}=r\mathbf{I}$ with $r\in \{0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1\}$ and the dimension of the observed process and the target process p to be equal to 20. Note that at the beginning we determine the matrix $\mathbf{\Sigma }$ and keep it constant for all of our simulations. Furthermore, we will always assume in this article that $\mathit{\mu }=\mathbf{0}$ is known and thus we do not estimate $\mathit{\mu };$ that is, $\mathit{\mu }={\mathit{\mu }}^{*}=\mathbf{0}.$

5.2. Behavior of the Control Statistics under Misspecification

The control procedures will be compared using the ARL as a performance criterion (see, e.g., Montgomery 2020). This is the most popular performance measure in SPC. Whereas the run length is equal to the first value at which the control chart gives a signal, the ARL is equal to the expected run length assuming that the change has occurred at the first position; that is, τ = 1.

Because the control statistics are dependent variables over time having a complicate correlation structure, no explicit expressions for the ARL are known. Note that even for univariate time series, explicit expressions are only known for special cases (see, e.g., Schmid 1995). In the case of an independent sequence of normally distributed variables, the ARL is obtained by solving an integral equation (e.g., Knoth 2021). This is why we use simulations to estimate the ARL.

In a first step, we generate a preliminary sample of the above VAR(1) process of size 500 and use it to estimate the parameters of the VAR(1) process following the three estimators presented in Section 3. It is assumed that the preliminary sample is in control. Then, in the Phase II analysis, we generate a realization of the Phase II process. Here we restrict ourselves to the case where the Phase II process is in control as well. The Phase I process and the Phase II process are assumed to be independent of each other. Now we apply one of the above-discussed control charts, fix a control limit, and determine the run length for the given data. This procedure is repeated 1,000 times and the ARL is used as an estimator of the true run length. Consequently, we get an estimator of the in-control ARL for a given control limit.

To compare the control charts, they have to be calibrated. This means that the control limit for each chart is determined in such way that the corresponding in-control ARL is equal to a fixed value ξ. In this article, we chose ξ = 200. This means that, on average, after 200 observations a wrong decision is made; that is, it is concluded that the process is out of control but in reality it is in control. The choice of ξ depends on the data frequency. In engineering, frequently the value 500 is chosen, whereas in finance smaller values are taken; for example, ξ = 60. A discussion on the choice of ξ can be found in Severin and Schmid (1999).

Now we have to determine the control limits for each chart. Therefore, the Regula falsi is applied to the estimated in-control ARL; that is, the estimated in-control ARL is calculated for various values of control limits based on 10⁴ independent observations. This procedure leads to estimators for the control limits. R. Bodnar, Bodnar, and Schmid (2023) noticed that the control charts $T_{6,t}$ and $T_{7,t}$ have the same behavior as $T_{8,t}$ and $T_{9,t},$ respectively. Thus, in the following we will consider only $T_{6,t}$ and $T_{7,t}.$ Moreover, we will omit the chart $T_{2,t}$ because even for the true process this control scheme performs very poorly in both the in-control and out-of-control states.

Tables 1 and 2 depict the in-control ARLs computed for the misspecified $T_{1,t}^{*},T_{3,t}^{*},T_{4,t}^{*},T_{6,t}^{*},T_{7,t}^{*},T_{\mathit{\text{Mah}},t}^{*},$ and $T_{\mathit{\text{MahInf}},t}^{*}$ MEWMA control charts for various values of $r\in \{0.1,0.2,\dots ,1.0\}$ when the in-control process is the 20-dimensional VAR(1) model in Table 1 and the 100-dimensional VAR(1) model in Table 2 as defined in Section 5.1. We observe that the estimation method has only a minor impact on the computed values of the estimated ARLs with the exception present for both MEWMA control schemes based on the Mahalanobis distance when r is small. In general, quite robust behavior toward the parameter misspecification is present when the control statistics are defined by using the Euclidean norm and the diagonalized Euclidean norm. In all of the cases, the computed ARLs depart from the target ARL by not more than 15% when p = 20. In the case of a larger dimension of the data-generating model, the control charts based on the Euclidean norm and the diagonalized Euclidean norm become more conservative with the deviation of the empirical ARLs to the target value of 200 bounded by around 65%. The results obtained for the control schemes based on the Mahalanobis distance are completely different. The ARLs are considerably smaller than the target value of 200, especially, when r is small. If p = 100, then the computed ARLs for the MEWMA control charts based on the Mahalanobis distance are very small, indicating high probabilities of the false alarm in this case due to the misspecification. These results are in line with our previous findings depicted in Figures 1 to 5, where it was noted that the high-dimensional asymptotic distributions of the control statistics in the case of Mahalanobis distance are drastically moved to the right, which results in small values of the estimated ARLs. Because the control statistics, which are based on the Euclidean distance and the diagonalized Euclidean distance, do not depend on the inverse of ${\mathbf{\Sigma }}_t^{*}$ and ${\mathbf{\Sigma }}_l^{*},$ these control charts are robust to large estimation errors that are present when the inverse of a covariance matrix is estimated. Finally, a similar performance of the MEWMA control charts based on the Euclidean distance and the diagonalized Euclidean distance is displayed in Table 2, which is again in line with the results depicted in Figures 3 and 4.

Table 1. ARLs of the $T_{1,t}^{*},T_{3,t}^{*},T_{4,t}^{*},T_{6,t}^{*},T_{7,t}^{*},T_{\mathit{\text{Mah}},t}^{*},$ and $T_{\mathit{\text{MahInf}},t}^{*}$ MEWMA control charts for $r\in \{0.1,0.2,\dots ,1.0\}$ when the in-control process is the 20-dimensional VAR(1) process.

r	Estimator	$T_{1,t}^{*}$	$T_{3,t}^{*}$	$T_{4,t}^{*}$	$T_{6,t}^{*}$	$T_{7,t}^{*}$	$T_{\mathit{\text{Mah}},t}^{*}$	$T_{\mathit{\text{MahInf}},t}^{*}$
	Nonparametric	223.0	218.6	220.4	201.5	201.5	53.7	56.8
0.1	VAR-based nonparametric	220.0	217.8	232.6	210.5	210.5	35.2	43.0
	Maximum likelihood	232.5	216.2	227.0	202.7	202.7	34.9	43.4
	Nonparametric	225.9	222.1	222.2	206.9	206.9	48.8	49.8
0.2	VAR-based nonparametric	207.7	225.0	225.3	204.8	204.8	35.2	39.4
	Maximum likelihood	225.6	225.9	224.0	203.5	203.5	36.7	39.9
	Nonparametric	217.4	222.9	221.3	207.3	207.3	50.5	49.9
0.3	VAR-based nonparametric	220.0	220.2	207.9	202.9	202.9	38.6	41.2
	Maximum likelihood	220.8	224.2	229.0	220.2	220.2	38.7	41.0
	Nonparametric	218.7	219.9	218.1	205.5	205.5	54.0	53.5
0.4	VAR-based nonparametric	227.4	211.2	227.9	183.5	183.5	40.9	42.7
	Maximum likelihood	220.5	217.5	218.8	197.9	197.9	46.1	44.7
	Nonparametric	217.4	217.9	217.4	200.9	200.9	56.5	57.0
0.5	VAR-based nonparametric	219.1	214.3	224.0	203.9	203.9	47.1	46.6
	Maximum likelihood	202.2	217.5	228.1	206.9	206.9	46.8	49.3
	Nonparametric	219.0	217.6	219.0	205.4	205.4	60.6	61.1
0.6	VAR-based nonparametric	220.5	195.7	208.2	196.4	202.2	52.7	53.5
	Maximum likelihood	226.3	219.9	214.6	201.7	201.7	51.3	54.3
	Nonparametric	220.0	220.9	220.0	207.9	207.9	66.7	65.7
0.7	VAR-based nonparametric	218.8	225.0	207.6	198.6	191.7	57.9	58.5
	Maximum likelihood	191.8	227.0	215.9	195.8	195.8	57.5	61.4
	Nonparametric	221.7	221.2	221.7	205.2	205.2	70.6	72.0
0.8	VAR-based nonparametric	225.3	221.5	219.0	200.8	206.9	63.5	62.6
	Maximum likelihood	214.4	209.5	213.8	207.8	207.8	65.7	67.7
	Nonparametric	216.8	221.6	216.8	202.9	202.9	75.3	75.5
0.9	VAR-based nonparametric	238.7	223.8	221.3	203.3	208.2	70.8	71.1
	Maximum likelihood	217.6	222.6	195.0	203.8	203.8	71.9	70.4
	Nonparametric	220.6	220.6	220.6	200.8	200.8	76.8	76.8
1.0	VAR-based nonparametric	221.4	221.4	221.4	199.3	200.6	81.8	76.8
	Maximum likelihood	213.8	213.8	213.8	201.3	201.3	78.0	78.0

The estimation in Phase 1 is performed by the three approaches described in Section 3.

Table 2. ARLs of the $T_{1,t}^{*},T_{3,t}^{*},T_{4,t}^{*},T_{6,t}^{*},T_{7,t}^{*},T_{\mathit{\text{Mah}},t}^{*},$ and $T_{\mathit{\text{MahInf}},t}^{*}$ MEWMA control charts for $r\in \{0.1,0.2,\dots ,1.0\}$ when the in-control process is the 100-dimensional VAR(1) process.

r	Estimator	$T_{1,t}^{*}$	$T_{3,t}^{*}$	$T_{4,t}^{*}$	$T_{6,t}^{*}$	$T_{7,t}^{*}$	$T_{\mathit{\text{Mah}},t}^{*}$	$T_{\mathit{\text{MahInf}},t}^{*}$
	Nonparametric	247.2	301.3	299.3	282.0	282.0	3.2	5.6
0.1	VAR-based nonparametric	328.9	315.2	213.5	302.6	302.6	1.4	4.2
	Maximum likelihood	334.7	327.0	215.7	315.4	315.4	1.4	4.2
	Nonparametric	311.4	310.1	315.9	297.1	297.1	2.9	3.8
0.2	VAR-based nonparametric	328.3	308.6	230.7	302.2	302.2	1.5	3.1
	Maximum likelihood	317.0	305.2	226.0	312.2	312.2	1.4	3.1
	Nonparametric	293.4	300.3	305.4	294.2	294.2	2.7	3.2
0.3	VAR-based nonparametric	309.4	309.5	246.1	310.6	310.6	1.5	2.7
	Maximum likelihood	328.4	303.3	245.0	315.7	315.7	1.5	2.7
	Nonparametric	298.1	297.3	300.1	295.5	295.5	2.6	2.8
0.4	VAR-based nonparametric	313.5	292.9	272.1	306.8	306.8	1.5	2.4
	Maximum likelihood	322.9	297.1	272.6	309.8	309.8	1.5	2.4
	Nonparametric	301.3	299.3	297.9	295.1	295.1	2.5	2.6
0.5	VAR-based nonparametric	275.8	288.0	292.2	292.0	292.0	1.6	2.3
	Maximum likelihood	297.0	297.3	289.4	303.7	303.7	1.5	2.3
	Nonparametric	286.5	287.2	274.2	271.4	271.4	2.5	2.5
0.6	VAR-based nonparametric	308.4	276.6	288.1	295.7	295.7	1.6	2.3
	Maximum likelihood	303.5	290.2	289.9	289.9	289.9	1.6	2.3
	Nonparametric	281.8	283.9	285.8	281.6	281.6	2.4	2.4
0.7	VAR-based nonparametric	291.2	288.7	239.7	243.2	243.2	1.6	2.2
	Maximum likelihood	298.3	264.7	235.1	258.0	258.0	1.6	2.2
	Nonparametric	281.4	285.6	283.2	277.3	277.3	2.3	2.3
0.8	VAR-based nonparametric	301.4	278.2	282.1	288.7	288.7	1.7	2.2
	Maximum likelihood	286.7	282.3	279.3	254.8	254.8	1.7	2.1
	Nonparametric	268.6	270.7	272.7	278.7	278.7	2.1	2.1
0.9	VAR-based nonparametric	289.3	283.5	278.0	272.6	272.6	1.8	2.0
	Maximum likelihood	275.8	270.9	269.4	274.5	274.5	1.8	2.0
	Nonparametric	272.1	272.1	272.1	275.2	275.2	2.0	2.0
1.0	VAR-based nonparametric	272.1	272.1	272.1	275.2	275.2	1.9	2.0
	Maximum likelihood	276.0	273.6	273.6	266.7	266.7	1.9	1.9

The estimation in Phase 1 is performed by the three approaches described in Section 3.

6. CONCLUSION

The commonly used multivariate control charts are derived under the assumption that the parameters of the target process are known before the control procedure starts. However, this assumption is not fulfilled in many situations of practical interest; for example, in economics and finance. In such situations, a detailed analysis is performed in Phase 1 of the surveillance procedure with the aim to fit the target process by estimating its parameters. The quality of the estimator is expected to have a considerable impact on the control procedure, especially when the data-generating model is of large dimension.

In this article, we analyze the in-control properties of the MEWMA control charts whose control statistics are based on the Mahalanobis distance, the Euclidean distance, and the diagonalized Euclidean distance. Misspecified MEWMA control schemes are proposed in which the unknown parameters of the in-control process are replaced by their estimators. The distributional properties of the control statistics of the misspecified control charts are investigated and their high-dimensional asymptotic distributions are derived. The established theoretical results are used to study the effect of the mispecification, and the finite-sample behavior of the control statistics is assessed via an intensive simulation study. Based on our findings, it is concluded that the MEWMA control schemes based on the Mahalanobis distance can suffer from considerable misspecification, whereas the control charts whose test statistics are defined by using the Euclidean distance and the diagonalized Euclidean distance are quite robust to misspecification, because we estimate the parameters of the target process in the Phase I analysis of the monitoring process.

DISCLOSURE

The authors have no conflicts of interest to report.

ACKNOWLEDGEMENT

The authors thank the Editor, the Associate Editor, and two anonymous reviewers for their constructive comments that improved the quality of this article.

Additional information

Funding

The first and third authors acknowledge financial support from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) Project No. 428472210.