Sequential Detection of Three-Dimensional Signals under Dependent Noise

TL;DR

This paper develops a sequential detection method for three-dimensional signals in space-time data affected by dependent noise, providing theoretical convergence results and extensive simulations.

Contribution

It introduces a novel sequential detection approach for 3D signals under dependent noise, with proven asymptotic properties and extensions for unknown reference signals.

Findings

Detector process converges to Brownian motion under null hypothesis

Limiting distribution derived for various local alternatives

Simulation study confirms effectiveness of the detection algorithm

Abstract

We study detection methods for multivariable signals under dependent noise. The main focus is on three-dimensional signals, i.e. on signals in the space-time domain. Examples for such signals are multifaceted. They include geographic and climatic data as well as image data, that are observed over a fixed time horizon. We assume that the signal is observed as a finite block of noisy samples whereby we are interested in detecting changes from a given reference signal. Our detector statistic is based on a sequential partial sum process, related to classical signal decomposition and reconstruction approaches applied to the sampled signal. We show that this detector process converges weakly under the no change null hypothesis that the signal coincides with the reference signal, provided that the spatial-temporal partial sum process associated to the random field of the noise terms disturbing the sampled signal con- verges to a Brownian motion. More generally, we also establish the limiting distribution under a wide class of local alternatives that allows for smooth as well as discontinuous changes. Our results also cover extensions to the case that the reference signal is unknown. We conclude with an extensive simulation study of the detection algorithm.

Figures

Tables

Table 1. Table 1: Simulated rejection rates for various sampling periods τ 1 , τ 2 subscript 𝜏 1 subscript 𝜏 2 \tau_{1},\tau_{2} and τ 3 subscript 𝜏 3 \tau_{3} for 1000 repetitions.

${\tau }_1$	$0.025$	$0.0\overline{3}$	$0.04$	$0.05$	$0.0\overline{6}$	$0.08$	$0.1$
${\tau }_2={\tau }_3=0.0\overline{3}$	0.074	0.064	0.067	0.079	0.076	0.070	0.055
${\tau }_2={\tau }_3=0.04$	0.074	0.066	0.067	0.061	0.064	0.076	0.062
${\tau }_2={\tau }_3=0.05$	0.074	0.069	0.084	0.070	0.063	0.066	0.058

Table 2. Table 2: Simulated rejection rates for various sampling periods τ 1 , τ 2 subscript 𝜏 1 subscript 𝜏 2 \tau_{1},\tau_{2} , and τ 3 subscript 𝜏 3 \tau_{3} for 10000 repetitions.

${\tau }_1$	$0.08$	$0.1$
${\tau }_2={\tau }_3=0.04$	0.0652	0.0598
${\tau }_2={\tau }_3=0.05$	0.0587	0.0567

Table 3. Table 3: Simulated rejection rates for 1000 repetitions for control limits calculated with the true σ 𝜎 \sigma ( n 1 = 50 , n 0 = 12 formulae-sequence subscript 𝑛 1 50 subscript 𝑛 0 12 n_{1}=50,n_{0}=12 ).

$\rho$	$0.1$	$0.2$	$0.3$	$0.4$	$0.5$
$a=-0.1$	0.099	0.086	0.065	0.053	0.036
$a=0.01$	0.063	0.063	0.059	0.056	0.039
$a=0.1$	0.058	0.053	0.047	0.041	0.035
$a=0.3$	0.033	0.033	0.036	0.038	0.033
$a=0.5$	0.031	0.028	0.027	0.029	0.028

Table 4. Table 4: Simulated power for the models in ( 5.1 ) for β 1 = 0.5 subscript 𝛽 1 0.5 \beta_{1}=0.5 and different sample sizes n 1 subscript 𝑛 1 n_{1} .

$n_1$	$80$	$160$	$240$	$320$
${\beta }_1=0.5$	0.689	0.760	0.778	0.786

Equations

$$y_{\mathbf{i}}=f\left(i\tau \right)+{\epsilon }_{\mathbf{i}},$$

$$y_{i_1,i_2,i_3}=f\left(i_1{\tau }_1,i_2{\tau }_2,i_3{\tau }_3\right)+{\epsilon }_{i_1,i_2,i_3}.$$

$$f(t,r_2,r_3)=$$

$$\text{sinc}\left({\Omega }_1\left(t-\frac{i_1\pi }{{\Omega }_1}\right)\right)\text{sinc}\left({\Omega }_2\left(r_2-\frac{i_2\pi }{{\Omega }_2}\right)\right)\text{sinc}\left({\Omega }_3\left(r_3-\frac{i_3\pi }{{\Omega }_3}\right)\right),$$

$$f(t,r_2,r_3)=$$

$$\text{sinc}(\pi {\tau }_1^{-1}(t-i_1{\tau }_1))\text{sinc}(\pi {\tau }_2^{-1}(r_2-i_2{\tau }_2))\text{sinc}(\pi {\tau }_3^{-1}(r_3-i_3{\tau }_3))$$

$$g(t,r_2,r_3)=\frac{{\Omega }_1{\Omega }_2{\Omega }_3}{{\pi }^3}\text{sinc}({\Omega }_{1}t)\text{sinc}({\Omega }_2{r}_2)\text{sinc}({\Omega }_3{r}_3),$$

$$\text{sinc}(\pi {\tau }_1^{-1}(t-i_1{\tau }_1))\text{sinc}(\pi {\tau }_2^{-1}(r_2-i_2{\tau }_2))\text{sinc}(\pi {\tau }_3^{-1}(r_3-i_3{\tau }_3))\ast g(t,r_2,r_3)$$

$$={\tau }_1{\tau }_2{\tau }_3\frac{{\Omega }_1{\Omega }_2{\Omega }_3}{{\pi }^3}\text{sinc}({\Omega }_1(t-i_1{\tau }_1))\text{sinc}({\Omega }_2(r_2-i_1{\tau }_2))\text{sinc}({\Omega }_3(r_3-i_1{\tau }_3)),$$

$${\hat{f}}_{n_1,n_2,n_3}(t,r_2,r_3):={\tau }_1{\tau }_2{\tau }_3\sum _{i_1=1}^{n_1}\sum _{i_2=1}^{n_2}\sum _{i_3=1}^{n_3}{y}_{i_1,i_2,i_3}\phi \left(t-i_1{\tau }_1,r_2-i_2{\tau }_2,r_3-i_3{\tau }_3\right),$$

$$\phi (t,r_2,r_3)=\frac{\sin ({\Omega }_{1}t)}{\pi t}\frac{\sin ({\Omega }_2{r}_2)}{\pi r_2}\frac{\sin ({\Omega }_3{r}_3)}{\pi r_3}=:{\tilde{\phi }}_1(t){\tilde{\phi }}_2(r_2){\tilde{\phi }}_3(r_3)$$

$$H_0:f(t,r_2,r_3)=f_0(t,r_2,r_3)$$

$$\left\{y_{i_1,i_2,i_3}:i_1\in \{1,\dots ,k\},i_2\in \{1,\dots ,n_2\},i_3\in \{1,\dots ,n_3\}\right\},$$

$${\mathcal{F}}_n(s,t,r_2,r_3)$$

$$=\sqrt{{\tau }_1{\tau }_2{\tau }_3}\sum _{l_1=1}^{\lfloor n_{1}s\rfloor }\sum _{l_2=1}^{n_2}\sum _{l_3=1}^{n_3}\left[y_{l_1,l_2,l_3}-f_0(l_1{\tau }_1,l_2{\tau }_2,l_3{\tau }_3)\right]$$

$$\phi \left(t-l_1{\tau }_1,r_2-l_2{\tau }_2,r_3-l_3{\tau }_3\right),$$

$${\mathcal{L}}_n:=\min \left\{n_0\le k\le n_1:\underset{\genfrac{}{}{0ex}{}{r_2\in [0,{\overline{\tau }}_2],}{r_3\in [0,{\overline{\tau }}_3]}}{\max }\left|{\mathcal{F}}_n\left(\frac{k}{n_1},\frac{{\overline{\tau }}_{1}k}{n_1},r_2,r_3\right)\right|>c_L\right\}$$

$${\mathcal{M}}_n:=\min \left\{n_0\le k\le n_1:\underset{0\le t\le {\overline{\tau }}_{1}k/n_1}{\max }\underset{\genfrac{}{}{0ex}{}{r_2\in [0,{\overline{\tau }}_2],}{r_3\in [0,{\overline{\tau }}_3]}}{\max }\left|{\mathcal{F}}_n\left(\frac{k}{n_1},t,r_2,r_3\right)\right|>c_M\right\}$$

$$f(t,r_2,r_3)=f_0(t,r_2,r_3),0\le t\le s_0{\overline{\tau }}_1,0<s_0<1.$$

$$\Delta (f;\mathbf{a},\mathbf{b}):=\sum _{v\subseteq \{1,\dots ,q\}}(-1{)}^{|v|}f\left({\mathbf{a}}_v:{\mathbf{b}}_{-v}\right).$$

$$W:=\prod _{j=1}^q(s_j,t_j],$$

$$X(W):=\Delta (X;\mathbf{s},\mathbf{t})=\sum _{v\subseteq \{1,\dots ,q\}}(-1{)}^{|v|}X\left({\mathbf{s}}_v:{\mathbf{t}}_{-v}\right).$$

$$B(W_1),\dots ,B(W_k)$$

$$\lambda (W_1),\dots ,\lambda (W_k),$$

