Robust online-surveillance of trend-coherence in multivariate data streams: the similar trend monitoring (STM) procedure

Abstract

When several data streams are observed simultaneously, it is often of great interest to monitor the coherences between all pairs of streams. We propose a new technique called Similar Trend Monitoring (STM) for this task: The current slopes of all univariate streams are estimated and compared pairwise at each time point. The STM statistic is the standardized slope difference, so that decisions about coherence can be made by means of the six-sigma-rule, for instance. The STM meets the high demands that come along with the online monitoring of multivariate data streams: it is fast to compute, robust against outliers, applicable when observations are missing, and does not require stationarity of the processes. We investigate the distribution and the performance of the STM and demonstrate its capabilities considering blood pressure time series from intensive care patient monitoring as an example.

Appendices

Appendix 1: The Orthogonalized Gnanadesikan- Kettenring covariance estimator based on $\mathrm{Q^\mathrm{adj}}$

We use the Orthogonalized Gnanadesikan-Kettenring covariance estimator (OGK, Maronna and Zamar 2002; Gnanadesikan and Kettenring 1972) to estimate the local error cross-covariance matrix $\varvec{\varSigma }_t$, where $\mathrm{Cov}[\varepsilon _t(i),\varepsilon _t(j)] = \sigma _t(i,j)$ and $\mathrm{Corr}[\varepsilon _t(i),\varepsilon _t(j)] = \sigma _t(i,j) / [\sigma _t(i) \cdot \sigma _t(j)]$. Let $\varvec{X}\in \mathbb {R}^{n \times K}$ denote a sample matrix with rows $\varvec{x}_i = (x_{i1},$ $\ldots ,$ $x_{in})$ and columns $X_k = (x_{1k},\ldots ,x_{nk})'$ with $i=1,\ldots ,n$ and $k=1,\ldots ,K$. Given a univariate scale estimator $\hat{\sigma }()$, the covariance estimation $\mathrm{OGK}(\varvec{X})$ is carried out as follows:

• 1. Determine $\varvec{D} = \mathrm{diag}(\hat{\sigma }(X_1),\ldots ,\hat{\sigma }(X_K))$, so that $\varvec{Y}= \varvec{X}\varvec{D}^{-1}$.

• 2. Estimate the correlation matrix $\varvec{R}$ of $\varvec{X}$:

  $$\begin{aligned} R_{ij} = \left( \hat{\sigma }(Y_i + Y_j)^2 - \hat{\sigma }(Y_i - Y_j)^2 \right) \bigg /4, \end{aligned}$$

  where $R_{ii}=1$, and $Y_i$ denotes the $i$-th column of $\varvec{Y}$.

• 3. Conduct an eigenvalue decomposition $\varvec{R} = \varvec{E}\varvec{\varLambda }\varvec{E}'$, where $\varvec{\varLambda }$ $=$ $\mathrm{diag}(\lambda _1,\ldots ,\lambda _K)$ contains the ordered eigenvalues and $\varvec{E}$ the referring eigenvectors of $\varvec{R}$.

• 4. Obtain the OGK covariance estimate by $\mathrm{OGK}(\varvec{X}) = \varvec{A}\varvec{\varGamma }\varvec{A}'$ with $\varvec{A}=\varvec{DE}$ and $\varvec{\varGamma } = \mathrm{diag}(\hat{\sigma }(Z_1)^2,\ldots ,\hat{\sigma }(Z_K)^2)$, where $Z_i$ denotes the $i$-th column of $\varvec{Z}= \varvec{X}(\varvec{A}')^{-1}$.

In order to estimate the error covariance matrix $\varvec{\varSigma }_t$, the obvious approach is to apply this algorithm to a matrix $\varvec{X}$ of residuals, which can be residuals of RM regression lines, for instance, that are fitted to univariate time windows. Lanius and Gather (2010) follow this approach and suggest to apply the $Q_n$ scale estimator (Rousseeuw and Croux 1993) to RM residuals. However, we do not use the OGK based on the $Q_n$ but on the $\mathrm{Q^\mathrm{adj}}$ scale estimator instead. One reason is that the $\mathrm{Q^\mathrm{adj}}$ estimator is model-free, i.e. it does not assume a certain shape of the signal. Therefore, the OGK-$\mathrm{Q^\mathrm{adj}}$ delivers good results when linearity is not given, in contrast to the OGK-${Q_n}$ (Chapter 4.2.2 Busse 2012). Furthermore, the $\mathrm{Q^\mathrm{adj}}$ estimations $\hat{\sigma }_t(k) = \hat{\sigma }(X_k)$, $k=1,\ldots ,K$, are delivered by the SCARM anyway, so that no additional computing time is needed to determine the matrix $\varvec{D}$ at step 1 of the OGK algorithm.