$$Z_n(v_1,v_2,v_3):=(n_1{n}_2{n}_3{)}^{-1/2}\sum _{i_1=1}^{\lfloor n_1{v}_1\rfloor }\sum _{i_2=1}^{\lfloor n_2{v}_2\rfloor }\sum _{i_3=1}^{\lfloor n_3{v}_3\rfloor }{\epsilon }_{i_1,i_2,i_3}\Rightarrow \sigma B$$

$$(v_1,v_2,v_3)\in [0,1{]}^3,$$

$${\sigma }^2:=\underset{{\min }_{1\le i\le 3}{n}_i\to \infty }{\lim }\mathbb{V}\text{ar}\left((n_1{n}_2{n}_3{)}^{-1/2}\sum _{i_1=1}^{n_1}\sum _{i_2=1}^{n_2}\sum _{i_3=1}^{n_3}{\epsilon }_{i_1,i_2,i_3}\right)=\sum _{\mathbf{k}\in {\mathbb{Z}}^3}\mathbb{E}({\epsilon }_0{\epsilon }_{\mathbf{k}}).$$

$$\mathbb{E}({\xi }_0)=0,\mathbb{E}\left({\xi }_0^2\right)<\infty \text{and}{\sigma }^2>0,$$

$${\mathcal{F}}_n(s,t,r_2,r_3)\Rightarrow \mathcal{F}(s,t,r_2,r_3),$$

$$\mathcal{F}(s,t,r_2,r_3):=\sqrt{{\overline{\tau }}_1{\overline{\tau }}_2{\overline{\tau }}_3}\sigma {\int }_0^s{\int }_0^1{\int }_0^1$$

$$\mathrm{d}B(z_1,z_2,z_3),$$

$$\mathbb{E}(\mathcal{F}(s,t,r_2,r_3))=0$$

$$\mathbb{C}\text{ov}\left(\mathcal{F}(s^{(1)},t^{(1)},r_2^{(1)},r_3^{(1)}),\mathcal{F}(s^{(2)},t^{(2)},r_2^{(2)},r_3^{(2)})\right)$$

$$={\sigma }^2{\overline{\tau }}_1{\overline{\tau }}_2{\overline{\tau }}_3{\int }_0^{\min \{s^{(1)},s^{(2)}\}}{\int }_0^1{\int }_0^1\phi (t^{(1)}-{\overline{\tau }}_1{z}_1,r_2^{(1)}-{\overline{\tau }}_2{z}_2,r_3^{(1)}-{\overline{\tau }}_3{z}_3)$$

Full text

 

**Keywords: Change-point problems; Correlated noise random fields; Image processing; Multivariate Brownian motion; Sampling theorems; Sequential detection. **

000Address correspondence to Annabel Prause or Ansgar Steland, Institute of Statistics, RWTH Aachen University, Wüllnerstraße 3, 52062 Aachen, Germany; E-mail: prause@stochastik.rwth-aachen.de or steland@stochastik.rwth-aachen.de000This is a preprint version of an article to appear in Sequential Analysis.

1 Introduction

Signal processing and signal transmission play an important role in many different areas. Typical problems include the reconstruction of a signal by its discretely sampled values as well as the detection of changes from a given reference signal. For univariate signals $f\colon\mathbb{R}\to\mathbb{R}$, sampled equidistantly using a sampling period $\tau>0$ and disturbed by additive noise, such that one obtains a block of noisy samples $\{y_{\bm{i}}:i\in\{1,\ldots,n\}\}$ where

[TABLE]

a nonparametric joint reconstruction/detection algorithm has been proposed in the paper of Pawlak and Steland (2013). Their approach has several appealing features. Firstly, the algorithm can detect changes while reconstructing the signal at the same time. Secondly, it is a nonparametric approach, i.e. no further information about the exact class to which the observed signal belongs is necessary. Lastly, the procedure works in a sequential way such that changes can be detected on-line, in contrast to off-line detection schemes which can first detect changes in retrospect, i.e. when the whole data set is already available.

A natural question arises whether this approach also works for high-dimensional signals. One answer to this problem is given in Prause and Steland (2015), where the authors treat matrix-valued signals and apply results from Pawlak and Steland (2013) by considering quadratic forms. In the present paper, however, we consider a more general framework by focusing our attention on signals $f\colon\mathbb{R}^{q}\to\mathbb{R}$ for $q>1$. Examples for such signals are multifaceted, including geographic and climatic data as well as image data, that are observed over a fixed time horizon. In order to simplify the notation we fix $q=3$ as this case also covers the most interesting applications. However, our results also hold true for $q=2$ and arbitrary $q>3$ and the corresponding proofs can easily be completed along the same lines. Thus, in the following we are interested in reconstructing three-dimensional signals and, at the same time, in detecting changes from a given reference signal. Here, one component represents the time and the other two the location. The application that we have in mind are video signals, i.e. sequences of image frames over time.

The basis on which we now want to establish our investigations is a finite block of noisy samples $\{y_{\bm{i}}=y_{i_{1},i_{2},i_{3}}:(i_{1},i_{2},i_{3})\in\{1,\ldots,n_{1}\}\times\{1,\ldots,n_{2}\}\times\{1,\ldots,n_{3}\}\}$ that – in accordance with model (1.1) – is obtained from the model

[TABLE]

Here, $f(t,r_{2},r_{3})$ is the unknown signal depending on time ($i_{1}$) and location ($i_{2}$ and $i_{3}$), $\left\{\varepsilon_{\bm{i}}=\varepsilon_{i_{1},i_{2},i_{3}}\right\}$ is a zero mean random field and $\tau_{j}=\tau_{n_{j}},\,j=1,2,3$, are the sampling periods. We assume that they fulfill $\tau_{j}\to 0$ and $n_{j}\tau_{j}\to\overline{\tau}_{j}$, $j=1,2,3$ as $\min_{1\leq i\leq 3}n_{i}\to\infty$. As in Pawlak and Steland (2013) we want to base our approaches on classical reconstruction procedures from the signal sampling theory, leading to sequential partial sum processes as detector statistics, see Section 2. In order to make these detector statistics applicable we need to determine proper control limits; thus, in Section 3 we will show that we can generalize the two main weak convergence results in Pawlak and Steland (2013) to our multidimensional context, i.e. we show weak convergence of the detection process towards Gaussian processes under different assumptions on the dependence structure of the noise processes where either the null hypothesis $f=f_{0}$ or the alternative $f\neq f_{0}$ holds true. In Section 4 we present extensions to weighting functions, which allow to detect the location of the change as well, and discuss how to treat the case of an unknown but time-constant reference signal. Finally, in Section 5 we present some simulation results concerning the rejection rates and the power of the detection algorithm.

2 The detection algorithm

We now want to extend the main results of Pawlak and Steland (2013) to signals $f\colon\mathbb{R}^{q}\to\mathbb{R}$ with $q=3$. As in Pawlak and Steland (2013) we base our estimator of $f(t,r_{2},r_{3})$ on results of the signal sampling theory like the Shannon-Whittaker theorem. This theorem has generalizations to signals with several variables. In three dimensions we have for band-limited functions on $[-\Omega_{1},\Omega_{1}]\times[-\Omega_{2},\Omega_{2}]\times[-\Omega_{3},\Omega_{3}]$ with $0<\Omega_{1},\Omega_{2},\Omega_{3}<\infty$ the representation

[TABLE]

where $\textnormal{sinc}(x)=\frac{\sin(x)}{x}$, cf. Jerri (1977), p. 1571. The most direct idea to construct an estimator of $f(t,r_{2},r_{3})$ would now be to just replace the values of $f$ at $\left(i_{1}\pi/\Omega_{1},i_{2}\pi/\Omega_{2},i_{3}\pi/\Omega_{3}\right)$ by the noisy observations $y_{\bm{i}}=y_{i_{1},i_{2},i_{3}}$. However, Pawlak and Stadtmüller (1996) have shown that this naive estimator is not even consistent in one dimension. Instead, they propose a post-filtering correction of the so-called oversampled version of the Shannon-Whittaker series to filter out high frequencies. This is the approach that we also adapt here. In three dimensions this oversampled version is of the form

[TABLE]

for $0<\tau_{j}\leq\pi/\Omega_{j},j=1,2,3$. If we convolve this version with

[TABLE]

use the fact that

[TABLE]

and replace $f\left(i_{1}\tau_{1},i_{2}\tau_{2},i_{3}\tau_{3}\right)$ by the noisy sample $y_{i_{1},i_{2},i_{3}}$, this finally leads to the following truncated convolution form as an estimator for $f(t,r_{2},r_{3}),$ namely

[TABLE]

cf. Pawlak and Stadtmüller (1996), p. 1427, and Pawlak and Stadtmüller (2007), p. 2527. Here,

[TABLE]

is a three-dimensional product reconstruction kernel with $\widetilde{\varphi}_{j}(0)=\Omega_{j}/\pi$ for $j=1,2,3$.

Given a reference signal $f_{0}(t,r_{2},r_{3})$, our aim now is to decide whether or not we can reject the null hypothesis

[TABLE]

for all $t\in[0,\overline{\tau}_{1}],r_{2}\in[0,\overline{\tau}_{2}],r_{3}\in[0,\overline{\tau}_{3}]$. As we receive our data in a sequential way over time as a sequence of image frames, we want to be able to detect changes as early as possible, i.e. we want to give an alarm as soon as we have enough evidence in our samples

[TABLE]

corresponding to the first $k$ image frames, to reject the null hypothesis. To achieve this aim we consider a sequential partial sum process over time which is defined as

[TABLE]

for $0<s_{0}\leq s\leq 1,t\in[0,\overline{\tau}_{1}],r_{2}\in[0,\overline{\tau}_{2}],r_{3}\in[0,\overline{\tau}_{3}]$. The subscript $n$ denotes the dependence on the sample sizes. Note that $\mathbb{E}_{0}$ denotes as usual the expectation taken under the null hypothesis. With this process we can easily define detectors such as the local detector

[TABLE]

or the global maximum norm detector

[TABLE]

with control limits $c_{L}>0$ and $c_{M}>0$ and $n_{0}\coloneqq\lfloor n_{1}s_{0}\rfloor$, $s_{0}\in(0,1)$. The reason to start monitoring in $n_{0}$ is that we assume that we have a kind of learning sample, i.e. we assume

[TABLE]

This assumption guarantees that no initial change in the signal occurs before the monitoring procedure starts which allows to estimate unknown parameters of the detection process such as the asymptotic variance $\sigma^{2}$, see below.

Now the question arises how to reasonably choose the control parameters $c_{L}$ and $c_{M}$. To answer this question we are interested in the limiting distribution of our detector statistics. These can be derived from the limiting distribution of our stochastic process $\mathcal{F}_{n}(s,t,r_{2},r_{3})$ which is subject of the next section.

3 Process limit distributions

The asymptotic results to be discussed now, under the null hypothesis of no change and a rich class of alternative hypotheses under which the true signal converges to the reference signal, are based on a weak assumption about the asymptotic distribution of the partial sum process of the random field $\left\{\varepsilon_{i_{1},i_{2},i_{3}}\right\}$. Before discussing that assumption and providing the asymptotic results, we introduce some preliminaries and notations used in the sequel.

3.1 Preliminaries

As usual, the cardinality of a set $u\subseteq\{1,\ldots,q\}$ is denoted as $|u|$. Moreover, for $v\subseteq\{1,\ldots,q\}$ we write $u-v$ for the complement of $v$ with respect to $u$. In particular, we just write $-v$ if we take the complement of $v$ with respect to the whole set $\{1,\ldots,q\}$. Sets of the form $\{j,j+1,\ldots,k\}$ for integers $j$ and $k$ with $j\leq k$ are abbreviated as $j:k$, such that $\{1,\ldots,q\}=1:q$.

To pick out the components of a vector $\bm{x}\in\mathbb{R}^{q}$ that correspond to a set $u\subseteq 1:q$, we write $\bm{x}_{u}$, i.e. $\bm{x}_{u}$ stands for a vector with $|u|$ components with selected entries of $\bm{x}$. Now, let $u,v\subseteq 1:q$ and $\bm{x},\bm{z}\in[\bm{a},\bm{b}]$ with $u\cap v=\emptyset$. The symbol $\bm{x}_{u}:\bm{z}_{v}$ then denotes the point $\bm{y}\in\left[\bm{a}_{u\cup v},\bm{b}_{u\cup v}\right]$, where $y_{j}=x_{j}$ for all $j\in u$ and $y_{j}=z_{j}$ for all $j\in v$.

Recall the concept of the $q$-fold alternating sum of a function $f$ over the hyperrectangle $[\bm{a},\bm{b}]$ which is defined as

[TABLE]

Now, let $W$ be a subset of $[0,1]^{q}$. For points in $[0,1]^{q}$ set $\bm{t}=(t_{1},\ldots,t_{q})$ and $\bm{s}=(s_{1},\ldots,s_{q})$ respectively. We call $W$ a block if it is of the form

[TABLE]

where each $(s_{j},t_{j}],j=1,\ldots,q$ is a left-open, right-closed subintervals of $[0,1]$. We now define the increment $X(W)$ of a stochastic process $X=\{X(\bm{t}):\bm{t}\in[0,1]^{q}\}$ around a block $W$ by means of the alternating sum (3.1) as

[TABLE]

We are now able to define the Brownian motion on $[0,1]^{q}$, cf. Deo (1975), p. 709, where $C_{q}$ stands for the space of all real-valued continuous functions on $[0,1]^{q}$.

Definition 3.1**.**

The Brownian motion $B=\{B(\bm{t}):\bm{t}\in[0,1]^{q}\}$ on $[0,1]^{q}$ is characterized by

• (a)

$P(B\in C_{q})=1$,

• (b)

if $W_{1},\ldots,W_{k}$ are pairwise disjoint blocks in $[0,1]^{q}$, then the increments

[TABLE]

are independent normal random variables with means zero and variances

[TABLE]

$\lambda$ being the $q$-dimensional Lebesgue measure on $[0,1]^{q}$.

We can now formulate the main assumption for the asymptotic theory of the detection process as follows, where $\Rightarrow$ – as usual – stands for weak convergence in an appropriately chosen function space, see Appendix A.1.

Assumption 1: Let $\{\varepsilon_{i_{1},i_{2},i_{3}}\}$ be a weakly stationary random field with $\mathbb{E}(\varepsilon_{\bm{i}})=0$ which satisfies a functional central limit theorem, i.e.

[TABLE]

in the Skorohod space $D[0,1]^{3}$ as $\min_{1\leq i\leq 3}n_{i}\to\infty$ for some constant $\sigma^{2}\in(0,\infty)$.

Here, the constant $\sigma^{2}$ equals the long-run variance of the random field $\{\varepsilon_{i_{1},i_{2},i_{3}}\}$, i.e.

[TABLE]

There exist several results in the literature about the weak invariance principle (3.2) under specific conditions on the random field. In particular, in the i.i.d. case we get the functional central limit theorem under the sole assumptions that

[TABLE]

see Corollary 1 in Wichura (1969). More generally, a functional central limit theorem for strictly stationary and $\varphi$-mixing random fields can be found in Deo (1975), cf. Theorem 1. Further results on weak invariance principles for random fields include weakly stationary associated as well as weakly stationary and $\alpha$-mixing random fields, cf. Bulinski and Kaene (1996), p. 2906, and Berkes and Morrow (1981), Theorem 1, respectively. The latter obtain a strong approximation of the partial sum field by a Brownian motion from which one can deduce a weak invariance principle quite directly. Other results on functional central limit theorems for random fields include the ones of Wang and Woodroofe (2013), cf. Theorem 1.1, and Machkouri et al. (2013), cf. Theorem 2. These authors consider random fields of the form $X_{\bm{i}}=g(\varepsilon_{\bm{i}-\bm{s}};\bm{s}\in\mathbb{Z}^{q})$ where $g$ is a measurable function and the $\{\varepsilon_{\bm{j}};\bm{j}\in\mathbb{Z}^{q}\}$ are i.i.d. random variables. Machkouri et al. (2013) introduce the notion of a $p$-stable random field and then obtain a weak invariance principle for the so-called smoothed partial sum process.

3.2 Asymptotics under the Null Hypothesis

With the help of Assumption 1 we are now in a position to formulate the following theorem stating the asymptotic behaviour of the process.

Theorem 3.1**.**

Suppose the noise process $\{\varepsilon_{\bm{i}}=\varepsilon_{i_{1},i_{2},i_{3}}\}$ meets Assumption 1. Assume that the sampling periods fulfill $n_{j}\tau_{j}\to\overline{\tau}_{j}$ for $j=1,2,3,$ as $\min_{1\leq i\leq 3}n_{i}\to\infty$. Then, under the null hypothesis $H_{0}$, we have

[TABLE]

as $\min_{1\leq i\leq 3}n_{i}\to\infty$ for $0<s_{0}\leq s\leq 1$, $t\in[0,\overline{\tau}_{1}]$, $r_{2}\in[0,\overline{\tau}_{2}]$ and $r_{3}\in[0,\overline{\tau}_{3}]$.

The limit stochastic process $\mathcal{F}(s,t,r_{2},r_{3})$ is of the form

[TABLE]

where $B(z_{1},z_{2},z_{3})$ is the standard Brownian motion on $[0,1]^{3}$.

The weak convergence takes place in a higher dimensional Skorohod space and the last integral is interpreted as multivariate Riemann-Stieltjes integral, see Appendix A.1 for more details.

The next lemma is a characterization of the correlation structure of the limit process $\mathcal{F}(s,t,r_{2},r_{3})$.

Lemma 3.1**.**

• (a)

The process $\mathcal{F}(s,t,r_{2},r_{3})$ is a nonstationary multivariable Gaussian process with

[TABLE]

and covariance function

[TABLE]

for $0<s_{0}\leq s^{(1)},s^{(2)}\leq 1$, $0\leq t^{(1)},t^{(2)}\leq\overline{\tau}_{1}$, $0\leq r^{(1)}_{2},r^{(2)}_{2}\leq\overline{\tau}_{2}$, and $0\leq r^{(1)}_{3},r^{(2)}_{3}\leq\overline{\tau}_{3}$.

• (b)

The process $\mathcal{F}(s,t,r_{2},r_{3})$ has continuous sample paths.

Now that we have the limit distribution of $\mathcal{F}_{n}(s,t,r_{2},r_{3})$ under the null hypothesis at our disposal, we can easily derive central limit theorems for the local and global maximum norm detector defined in (2.2) and (2.3).

Lemma 3.2**.**

Assume that condition (2.4) holds true. Then, under the conditions of Theorem 3.1 the detectors satisfy the following central limit theorems:

[TABLE]

as $\min_{1\leq i\leq 3}n_{i}\to\infty$.

3.3 Asymptotics under the Alternative