The OGK-$\mathrm{Q^\mathrm{adj}}$ estimates the error covariance directly from the data, i.e. the matrix $\varvec{X}$ in the OGK-$\mathrm{Q^\mathrm{adj}}$-algorithm contains the $n$ most recent observation vectors. (In contrast, the matrix $\varvec{X}$ in the OGK-$Q_n$-algorithm contains residuals.) However, we have to decide how many observations are used for the OGK-$\mathrm{Q^\mathrm{adj}}$ estimation. The STM uses the SCARM to obtain an individual window width $n_t(k)$ for all of the $K$ univariate components of the $K$-variate time series. That is, we do not have a sample matrix $\varvec{X}\in \mathbb {R}^{K \times n}$ but $K$ sample vectors $X_k$ of different lengths $n_t(k)$. We therefore estimate $\sigma _t(i,j)$ from the matrix $\varvec{X}^* = (X^*_i,X^*_j) \in \mathbb {R}^{n \times 2}$, where $X^*_i$ and $X^*_j$ are the recent $n = \max \ \{ n_t(i), n_t(j) \}$ observations of component $i$ and $j$. Choosing $n$ as the minimum of $n_t(i)$ and $n_t(j)$ would be the more comprehensible approach at first glance, since both $X^*_i$ and $X^*_j$ would have a linear structure then. However, this would mean discarding observations, and linearity is not required for the $\mathrm{Q^\mathrm{adj}}$ estimator.

Note that the OGK is not affine equivariant, but it is approximately affine equivariant, see Maronna and Zamar (2002). Furthermore, given a sample matrix $\varvec{X}=(X_1,X_2) \in \mathbb {R}^{n \times 2}$, the OGK is regression invariant if the same trend is added to both univariate samples $X_2$ and $X_2$ and a regression invariant scale estimator $\hat{\sigma }()$ is used:

Lemma 1

Let $\varvec{X}{=}(X_1,X_2) {\in } \mathbb {R}^{n \times 2}$ and $\varvec{X}^* {=} (X_1^*,X_2^*) \in \mathbb {R}^{n \times 2}$, $k=1,2$, where $X_k$ = $(x_{ik})_{i=1,\ldots ,n}'$ and $X_k^*$ = $(x_{ik}^*)_{i=1,\ldots ,n}'$ with $x_{ik}^* = x_{ik} + i \cdot a + b$, $a,b \in \mathbb {R}$, $i=1,\ldots ,n$. If the OGK uses a regression invariant scale estimator $\hat{\sigma }()$, then OGK$(\varvec{X})$ = OGK$(\varvec{X}^*)$.

Proof

In the following, $\varvec{D}, \varvec{Y}, \varvec{R}, \varvec{E}, \varvec{A}, \varvec{Z}, \varvec{\varGamma }$ denote the matrices of the OGK-algorithm from above. Obviously,

$$\begin{aligned} \varvec{D}= \mathrm{diag}\left( \hat{\sigma }(\varvec{X}_1), \hat{\sigma }(\varvec{X}_2) \right) = \mathrm{diag}\left( \hat{\sigma }(\varvec{X}_1^*), \hat{\sigma }(\varvec{X}_2^*) \right) =: \varvec{D}^*, \end{aligned}$$

since $\hat{\sigma }()$ is regression invariant. Let $\varvec{R}^* = (R_{ij})_{i,j=1,2}$ with $R^*_{ii}=1$ and $R_{ij}^* = \left( \hat{\sigma }(Y_1^* + Y_2^*)^2 - \hat{\sigma }(Y_1^* - Y_2^*)^2 \right) / 4$. Let $\varvec{Y}^* = \varvec{X}^* \varvec{D}^{-1}$ with $\varvec{Y}^* = (Y_1^*, Y_2^*) \in \mathbb {R}^{n \times 2}$. It is not difficult to show that