We now investigate the behaviour of our statistic $\mathcal{F}_{n}(s,t,r_{2},r_{3})$ under a general class of alternatives $H_{1}:f_{n_{1},n_{2},n_{3}}(t,r_{2},r_{3})\neq f_{0}(t,r_{2},r_{3})$, i.e. in the situation when the observed signal and the reference signal differ. We assume that our observed data $\{y_{\bm{i}}=y_{i_{1},i_{2},i_{3}}:(i_{1},i_{2},i_{3})\in\{1,\ldots,n_{1}\}\times\{1,\ldots,n_{2}\}\times\{1,\ldots,n_{3}\}\}$ obey the following model:

[TABLE]

with the true signal $f_{n_{1},n_{2},n_{3}}(t,r_{2},r_{3})$ depending on the sample size $(n_{1},n_{2},n_{3})$ and $f_{n_{1},n_{2},n_{3}}(t,r_{2},r_{3})\to f_{0}(t,r_{2},r_{3})$ as $\min_{1\leq i\leq 3}n_{i}\to\infty$. The process $\{\varepsilon_{\bm{i}}=\varepsilon_{i_{1},i_{2},i_{3}}\}$ is again the zero mean noise random field fulfilling Assumption 1.

It turns out that the process $\mathcal{F}_{n}(s,t,r_{2},r_{3})$ converges to a well-defined and non-degenerate limit process under general conditions on the variation of the difference $f_{n_{1},n_{2},n_{3}}(t,r_{2},r_{3})-f_{0}(t,r_{2},r_{3})$, similar as in Pawlak and Steland (2013). However, whereas in dimension $q=1$ the Vitali variation suffices, in higher dimensions one has to consider the variation in the sense of Hardy and Krause.

Let $a,b\in\mathbb{R}$ with $a\leq b$. A ladder on $[a,b]$ is a set $\mathcal{Y}$ containing $a$ and finitely many, possibly zero, values from $(a,b)$, see Owen (2005), p. 2. The successor of an element $y\in\mathcal{Y}$ is denoted by $y^{+}$. For $(y,\infty)\cap\mathcal{Y}=\emptyset$ we set $y^{+}=b$ and otherwise $y^{+}$ is the smallest element of $(y,\infty)\cap\mathcal{Y}$. In particular, if we consider a classical partition of $[a,b]$, we have $\mathcal{Y}=\{y_{0},y_{1},\ldots,y_{m}\}$ with $a=y_{0}<y_{1}<\ldots<y_{m}$ such that $y_{k+1}$ is the successor of $y_{k}$ for all $k<m$ and it is $b$ for $k=m$.

If we now define $\mathbb{Y}$ as the set of all ladders on $[a,b]$, we can write the total variation of a univariate function $f$ defined on $[a,b]$ as

[TABLE]

In order to generalize the one-dimensional variation to the multidimensional case we need the concept of multidimensional ladders. We now consider a hyperrectangle $[\bm{a},\bm{b}]$ with $\bm{a},\bm{b}\in\mathbb{R}^{q}$ and $\bm{a}\leq\bm{b}$. We define a ladder $\mathcal{Y}$ on $[\bm{a},\bm{b}]$ as $\mathcal{Y}\coloneqq\prod_{j=1}^{q}\mathcal{Y}_{j}$, where $\mathcal{Y}_{j}$ is a ladder on $\left[a_{j},b_{j}\right]$ for $j=1,\ldots,q$. Similarly, we say that $\bm{y}^{+}$ is the successor of $\bm{y}$ if $y^{+}_{j}$ is the successor of $y_{j}$ for each $j=1,\ldots,q$. Finally, we set $\mathbb{Y}\coloneqq\prod_{j=1}^{q}\mathbb{Y}_{j}$ for the set of all ladders on $[\bm{a},\bm{b}]$, where $\mathbb{Y}_{j}$ denotes the set of all ladders on $\left[a_{j},b_{j}\right]$ for $j=1,\ldots,q$. With these alternating sums we are now able to define the variation of $f$ over $\mathcal{Y}$ as

[TABLE]

This leads to the following definition, see Owen (2005), Definition 1 and 2.

Definition 3.2**.**

The variation of $f$ on the hyperrectangle $[\bm{a},\bm{b}]$, in the sense of Vitali, is

[TABLE]

The variation of $f$ on the hyperrectangle $[\bm{a},\bm{b}]$, in the sense of Hardy and Krause, is

[TABLE]

Likewise as in the one-dimensional case, we say that the function $f$ is of bounded variation in the sense of Vitali (and we write $f\in BV$ or $f\in BV[\bm{a},\bm{b}]$) if $V(f)<\infty$. Note that for $q=1$ the variation in the sense of Vitali corresponds with the common definition of the variation of a univariate function. We say that the function $f$ is of bounded variation in the sense of Hardy and Krause (and we write $f\in BVHK$ or $f\in BVHK[\bm{a},\bm{b}]$) if $V_{HK}(f)<\infty$. Note that the summand for $u=1:q$ in the above sum equals $V(f)$, such that $V(f)<\infty$ if $V_{HK}(f)<\infty$. Moreover, it was shown by Young (1917), that $V_{HK}(f)<\infty$ if and only if $V_{[\bm{a},\bm{b}]}(f)<\infty$ and $V_{\left[\bm{a}_{u},\bm{b}_{u}\right]}f\left(\bm{x}_{u}:\bm{z}_{-u}\right)<\infty$ for all $0<|u|<q$ and all $\bm{z}_{-u}\in\left[\bm{a}_{-u},\bm{b}_{-u}\right]$, which is the original definition of bounded variation of Hardy, see Hardy (1906). This means, that in the above definition $\bm{b}_{-u}$ could be replaced by an arbitrary fixed point of the hyperrectangle $\left[\bm{a}_{-u},\bm{b}_{-u}\right]$ for $\emptyset\neq u\subseteq 1:q$.

For more details on these two notions of multidimensional variation we refer the reader to Owen (2005).

3.3.1 Deterministic Disturbances

We first consider the case where $f_{n_{1},n_{2},n_{3}}(t,r_{2},r_{3})$ is deterministic and start with the following illustrative example, where $n=n_{1}=n_{2}=n_{3}$. Consider

[TABLE]

for some $\theta\in(0,1)$, $\gamma\geq 0$, $\beta>0$, and a non-zero ‘location’ function $\widetilde{\delta}(r_{2},r_{3})$ defined on $[0,\overline{\tau}_{2}]\times[0,\overline{\tau}_{3}]$. This is a local alternative with a change point at time $t=\theta\overline{\tau}_{1}$. This means, that up to time $\theta\overline{\tau}_{1}$ the observed data obey $f_{0}$ and after this point in time they get disturbed by $\left(t-\theta\overline{\tau}_{1}\right)^{\gamma}\widetilde{\delta}(r_{2},r_{3})/n^{\beta}$. This disturbance depends on the function $\widetilde{\delta}(r_{2},r_{3})$ which assigns different weights at the locations $(r_{2},r_{3})$ changing with the time. In the following we require a more general model for local alternatives, namely we consider

[TABLE]

for some $\beta_{1},\beta_{2},\beta_{3}>0$ and a deterministic function $\delta(t,r_{2},r_{3})$. We assume that $\delta$ meets the following assumption.

Assumption 2: Let $\delta(t,r_{2},r_{3})$ be a nonzero function defined on $[0,\overline{\tau}_{1}]\times[0,\overline{\tau}_{2}]\times[0,\overline{\tau}_{3}]$ which is

• (a)

continuous

• (b)

of bounded variation in the sense of Hardy and Krause.

Now we are able to describe the asymptotic behaviour of the statistic $\mathcal{F}_{n}(s,t,r_{2},r_{3})$ under local alternatives.

Theorem 3.2**.**

Assume the sampling model in (3.4) with the local alternative given in (3.6) where $\beta=3/2$. Let Assumption 1 hold and suppose that either Assumption 2 (a) and $n\tau_{j}\to\overline{\tau}_{j}$, $n\to\infty$, $j=1,2,3,$ or Assumption 2 (b) and $n\tau_{j}=\overline{\tau}_{j}$, $j=1,2,3,$ is satisfied. Then, we have

[TABLE]

as $\min_{1\leq i\leq 3}n_{i}\to\infty$ for $0<s_{0}\leq s\leq 1$, $t\in[0,\overline{\tau}_{1}]$, $r_{2}\in[0,\overline{\tau}_{2}]$ and $r_{3}\in[0,\overline{\tau}_{3}]$.

The limit stochastic process $\mathcal{F}^{\delta}(s,t,r_{2},r_{3})$ is given by

[TABLE]

Similarly as above we obtain central limit theorems for our detectors, defined in (2.2) and (2.3), under the alternative.

Lemma 3.3**.**

Let the condition in (2.4) hold. Then, under the assumptions of Theorem 3.2 we obtain the asymptotic distribution of the local and global maximum norm detector by replacing $\mathcal{F}(s,t,r_{2},r_{3})$ by $\mathcal{F}^{\delta}(s,t,r_{2},r_{3})$ in Corollary 3.2.

3.3.2 Random Disturbances

We now consider random disturbances, i.e. we require that our data obey the model

[TABLE]

where now $f_{n}(t,r_{2},r_{3};\omega)\to f_{0}(t,r_{2},r_{3})$ a.s. as $n\to\infty$. To be more precise we require that

[TABLE]

for $\beta>0$, and $\Delta(t,r_{2},r_{3};\omega)$ being a random function that is independent of the random field $\{\varepsilon_{i_{1},i_{2},i_{3}}\}$. Moreover, we assume that $\Delta(t,r_{2},r_{3};\omega)\neq 0$ a.s. We require that $\Delta(t,r_{2},r_{3};\omega)$ meets the following assumption.

Assumption 3: Let $\Delta(t,r_{2},r_{3};\omega)$ be an a.s. nonzero random function defined on $[0,\overline{\tau}_{1}]\times[0,\overline{\tau}_{2}]\times[0,\overline{\tau}_{3}]\times\Omega$ that is independent of the random field $\{\varepsilon_{i_{1},i_{2},i_{3}}\}$ and whose sample paths are

• (a)

continuous a.s.

• (b)

of bounded variation in the sense of Hardy and Krause a.s.

Then we can describe the asymptotic behaviour of the statistic $\mathcal{F}_{n}(s,t,r_{2},r_{3})$ under random local alternatives.

Theorem 3.3**.**

Assume the sampling model in (3.7) with the local alternative given in (3.8) where $\beta=3/2$. Let Assumption 1 hold and suppose that either Assumption 3 (a) and $n\tau_{j}\to\overline{\tau}_{j}$, $n\to\infty$, $j=1,2,3,$ or Assumption 3 (b) and $n\tau_{j}=\overline{\tau}_{j}$, $j=1,2,3,$ is satisfied. Then, we have

[TABLE]

as $\min_{1\leq i\leq 3}n_{i}\to\infty$ for $0<s_{0}\leq s\leq 1$, $t\in[0,\overline{\tau}_{1}]$, $r_{2}\in[0,\overline{\tau}_{2}]$ and $r_{3}\in[0,\overline{\tau}_{3}]$.

The limit stochastic process $\mathcal{F}^{\Delta}(s,t,r_{2},r_{3})$ is given by

[TABLE]

4 Extensions: Weighting functions and unknown reference signals

In this section we want to demonstrate that we can easily extend some of the results of the previous section into several directions with only slight modifications. We give two examples that shall show the great flexibility and applicability of our results. We point out, among other things, that the detector statistic in (2) can serve with small changes not only as a detector for changes in time, but also as a detector for the concrete location where a change takes place. Moreover, we can extend the result in Theorem 3.1 to allow for unknown reference signals, by appropriate centering of the observations.

4.1 Additional Weighting Functions

We begin with a generalization of the detector statistic in (2) in order to be able to detect the position of a change. This can be achieved by adding a suitable weighting function $w$ for the different pixels of the image and leads to the sequential monitoring process

[TABLE]

for $0<s_{0}\leq s\leq 1,t\in[0,\overline{\tau}_{1}],r_{2}\in[0,\overline{\tau}_{2}],r_{3}\in[0,\overline{\tau}_{3}]$. Before we get more specific about possible forms of $w$, we first want to reformulate Theorem 3.1 for the new detection process $\mathcal{F}^{w}_{n}$.

Theorem 4.1**.**

Let $w$ be a continuous function with

[TABLE]

Suppose the noise process $\{\varepsilon_{\bm{i}}=\varepsilon_{i_{1},i_{2},i_{3}}\}$ meets Assumption 1. We assume that the sampling periods fulfill $n_{j}\tau_{j}\to\overline{\tau}_{j}$, $j=1,2,3,$ as $\min_{1\leq i\leq 3}n_{i}\to\infty$. Then, under the null hypothesis $H_{0}$, we have

[TABLE]

as $\min_{1\leq i\leq 3}n_{i}\to\infty$ for $0<s_{0}\leq s\leq 1$, $t\in[0,\overline{\tau}_{1}]$, $r_{2}\in[0,\overline{\tau}_{2}]$, and $r_{3}\in[0,\overline{\tau}_{3}]$.

The limit stochastic process $\mathcal{F}^{w}(s,t,r_{2},r_{3})$ is of the form

[TABLE]

where $B(z_{1},z_{2},z_{3})$ is the standard Brownian motion on $[0,1]^{3}$.

Now the question arises which forms would be suitable for $w$. To answer this question it is useful to know the approximate form of the change that one wants to detect. If the aim is, for example, to detect a simple rectangle with length $c_{2}>0$ and width $c_{3}>0$ that occurs at a certain point in time, one could define $w$ for $\delta>0$ as

[TABLE]

Then, one could define suitable detectors as before, e.g. the local detector as

[TABLE]

and the global maximum norm detector as

[TABLE]

for control limits $c^{w}_{L}>0$ and $c^{w}_{M}>0$. Again, the application of the continuous mapping theorem leads to central limit theorems for these detectors. If, for example, the local detector $\mathcal{L}_{n}$ exceeds its control limit $c^{w}_{L}$ for $(k^{\star},l^{\star},m^{\star})$, we know that the rectangle occured on the $\overline{\tau}_{1}k^{\star}/n_{1}$th image frame at position $\left(\overline{\tau}_{2}l^{\star}/n_{2},\overline{\tau}_{3}m^{\star}/n_{3}\right)$.

Remark 4.1**.**

It is important that one chooses the weight function $w$ not only as a characteristic function, as then the continuity assumption of Theorem 4.1 on $w$ is not fulfilled. Moreover, characteristic functions with a domain that is not parallel to the axes have infinite variation which encourages the ‘smoothing’ of $w$ as well, cf. Owen (2005), p. 14. If one also wants to allow for characteristic functions without smoothing one has to consider the integrals as Itô integrals instead of considering them as Riemann Stieltjes integrals as it is done in this work.

The detection of more complex forms for changes than rectangles is possible as well by choosing corresponding domains for the smoothed characteristic function in (4.1).

4.2 The Case of an Unknown Reference Signal

If the reference signal is unknown to us, the above procedures are not applicable. To overcome this drawback, we shall assume in the sequel that the reference signal is time-constant, i.e.

[TABLE]

holds true for all $t\in[0,\overline{\tau}_{1}]$. In this case one may center the spatial-temporal observations at appropriately defined averages of previous observations. Here one can either use the learning sample

[TABLE]

or include all observations available at the current time instant. For $(l_{2},l_{3})\in\left\{1,\ldots,n_{2}\right\}\times\left\{1,\ldots,n_{3}\right\}$ let us define

[TABLE]

and consider

[TABLE]

We then get the following theorem.

Theorem 4.2**.**

Under the conditions of Theorem 3.1 and if the reference signal fulfills (4.2), we have under the null hypothesis

[TABLE]

as $\min_{1\leq i\leq 3}n_{i}\to\infty$ for $0<s_{0}\leq s\leq 1$, $t\in[0,\overline{\tau}_{1}]$, $r_{2}\in[0,\overline{\tau}_{2}]$, and $r_{3}\in[0,\overline{\tau}_{3}]$.

The limit stochastic process $\mathcal{F}^{c}(s,t,r_{2},r_{3})$ is of the form

[TABLE]

where

[TABLE]

for $z_{1}\in[s_{0},1]$, $(z_{2},z_{3})\in[0,1]^{2}$.

5 Simulation of the detection process

In this section we want to investigate the performance of the global maximum norm detector defined in (2.3). Before we do so, we have to explain how we calculate an appropriate control limit $c_{M}$ such that we can guarantee that the asymptotic false alarm probability is smaller than a given $\alpha$. We first note that for the global maximum norm detector $\mathcal{M}_{n}$ a type one error occurs if $\mathcal{M}_{n}/n_{1}<1$. Next, we adapt Theorem 2 of Pawlak and Steland (2013) to our situation in order to obtain a more explicit formula for the type one error. For that we assume condition (2.4), namely that

[TABLE]

i.e. we assume that there is no initial change in the signal.

Theorem 5.1**.**

Let the condition in (2.4) hold. Under the assumptions of Theorem 3.1 we have

[TABLE]

where $\mathcal{F}(s,t,r_{2},r_{3})$ is a zero mean Gaussian process with the covariance function defined in ((a)).

Because of this theorem we can ensure that the asymptotic false alarm probability is not greater than $\alpha\in(0,1)$, if we choose the control limit $c_{M}=c_{M}(\alpha)$ as the smallest $c_{M}$ such that

[TABLE]

As it is, however, not easy to obtain a concrete formula for the distribution of

[TABLE]

we propose the following Monte Carlo algorithm to simulate $X$ and the control limit $c_{M}$. The algorithm is an adaption of the proposed algorithm in Pawlak and Steland (2013), p. 8, to the multivariable process $\mathcal{F}(s,t,r_{2},r_{3})$.

Step 1: Generate trajectories of the Gaussian process $\mathcal{F}(s,t,r_{2},r_{3})$ on a grid $\{(s_{i},t_{j},(r_{2})_{k},(r_{3})_{l}):i,j,k,l=1,\ldots,N\}$ where $0\leq s_{1}<\ldots<s_{N}\leq 1$, $0\leq t_{1}<\ldots<t_{N}\leq\overline{\tau}_{1}$, $0\leq(r_{2})_{1}<\ldots<(r_{2})_{N}\leq\overline{\tau}_{2}$ and $0\leq(r_{3})_{1}<\ldots<(r_{3})_{N}\leq\overline{\tau}_{3}$ for some $N\in\mathbb{N}$.

Step 2: Return $X$ by calculating the maximum of the values $|\mathcal{F}(s_{i},t_{j},(r_{2})_{k},(r_{3})_{l})|$ for all $(i,j,k,l)$ such that the constraints $s_{0}\leq s_{i}\leq 1$ and $0\leq t_{j}\leq s_{i}\overline{\tau}_{1}$ are satisfied.