$$\begin{aligned} \hat{\sigma }(Y_1^* + Y_2^*) = \hat{\sigma }(Y_1 + Y_2) ~\hbox {and}~ \hat{\sigma }(Y_1^* - Y_2^*) = \hat{\sigma }(Y_1 - Y_2), \end{aligned}$$

since $\hat{\sigma }(Y_1^*) = \hat{\sigma }(Y_1)$ and $\hat{\sigma }(Y_2^*) = \hat{\sigma }(Y_2)$ due to the regression invariance of $\hat{\sigma }()$. Hence, $R_{12} {=} R^*_{12}$ and obviously also $R_{21} {=} R^*_{21}$. Thus, $\varvec{E}{=} \varvec{E}^*$, where $\varvec{E}^*$ denotes the matrix $\varvec{E}$ computed on $\varvec{X}^*$ instead of $\varvec{X}$. Thus, $\varvec{A}= \varvec{D}\varvec{E}= \varvec{D}^* \varvec{E}^* =: \varvec{A}^*$. Furthermore, it is $\varvec{Z}= \varvec{X}(\varvec{A}')^{-1}$, see the OGK-algorithm, where $\varvec{Z}= (Z_1,Z_2) \in \mathbb {R}^{n \times 2}$. Let $\varvec{Z}^* := \varvec{X}^* ((\varvec{A}^*)')^{-1} = \varvec{X}^* (\varvec{A}')^{-1}$ with $\varvec{Z}^* = (Z_1^*,Z_2^*) \in \mathbb {R}^{n \times 2}$. It is not difficult to show that

$$\begin{aligned} \hat{\sigma }(Z_1^*) = \hat{\sigma }(Z_1) ~\hbox {and}~ \hat{\sigma }(Z_2^*) = \hat{\sigma }(Z_2), \end{aligned}$$

because of the regression invariance of $\hat{\sigma }()$. Hence $\varvec{\varGamma }^* := \mathrm{diag}\left( \hat{\sigma }(Z_1^*)^2, \hat{\sigma }(Z_2^*)^2 \right) = \varvec{\varGamma }$ and, thus,

$$\begin{aligned} \hbox {OGK}(\varvec{X}) = \varvec{A}\varvec{\varGamma }\varvec{A}' = \varvec{A}^* \varvec{\varGamma }^* (\varvec{A}^*)' = \hbox {OGK}(\varvec{X}^*). \end{aligned}$$

Appendix 2: Distribution of the STM statistic: further simulations

Here we investigate the distribution of the STM statistic given AR(1) time series from (16) with innovations $e_t$ following

• (i) a $t$-distribution with three degrees of freedom (heavy-tailed errors),

• (ii) a Weibull distribution with shape parameter two and scale parameter one (skewed errors),

• (iii) a $0.95N(0,1) + 0.05N(5,1)$-distribution (contaminated normal errors).

For each error type (i)–(iii), we consider each combination of $\varphi (1)$ and $\varphi (2)$ with $\varphi (1) \le \varphi (2)$ and $\varphi (1),\varphi (2) \in \{-0.9,-0.6,\ldots ,0.9\}$, and compute 100.000 STM statistics in each case. Histograms of the STM statistics can be found on the website mentioned in Sect. 2.5. The empirical quantiles computed from these statistics are listed in Tables 2, 3 and 4. The quantiles in Table 3 are quite similar to those in Table 1. The quantiles in Tables 2 and 4 are larger than the corresponding standard normal quantiles in almost all cases, in particular if the autocorrelations are strong. Furthermore, the quantiles in Tables 2 and 4 are larger than the quantiles in Table 1 if $\varphi (1)$ or $\varphi (2)$ is equal to $0.6$ or $0.9$. Obviously, the combination of extreme observations in the data stream together with strong positive autocorrelations leads to an increased probability for huge absolute values of the STM statistic.

Table 2 Empirical $(1-\gamma )$-quantiles of 100.000 STM statistics computed from two AR(1) time series with autocorrelation parameters $\varphi (k)$, $k=1,2$, and innovations from a $t_3$-distribution

Table 3 Empirical $(1-\gamma )$-quantiles of 100.000 STM statistics computed from two AR(1) time series with autocorrelation parameters $\varphi (k)$, $k=1,2$, and innovations from a Weibull distribution

Table 4 Empirical $(1-\gamma )$-quantiles of 100.000 STM statistics computed from two AR(1) time series with autocorrelation parameters $\varphi (k)$, $k=1,2$, and innovations from a contaminated normal distribution