Step 3: Using a large numer of repetitions of Step 1 and Step 2 produce realizations of $X$ to determine the empirical $(1-\alpha)-$quantile as an approximation for $c_{M}(\alpha)$.

We now begin our investigations with an illustrative example of the detection scheme. For that we assume that our reference signal is given by

[TABLE]

on $[0,2]^{3}$. Moreover, we assume that at the point in time $t=1$ a jump occurs over the whole image sequence which leads to an alternative signal of the form

[TABLE]

with $(t,r_{2},r_{3})\in[0,2]^{3}$. Thus, we obtain our noisy sample by the model

[TABLE]

for $j=0$ (null hypothesis) resp. $j=1$ (alternative), where $\{\varepsilon_{i_{1},i_{2},i_{3}}\}$ is an i.i.d. $\mathcal{N}(0,1)$-distributed random field. In the following we take $\tau_{1}=0.04$ and $\tau_{2}=\tau_{3}=0.05$ corresponding to $n_{1}=\overline{\tau}_{1}/\tau_{1}=50$ observations in the time domain, and to $n_{2}=\overline{\tau}_{2}/\tau_{2}=40$ and $n_{3}=\overline{\tau}_{3}/\tau_{3}=40$ observations in the spatial domain. Moreover, we take the bandwidths as $\Omega_{1}=\Omega_{2}=\Omega_{3}=10$ and choose $s_{0}=0.05$ leading to $n_{0}=\left\lfloor s_{0}n_{1}\right\rfloor=2$.

We now consider the global maximum norm detector

[TABLE]

for $k=n_{0},\ldots,n_{1}$. If we take $\alpha=0.05$ and apply the Monte Carlo algorithm from above, we obtain as value for the control limit $c_{M}=23.3147$ which is the horizontal line in Figure 1. Furthermore, the solid line corresponds to the detection process under the alternative whereas the dashed line corresponds to the detection process under the null hypothesis. We can see that the partial sum process stays below the control limit for the whole observation period $[0,2]$ if there is no change in the signal. If we have, however, a change-point at $t=1$ the detection process directly reacts and crosses the control limit a short while later, namely for $k=30$ corresponding to the point in time $t_{0}=30\cdot 0.04=1.2$.

In the following simulation study we investigate the accuracy of the global maximum norm detector. Moreover, we evaluate the influence of different sampling periods and different correlation structures of the noise process. We also want to find out the influence of the asymptotic variance and its estimator developed in Prause and Steland (2016) on the proper selection of the control limit $c_{M}$.

5.1 Influence of the Sampling Periods

We begin by analyzing the influence of different sampling periods in the spatial and time domain with respect to the rejection rates. For that, we calculate the corresponding control limit $c_{M}$ with the help of the Monte Carlo algorithm described above. Thus, we evaluate the process $\mathcal{F}(s,t,r_{2},r_{3})$ on the grid $\{(s_{i},t_{j},(r_{2})_{k},(r_{3})_{l}):i,j,k,l=1,\ldots,N\}$ with $N=15$. After calculating the required maxima of $|\mathcal{F}(s_{i},t_{j},(r_{2})_{k},(r_{3})_{l})|$ we use the $95\%$-quantile of 10000 simulation replicates to estimate $c_{M}$.

In the following we adapt the setting of the illustrative example. As the noise process is modelled by an i.i.d. $\mathcal{N}(0,1)$-distributed random field, we obtain for the asymptotic variance $\sigma^{2}=1$.

Table 1 shows the simulated type one errors for various sampling periods $\tau_{1},\tau_{2}$, and $\tau_{3}$ for 1000 repetitions. We can see that the simulated rejection rates lie between 0.055 and 0.084 and thus that there is only a small influence of the sampling periods on the accuracy of the detector.

If we use 10000 instead of 1000 repetitions, we get even more accurate results. On account of very high computational costs, however, we only examined the rejection rates for 10000 repetitions in four cases.

5.2 Influence of Noise Correlations

In this subsection we investigate how the rejection rates behave when using model (M4) as in Prause and Steland (2016) for the noise process instead of taking i.i.d. errors, i.e. we put for $\rho\in(-1,1)$

[TABLE]

where $u_{t,i,j}\stackrel{{\scriptstyle i.i.d.}}{{\sim}}\mathcal{N}(0,1)$ for all $t,i,j\in\mathbb{Z}$ and the $v_{t,i,j}$ follow for each fixed $t$ their model (M1) which means that

[TABLE]

for i.i.d. innovations $\eta_{i,j}$ with $\eta_{i,j}\sim\mathcal{N}(0,1)$ for all $i$ and $j$ and real weights $a_{k},\ k=1,\ldots,9$, where $a_{5}=1$ and $a_{k}=a$ for $k\neq 5$. Moreover, we suppose that the $v_{t,i,j}$ are uncorrelated for different values of $t$ and that $u_{t_{1},i_{1},j_{1}}$ and $v_{t_{2},i_{2},j_{2}}$ are uncorrelated for all $t_{1},t_{2},i_{1},i_{2},j_{1},j_{2}\in\mathbb{Z}$.

We now fix $\tau_{1}=0.04,\tau_{2}=\tau_{3}=0.05$, and $u_{0}=0.25$ leading to $n_{0}=12$ as size for the learning sample in the time domain. The rest of the setting stays the same as in the illustrative example. We allow the autoregressive parameter $\rho$ to vary over the set $\{0.1,0.2,0.3,0.4,0.5\}$ while the moving average parameter $a$ lies in the set $\{\-0.1,0.01,0.1,0.3,0.5\}$. By Theorem 3.1 we now obtain the proper control limit via the formula $c_{M,\sigma}=\sigma c_{M}$ where $\sigma$ is the asymptotic standard deviation of the noise process and $c_{M}$ the control limit for i.i.d. error terms. If the dependence structure of the noise process is unknown, one has to replace $\sigma$ by a proper estimator, see Prause and Steland (2016), leading to a control limit of the form $c_{M,\widehat{\sigma}}=\widehat{\sigma}c_{M}$.

Table 3 shows the rejection rates for 1000 repetitions for control limits calculated with the true $\sigma$. We see that the rejection rates of the detector are quite accurate over the whole set of considered parameters. Smaller values of $\rho$ and $a$, reflecting a weak dependence structure of the noise process, lead to higher rejection rates, while larger values of these parameters, reflecting a strong dependence of the error terms, lead to lower rejection rates. Moreover, we can see that in most cases the rejection rates decrease for fixed $a$ and growing $\rho$ as well as for fixed $\rho$ and growing $a$, where this decrease is greater for smaller than for larger values of $a$ and $\rho$ respectively.

When the asymptotic variance is unknown we can use the estimators proposed in Prause and Steland (2016) in order to calculate the proper control limits. In this paper the accuracy of the estimators is shown in an extensive simulation study for various models of the underlying random field.

5.3 Power Study

We now want to analyse the power of our detection scheme when allowing for alternatives of the form (3.6). For that we take the local departure models as

[TABLE]

with

[TABLE]

for $(t,r_{2},r_{3})\in[0,2]^{3}$ and $(g_{1},g_{2},g_{3})=(68,125,559)$. Moreover, we put $\beta_{2}=\beta_{3}=0.5$ and choose $\beta_{1}\in\{0.3,0.5,1\}$.

Note that the alternative signal $f_{n_{1},n_{2},n_{3}}(t,r_{2},r_{3})$ displays an amplitude, a frequency as well as a phase distortion in the time domain compared to the reference signal $f_{0}(t,r_{2},r_{3})=\sin(6t)\sin(4r_{2})\sin(4r_{3})$.

By the asymptotic theory of Theorem 3.2 we expect that changes with $\beta_{1}>1/2$ cannot be detected for large sample sizes in the time domain, whereas changes with $\beta_{1}\leq 1/2$ can easily be detected. In particular, the case $\beta_{1}=1/2$ corresponds to the case where the process $\mathcal{F}_{n}(s,t,r_{2},r_{3})$ converges to the process $\mathcal{F}^{\delta}(s,t,r_{2},r_{3})$ defined in Theorem 3.2.

The values of $g_{i},i=1,2,3,$ were chosen in such a way that the initial power for $n_{1}=20$ observations in the time domain is reasonably high and the same for all three different values of $\beta_{1}$. Here, $g_{1}=68$ corresponds to $\beta_{1}=0.3$, $g_{2}=125$ to $\beta_{1}=0.5$, and $g_{3}=559$ to $\beta_{1}=1$.

The rest of the simulation setting stays the same as in the illustrative example; in particular, we suppose that the error terms are i.i.d. and that this is known. The resulting power curves are shown in Figure 2.

Due to very high computational costs the sample size $n_{1}$ only varies between $20$ and $80$. We see that the curves indeed confirm what is predicted by the theory: For $\beta_{1}<1/2$ the power increases towards one and for $\beta_{1}>1/2$ it decreases towards the type one error rate of $\alpha=5\%$. For $\beta_{1}=1/2$ one could assume that the power increases as well when looking at Figure 2; Table 4 shows, however, that for larger sample sizes in the time domain it levels off between 0.7 and 0.8 which thus confirms the theory as well.

Appendix A Proofs

A.1 Preliminaries

Before proving the main results of this paper we briefly want to review the Skorohod space for functions defined on $[0,1]^{q}$ as well as the multivariate Riemann-Stieltjes integral.

The functional $\mathcal{F}_{n}(s,t,r_{2},r_{3})$ and its variations that were introduced in the previous sections can be viewed as an element of the Skorohod space $D_{q}=D[0,1]^{q}$ for functions defined on $[0,1]^{q}$ with $q\in\mathbb{N}$. This function space is the generalisation of the well-known Skorohod space $D=D[0,1]$ for functions with a single time parameter to functions with several time parameters and thus allows for certain discontinuities.

Informally speaking, the space $D[0,1]^{q}$ contains all real-valued functions that are ‘continuous from above, with limits from below’. To be more precise, let $\bm{t}\in[0,1]^{q}$ and define $R_{p}$, for all $1\leq p\leq q$, as one of the relations ‘$<$’ and ‘$\geq$’. Denote by $Q_{R_{1},\ldots,R_{q}}(\bm{t})$ the quadrant

[TABLE]

We say that $f$ is an element of $D_{q}$ if the following two conditions hold:

• (a)

The limit

[TABLE]

exists for each of the $2^{q}$ quadrants $Q=Q_{R_{1},\ldots,R_{q}}(\bm{t})$.

• (b)

We have

[TABLE]

These conditions are analogously to the ones for $D[0,1]$. For more details, including the notion of weak convergence in $D[0,1]^{q}$, we refer to Straf (1972) and Bickel and Wichura (1971) respectively.

The generalization of the univariate Riemann-Stieltjes integral to the multivariate Riemann-Stieltjes integral in $q$ dimensions is now done by means of the $q$-fold alternating sum of a function in (3.1), cf. also Sard (1963).

Let $\mathcal{Y}=\prod_{j=1}^{q}\mathcal{Y}_{j}$ be a ladder on the hyperrectangle $[\bm{a},\bm{b}]\subset\mathbb{R}^{q}$, where $\mathcal{Y}_{j}$ is a ladder on $\left[a_{j},b_{j}\right]$ for all $j=1,\ldots,q$. Let $\widetilde{\bm{y}}\in\left[\bm{y},\bm{y}^{+}\right]$ for each $\bm{y}\in\mathcal{Y}$, where $\bm{y}^{+}$ denotes the component-by-component successor of $\bm{y}$, see above. Define a norm $\|\mathcal{Y}\|$ on $[\bm{a},\bm{b}]$ as

[TABLE]

Suppose that $f$ and $h$ are real-valued functions defined on the hyperrectangle $[\bm{a},\bm{b}]$. Now, for each ladder $\mathcal{Y}\in\mathbb{Y}=\prod_{j=1}^{q}\mathbb{Y}_{j}$ consider the so-called Riemann-Stieltjes sum

[TABLE]

Analogously to the one-dimensional case we now define the Riemann-Stieltjes integral of $h$ with respect to $f$ as

[TABLE]

if the latter exists. The integral on the left is understood as a multivariate integral, namely as

[TABLE]

Similar to the one-dimensional case this integral and all ‘lower dimensional’ integrals exist if $h$ is a continuous function on the hyperrectangle $[\bm{a},\bm{b}]$, $\bm{a},\bm{b}\in\mathbb{R}^{q}$ and if $f$ is a function of bounded variation in the sense of Hardy and Krause on $[\bm{a},\bm{b}]$.

Moreover, there exists a generalization of the integration by parts formula for multivariate Riemann-Stieltjes integrals, see Young (1917), p. 287. This allows us to define the multivariate Riemann-Stieltjes integral even with respect to functions $h$ that are not of bounded variation in the sense of Vitali or in the sense of Hardy and Krause. In this case we take the integration by parts formula as a definition for the integral, i.e. we put

[TABLE]

whenever the right-hand side exists. The square bracket notation is used for the evaluation, respectively, the increment of a multivariate antiderivative $f$ over a hyperrectangle $[\bm{a},\bm{b}]\subset\mathbb{R}^{q}$. Thus, we have

[TABLE]

If the evaluation of $f$ only takes place over a hyperrectangle $\left[\bm{a}_{-v},\bm{b}_{-v}\right],\ v\subset 1:q$, we write $\left[f(\bm{x})\right]_{\bm{a}_{-v}}^{\bm{b}_{-v}}$, and this is defined as

[TABLE]

Thus, we can define the Riemann-Stieltjes integral with respect to a (multivariable) Brownian motion, where the sample paths are of unbounded variation almost surely.

A.2 Proofs

This section is devoted to rigorous proofs of the presented theorems. Note that parts of the material are taken from Prause (2015). We start with proving Theorem 3.1 and introduce the following abbreviating notation. Let

[TABLE]

for $v_{1},v_{2},v_{3}\in[0,1]$, $t\in[0,\overline{\tau}_{1}]$, $r_{2}\in[0,\overline{\tau}_{2}]$ and $r_{3}\in[0,\overline{\tau}_{3}]$. Here,

[TABLE]

with $\varphi_{j}(v_{j},\overline{\tau}_{j}v_{j})=\Omega_{j}/\pi$ for $j=1,2,3$.

For the proof of Theorem 3.1 we need the following lemmas.

Lemma A.1**.**

Let $\varphi_{w}$ be defined as in (A.4). Then

[TABLE]

Proof..

By Proposition 11 in Owen (2005) we know that for f,g\in\BVHK we also have fg\in\BVHK; thus, we can consider each factor of $\varphi_{w}$ separately and the assertion follows by the one-dimensional theory for functions of bounded variation. ∎

Lemma A.2**.**

Let $h$ be a real-valued function on $[\bm{a},\bm{b}]$ and let $f$ be a function of bounded variation in the sense of Hardy and Krause on $[\bm{a},\bm{b}]$. Suppose that

[TABLE]

exists for all $\emptyset\neq v\subseteq 1:q$. Then there exists a constant $c\in\mathbb{R}$ such that

[TABLE]

Proof..

The assertion follows by an application of the integration by parts formula combined with the definition of the multivariate Riemann-Stieltjes integral and the variation in the sense of Hardy and Krause. ∎

Lemma A.3**.**

Let $\bm{a},\bm{b}\in\mathbb{R}^{q}$, $\widetilde{\bm{a}},\widetilde{\bm{b}}\in\mathbb{R}^{p}$ with $\bm{a}\leq\bm{b},\widetilde{\bm{a}}\leq\widetilde{\bm{b}}$, $p,q\in\mathbb{N}$. Let $\psi$ be a bounded function on $[\bm{a},\bm{b}]\times[\widetilde{\bm{a}},\widetilde{\bm{b}}]$ with

[TABLE]

Let $\{f_{n}\}$ be a sequence of functions such that $f_{n}\in D[\bm{a},\bm{b}]$ for all $n\in\mathbb{N}$. Suppose that $\lim_{n\to\infty}f_{n}=f$ in the Skorohod topology for a function $f\in C[\bm{a},\bm{b}]$. Moreover, suppose that

[TABLE]

exists for $\bm{s}_{v}\in[\bm{a}_{v},\bm{b}_{v}]$, $\bm{y}\in[\widetilde{\bm{a}},\widetilde{\bm{b}}]$ and for all $n\in\mathbb{N}$ and $\emptyset\neq v\subseteq 1:q$. Then,

[TABLE]

Proof..

An application of Lemma A.2 and the fact that Skorohod convergence towards a continuous limit implies uniform convergence show the assertion. ∎

We now prove Theorem 3.1. To simplify the notation we put $n=n_{1}=n_{2}=n_{3}$, but the proof can be completed along the same lines if the $n_{i},i=1,2,3,$ differ.

Proof of Theorem 3.1.

We consider the function space $D[0,1]^{3}$ and define the functional

[TABLE]

for $0<s_{0}\leq s\leq 1$, $t\in[0,\overline{\tau}_{1}]$, $r_{2}\in[0,\overline{\tau}_{2}]$, and $r_{3}\in[0,\overline{\tau}_{3}]$ on it, whenever the integral exists. Here $\varphi_{w}$ is defined as in (A.4). Let $\{f_{n}\}$ be a sequence of functions such that $f_{n}\in D[0,1]^{3}$ for all $n\in\mathbb{N}$. We know by Lemma A.3, that $\lim_{n\to\infty}f_{n}=f$ in the Skorohod topology with $f\in C[0,1]^{3}$ implies

[TABLE]

in $D([s_{0},1]\times[0,\overline{\tau}_{1}]\times[0,\overline{\tau}_{2}]\times[0,\overline{\tau}_{3}])$ as $n\to\infty$, since

[TABLE]

by Lemma A.1. Assumption 1 and the continuous mapping theorem now lead to

[TABLE]

as $n\to\infty$, since $P(B(v_{1},v_{2},v_{3})\in C\left([0,1]^{3}\right))=1$.

On the other hand we have

[TABLE]

We now interpret this last integral as a multivariate Riemann-Stieltjes integral. Let $\mathcal{Y}_{1}$ be a ladder on $[0,s]$ for $s\in[s_{0},1]$ and let $\mathcal{Y}_{2}$ and $\mathcal{Y}_{3}$ be ladders on $[0,1]$. Then we can write the triple integral as

[TABLE]

where $\bm{y}_{+}=\left(y^{+}_{1},y^{+}_{2},y^{+}_{3}\right)$ is the component-by-component successor of the point $\bm{y}=\left(y_{1},y_{2},y_{3}\right)$, $\widetilde{\bm{y}}=\left(\widetilde{y}^{1},\widetilde{y}^{2},\widetilde{y}^{3}\right)$ is an arbitrary point in the cube $[\bm{y},\bm{y}_{+}]$ and

[TABLE]

Consider now, without loss of generality, a ladder $\mathcal{Y}=\prod_{i=1}^{3}\mathcal{Y}_{i}$ with $y_{i}<k_{i}/n\leq y^{+}_{i}$ for all $k_{1}\in\{1,\ldots,\left\lfloor ns\right\rfloor\}$ and all $k_{2},k_{3}\in\{1,\ldots,n\}$, respectively, and write $\widetilde{y}_{k_{i}}$ for points $\widetilde{y}_{k_{i}}\in(y_{i},y_{i}^{+}],i=1,2,3$. As the floor function $\lfloor ny\rfloor$ is constant on intervals of the form $\left[(k-1)/n,k/n\right)$, $k\in\{1,\ldots,n\}$, the triple sums in the last expression can be combined and we obtain

[TABLE]

This leads to

[TABLE]

where the last equality holds as

[TABLE]

Now we can write (A.7) as

[TABLE]

If we recall the definition of $\varphi_{w}$ from (A.4) and the fact that $\tau_{j}\approx\overline{\tau}_{j}/n$, $j=1,2,3$, we obtain that

[TABLE]

This completes the proof. ∎

Proof of Lemma 3.1.

The proof of this lemma is omitted here as it is analogously to the proof of Corollary 1 in Pawlak and Steland (2013), since the properties in Theorem 5.1.4 in Ash and Gardner (1975) also hold true for multivariate Riemann-Stieltjes integrals. ∎

Proof of Lemma 3.2.

The assertion directly follows with Theorem 3.1 and the continuous mapping theorem. ∎

Before we prove Theorem 3.2 we review the Hwlaka-Koksma inequality which gives an error bound for the discrete approximation of a Riemann integral for functions of bounded variation in the sense of Hardy and Krause, see for example Niederreiter (1992), Theorem 2.11. For this we need the notion of discrepancy which measures the deviation of an arbitrary point set in the unit cube from a uniformly distributed point set. Hence, let $P$ be a point set consisting of $\bm{x}^{(1)},\ldots,\bm{x}^{(N)}\in[0,1]^{q}$. Now, for an arbitrary subset $B$ of $[0,1]^{q}$, we define

[TABLE]

i.e. $A(B;P)$ counts the number of points of $\bm{x}^{(1)},\ldots,\bm{x}^{(N)}$ lying in $B$. We then define the general discrepancy of the point set $P$ for a nonempty family $\mathcal{B}$ of Lebesgue-measurable subsets of $[0,1]^{q}$ as

[TABLE]

This leads to the following notion of discrepancy, cf. Definition 2.1. in Niederreiter (1992).

Definition A.1**.**

The star discrepancy $D_{N}^{\star}(P)=D_{N}^{\star}\left(\bm{x}^{(1)},\ldots,\bm{x}^{(N)}\right)$ of the point set $P$ is defined by $D_{N}^{\star}(P)\coloneqq D_{N}(\mathcal{J}^{\star};P)$, where $\mathcal{J}^{\star}$ is the family of all subintervals of $[0,1]^{q}$ of the form $\prod_{i=1}^{q}[0,u_{i})$.

We can now state the well-known Hwlaka-Koksma inequality for multivariable functions.

Lemma A.4**.**

If $f$ has bounded variation $V_{HK}(f)$ on $[0,1]^{q}$ in the sense of Hardy and Krause, then, for any $\bm{x}^{(1)},\ldots,\bm{x}^{(N)}\in[0,1]^{q}$, we have

[TABLE]

We now prove Theorem 3.2. To simplify the notation we put again $n=n_{1}=n_{2}=n_{3}$, but the proof can be completed along the same lines if the $n_{i},i=1,2,3,$ differ.

Proof of Theorem 3.2.

The local alternative is given by

[TABLE]

with

[TABLE]

for $\beta>0$. Our test statistic $\mathcal{F}_{n}(s,t,r_{2},r_{3})$ is therefore defined as

[TABLE]

$T_{n}^{(1)}(s,t,r_{2},r_{3})$ equals $\mathcal{F}_{n}(s,t,r_{2},r_{3})$ under the null hypothesis which converges to the process $\mathcal{F}(s,t,r_{2},r_{3})$ of Theorem 3.1 for $n\to\infty$.

By assumption on the sampling periods we obtain for the second process

[TABLE]

We now fix $\beta=3/2$ and set

[TABLE]

If we can show that

[TABLE]

tends to zero as $n\to\infty$ the assertion follows, since uniform convergence always implies convergence in the Skorohod topology.

If $\delta(t,r_{2},r_{3})$ is continuous we can proceed in an analogous way as in the proof of Theorem 3 of Pawlak and Steland (2013). Thus, it remains to treat the case that $\delta(t,r_{2},r_{3})$ is of bounded variation in the sense of Hardy and Krause. Our aim is to apply the Hwlaka-Koksma inequality of Lemma A.4. As this inequality is formulated for integrals over the unit cube we first observe that

[TABLE]

Put $\left(x_{k_{1}},x_{k_{2}},x_{k_{3}}\right)=\left(\frac{k_{1}\overline{\tau}_{1}}{n},\frac{k_{2}\overline{\tau}_{2}}{n},\frac{k_{3}\overline{\tau}_{3}}{n}\right)$ and $\left(\widetilde{x}_{k_{1}},\widetilde{x}_{k_{2}},\widetilde{x}_{k_{3}}\right)=\left(\frac{k_{1}}{ns},\frac{k_{2}}{n},\frac{k_{3}}{n}\right),$ for $\left(k_{1},k_{2},k_{3}\right)\in\{1,\ldots,\lfloor ns\rfloor\}\times\{1,\ldots,n\}\times\{1,\ldots,n\}$ and $s\geq s_{0}>0$. Then we obtain

[TABLE]

Thus we can reformulate (A.2) as

[TABLE]

An upper bound for this expression without the suprema is

[TABLE]

We first consider $S_{2}$. As both $\delta$ and $\varphi$ are of bounded variation in the sense of Hardy and Krause, $\varphi_{\delta}$ is bounded. Thus, for some constant $C\in\mathbb{R}$ we have

[TABLE]

which tends to zero as $n\to\infty$, uniformly for all $s\in[s_{0},1]$ and $(t,r_{2},r_{3})\in[0,\overline{\tau}_{1}]\times[0,\overline{\tau}_{2}]\times[0,\overline{\tau}_{3}]$.

To estimate $S_{1}$ we apply the Hwlaka-Koksma inequality. If we put $K\coloneqq\{1,\ldots,\left\lfloor ns\right\rfloor\}\times\{1,\ldots,n\}\times\{1,\ldots,n\}$ and

[TABLE]

where $N_{s}\coloneqq\lfloor ns\rfloor n^{2}$, we have

[TABLE]

By Proposition 11 in Owen (2005) we can estimate the variation by

[TABLE]

Since $\delta(t,r_{2},r_{3})$, by assumption, and $(z_{1},z_{2},z_{3})\mapsto\varphi_{\delta}\left(t,r_{2},r_{3},s\overline{\tau}_{1}z_{1},\overline{\tau}_{2}z_{2},\overline{\tau}_{3}z_{3}\right)$, by a similar argument as in Lemma A.1, are of bounded variation in the sense of Hardy and Krause uniformly in $t,r_{2},r_{3}$, we obtain

[TABLE]

It remains to verify that the discrepancy $D^{\star}_{N_{s}}\left(\widetilde{x}_{k_{1}},\widetilde{x}_{k_{2}},\widetilde{x}_{k_{3}}\right)$ is $o(1)$. As for arbitrary $u_{1},u_{2},u_{3}\in[0,1)$ we have $\lfloor nsu_{1}\rfloor$ points $\widetilde{x}_{k_{1}}$ with $\widetilde{x}_{k_{1}}<u_{1}$, $\lfloor nu_{2}\rfloor$ points $\widetilde{x}_{k_{2}}$ with $\widetilde{x}_{k_{2}}<u_{2}$, and $\lfloor nu_{3}\rfloor$ points $\widetilde{x}_{k_{3}}$ with $\widetilde{x}_{k_{3}}<u_{3}$ we obtain

[TABLE]

for appropiately chosen $\varepsilon,\varepsilon_{1},\varepsilon_{2},\varepsilon_{3}\in[0,1)$. This finally leads to

[TABLE]

uniformly for all $s\in[s_{0},1]$. Now, combining (A.11) and (A.12) with (A.10) it follows that $S_{1}\to 0$ as $n\to\infty$, uniformly in $s,t,r_{2},r_{3}$. Thus, assertion (A.2) also follows for the case that $\delta$ is a function of bounded variation in the sense of Hardy and Krause which finally completes the proof. ∎

Proof of Lemma 3.3.

The assertion directly follows with Theorem 3.2 and the continuous mapping theorem. ∎

Proof of Theorem 3.3.

The proof is analogoulsy to the one of Theorem 3.2. ∎

Proof of Theorem 4.1.

As $w$ fulfills the same assumptions as $\varphi$ one can adopt the proof of Theorem 3.1. ∎

Proof of Theorem 4.2.

The result easily follows from the fact that under Assumption 1 $y_{i,j,k}-\overline{y}_{\cdot,j,k}=\varepsilon_{i,j,k}-\overline{\varepsilon}_{\cdot,j,k}$ and thus

[TABLE]

as $\min_{1\leq i\leq 3}n_{i}\to\infty$. Now one may argue as in the proof of Theorem 3.1. ∎