A Wild Bootstrap for Degenerate Kernel Tests

Introduction

Statistical tests based on distribution embeddings into reproducing kernel Hilbert spaces have been applied in many contexts, including two sample testing , tests of independence , tests of conditional independence , and tests for higher order (Lancaster) interactions . For these tests, consistency is guaranteed if and only if the observations are independent and identically distributed. Much real-world data fails to satisfy the i.i.d. assumption: audio signals, EEG recordings, text documents, financial time series, and samples obtained when running Markov Chain Monte Carlo, all show significant temporal dependence patterns.

The asymptotic behaviour of kernel test statistics becomes quite different when temporal dependencies exist within the samples. In recent work on independence testing using the Hilbert-Schmidt Independence Criterion (HSIC) , the asymptotic distribution of the statistic under the null hypothesis is obtained for a pair of independent time series, which satisfy an absolute regularity or a $\phi$-mixing assumption. In this case, the null distribution is shown to be an infinite weighted sum of dependent $\chi^{2}$-variables, as opposed to the sum of independent $\chi^{2}$-variables obtained in the i.i.d. setting . The difference in the asymptotic null distributions has important implications in practice: under the i.i.d. assumption, an empirical estimate of the null distribution can be obtained by repeatedly permuting the time indices of one of the signals. This breaks the temporal dependence within the permuted signal, which causes the test to return an elevated number of false positives, when used for testing time series. To address this problem, an alternative estimate of the null distribution is proposed in , where the null distribution is simulated by repeatedly shifting one signal relative to the other. This preserves the temporal structure within each signal, while breaking the cross-signal dependence.

A serious limitation of the shift procedure in is that it is specific to the problem of independence testing: there is no obvious way to generalise it to other testing contexts. For instance, we might have two time series, with the goal of comparing their marginal distributions - this is a generalization of the two-sample setting to which the shift approach does not apply.

We note, however, that many kernel tests have a test statistic with a particular structure: the Maximum Mean Discrepancy (MMD), HSIC, and the Lancaster interaction statistic, each have empirical estimates which can be cast as normalized $V$-statistics, $\frac{1}{n^{m-1}}\sum_{1\leq i_{1},...,i_{m}\leq n}h(Z_{i_{1}},...,Z_{i_{m}})$, where $Z_{i_{1}},...,Z_{i_{m}}$ are samples from a random process at the time points $\{i_{1},\ldots,i_{m}\}$. We show that a method of external randomization known as the wild bootstrap may be applied to simulate from the null distribution. In brief, the arguments of the above sum are repeatedly multiplied by random, user-defined time series. For a test of level $\alpha$, the $1-\alpha$ quantile of the empirical distribution obtained using these perturbed statistics serves as the test threshold. This approach has the important advantage over that it may be applied to all kernel-based tests for which $V$-statistics are employed, and not just for independence tests.

The main result of this paper is to show that the wild bootstrap procedure yields consistent tests for time series, i.e., tests based on the wild bootstrap have a Type I error rate (of wrongly rejecting the null hypothesis) approaching the design parameter $\alpha$, and a Type II error (of wrongly accepting the null) approaching zero, as the number of samples increases. We use this result to construct a two-sample test using MMD, and an independence test using HSIC. The latter procedure is applied both to testing for instantaneous independence, and to testing for independence across multiple time lags, for which the earlier shift procedure of cannot be applied.

We begin our presentation in Section 2, with a review of the $\tau$-mixing assumption required of the time series, as well as of $V$-statistics (of which MMD and HSIC are instances). We also introduce the form taken by the wild bootstrap. In Section 3, we establish a general consistency result for the wild bootstrap procedure on $V$-statistics, which we apply to MMD and to HSIC in Section 4. Finally, in Section 5, we present a number of empirical comparisons: in the two sample case, we test for differences in audio signals with the same underlying pitch, and present a performance diagnostic for the output of a Gibbs sampler (the MCMC M.D.); in the independence case, we test for independence of two time series sharing a common variance (a characteristic of econometric models), and compare against the test of in the case where dependence may occur at multiple, potentially unknown lags. Our tests outperform both the naive approach which neglects the dependence structure within the samples, and the approach of , when testing across multiple lags. Detailed proofs are found in the appendices.

Background

The main results of the paper are based around two concepts: $\tau$-mixing , which describes the dependence within the time series, and $V$-statistics , which constitute our test statistics. In this section, we review these topics, and introduce the concept of wild bootstrapped $V$-statistics, which will be the key ingredient in our test construction.

and $\Lambda$ is the set of all one-Lipschitz continuous real-valued functions on the domain of $X$. $\tau(\mathcal{M},X)$ can be interpreted as the minimal $L_{1}$ distance between $X$ and $X^{*}$ such that $X\overset{d}{=}X^{*}$ and $X^{*}$ is independent of $\mathcal{M}\subset\mathcal{F}$. Furthermore, if $\mathcal{F}$ is rich enough, this $X^{*}$ can be constructed.

Note that this mixing definition differs from commonly used notion of $\beta$ mixing (or $\phi$ mixing), which was required in the previous work . We describe in more detail how these notions of dependence are related in Appendix B.

V𝑉V-statistics.

The test statistics considered in this paper are always $V$-statistics. Given the observations $Z=\left\{Z_{t}\right\}_{t=1}^{n}$, a $V$-statistic of a symmetric function $h$ taking $m$ arguments is given by

where $N^{m}$ is a Cartesian power of a set $N=\{1,...,n\}$. For simplicity, we will often drop the second argument and write simply $V(h)$.

We will refer to the function $h$ as to the core of the $V$-statistic $V(h)$. While such functions are usually called kernels in the literature, in this paper we reserve the term kernel for positive-definite functions taking two arguments. A core $h$ is said to be $j$-degenerate if for each $z_{1},\ldots,z_{j}$ $Eh(z_{1},\ldots,z_{j},Z_{j+1}^{*},\ldots,Z_{m}^{*})=0,$ where $Z_{j+1}^{*},\ldots,Z_{m}^{*}$ are independent copies of $Z_{1}$. If $h$ is $j$-degenerate for all $j\leq m-1$, we will say that it is canonical. For a one-degenerate core $h$, we define an auxiliary function $h_{2}$, called the second component of the core, and given by $h_{2}(z_{1},z_{2})=Eh(z_{1},z_{2},Z_{3}^{*},\ldots,Z_{m}^{*}).$ Finally we say that $nV(h)$ is a normalized $V$-statistic, and that a $V$-statistic with a one-degenerate core is a degenerate $V$-statistic. This degeneracy is common to many kernel statistics when the null hypothesis holds .

Our main results will rely on the fact that $h_{2}$ governs the asymptotic behaviour of normalized degenerate $V$-statistics. Unfortunately, the limiting distribution of such $V$-statistics is quite complicated - it is an infinite sum of dependent $\chi^{2}$-distributed random variables, with a dependence determined by the temporal dependence structure within the process $\{Z_{t}\}$ and by the eigenfunctions of a certain integral operator associated with $h_{2}$ . Therefore, we propose a bootstrapped version of the $V$-statistics which will allow a consistent approximation of this difficult limiting distribution.

Bootstrapped V𝑉V-statistic.

We will study two versions of the bootstrapped $V$-statistics

As noted in in [19, Remark 2], a simple realization of a process that satisfies this assumption is $W_{t,n}=e^{-1/l_{n}}W_{t-1,n}+\sqrt{1-e^{-2/l_{n}}}\epsilon_{t}$ where $W_{0,n}$ and $\epsilon_{1},\ldots,\epsilon_{n}$ are independent standard normal random variables. For simplicity, we will drop the index $n$ and write $W_{t}$ instead of $W_{t,n}$. A process that fulfils the bootstrap assumption will be called bootstrap process. Further discussion of the wild bootstrap is provided in the Appendix A. The versions of the bootstrapped $V$-statistics in (2) and (3) were previously studied in for the case of canonical cores of degree $m=2$. We extend their results to higher degree cores (common within the kernel testing framework), which are not necessarily one-degenerate. When stating a fact that applies to both $B_{1}$ and $B_{2}$, we will simply write $B$, and the argument $h$ or index $n$ will be dropped when there is no ambiguity.

Asymptotics of wild bootstrapped V𝑉V-statistics

In this section, we present main Theorems that describe asymptotic behaviour of $V$-statistics, all the poofs are in the Appendix C. In the next section, these results will be used to construct kernel-based statistical tests applicable to dependent observations. Tests are constructed so that the $V$-statistic is degenerate under the null hypothesis and non-degenerate under the alternative. Theorem 1 guarantees that the bootstrapped $V$-statistic will converge to the same limiting null distribution as the simple $V$-statistic.

Throughout this paper we will make one mild assumption

where $Z_{i}=(Z_{i_{1}},\cdots Z_{i_{m}})$. This assumption is almost always automatically satisfied, since most of the kernels used in practice are bounded.

Assume that the stationary process $Z_{t}$ is $\tau$-dependent with $\sum_{r=1}^{\infty}r^{2}\sqrt{\tau(r)}<\infty$. If the core $h$ is a Lipschitz continuous, one-degenerate and its $h_{2}$-component is a positive definite kernel, such that $Eh_{2}(Z_{0},Z_{0})<\infty$, then $nB_{n}$ (2), (3), and $nV_{n}$ (1) converge weakly to the same distribution $V$. Moreover $nB_{n}(h_{2})$ and $nV_{n}(h_{2})$ converge weakly to $\binom{m}{2}^{-1}V$.

On the other hand, if the $V$-statistic is not degenerate, which is usually true under the alternative, it converges to some non-zero constant.

Assume that the stationary process $Z_{t}$ is $\tau$-dependent with $\tau(r)=o(r^{-4})$. If the core $h$ is a Lipschitz continuous, and $h_{0}$ component is positive then $V_{n}$ converges in mean squared to $h_{0}$.

In this setting, Theorem 3 guarantees that the bootstrapped $V$-statistic will converge to zero in probability. This property is necessary in testing, as it implies that the test thresholds computed using the bootstrapped $V$-statistics will also converge to zero, and so will the corresponding Type II error.

Assume that the stationary process $\{Z_{t}\}$ is $\tau$-dependent with a coefficient $\tau(r)=o(r^{-4})$. If the core $h$ is a function of $m>1$ arguments then $B_{1}(h)$ and $o(n)B_{2}(h)$ converge to zero in mean squared.

Although both $B_{2}$ and $B_{1}$ converge to zero, the rate does not seem to be that same. As a consequence, tests that utilize $B_{2}$ usually give lower Type II error then the ones that use $B_{1}$. On the other hand, $B_{1}$ seems to better approximate $V$-statistic distribution under the null hypothesis. This agrees with our experiments in Section 5 as well as with those in [19, Section 5]). These results a sufficient for adopting kernel tests developed for i.i.d. data to tests that work on random processes. In particular Theorem 1 justifies usage of bootstraped $V$-statistics for estimating quantiles of the null distribution, while Theorems 23 guarantee consistency.

Calculate the test statistic $nV_{n}(h)$.

Obtain wild bootstrap samples $\{B_{n}(h)\}_{i=1}^{D}$ and estimate the $1-\alpha$ empirical quantile of these samples.

If $nV_{n}(h)$ exceeds the quantile, reject.

Applications to Kernel Tests

In this section, we describe how the wild bootstrap for $V$-statistics can be used to construct kernel tests for independence and the two-sample problem, which are applicable to weakly dependent observations. We start by reviewing the main concepts underpinning the kernel testing framework.

In this section, we describe how the wild bootstrap for $V$-statistics can be used to construct kernel tests for independence and the two-sample problem, in presence of weakly dependent observations.

Denote the observations by $\{X_{i}\}_{i=1}^{n_{x}}\sim P_{x}$, and $\{Y_{j}\}_{j=1}^{n_{y}}\sim P_{y}$. Our goal is to test the null hypothesis $\mathbf{H}_{0}:P_{x}=P_{y}$ vs. the alternative $\mathbf{H}_{1}:P_{x}\neq P_{y}$. In the case where samples have equal sizes, i.e., $n_{x}=n_{y}$, application of the wild bootstrap to MMD-based tests on dependent samples is straightforward: the empirical MMD can be written as a $V$-statistic with the core of degree two on pairs $z_{i}=(x_{i},y_{i})$ given by $h(z_{1},z_{2})=k(x_{1},x_{2})-k(x_{1},y_{2})-k(x_{2},y_{1})+k(y_{1},y_{2})$. It is clear that whenever $k$ is Lipschitz continuous and bounded, so is $h$. Moreover, $h$ is a valid positive definite kernel, since it can be represented as an RKHS inner product $\left\langle k(\cdot,x_{1})-k(\cdot,y_{1}),k(\cdot,x_{2})-k(\cdot,y_{2})\right\rangle_{\mathcal{H}_{k}}$. Under the null hypothesis, $h$ is also one-degenerate, i.e., $Eh\left((x_{1},y_{1}),(X_{2},Y_{2})\right)=0$. Therefore, we can use the bootstrapped statistics in (2) and (3) to approximate the null distribution and attain a desired test level.

When $n_{x}\neq n_{y}$, however, it is no longer possible to write the empirical MMD as a one-sample $V$-statistic. We will therefore require the following bootstrapped version of MMD

Let $k$ be bounded and Lipschitz continuous, and let $\left\{X_{t}\right\}$ and $\left\{Y_{t}\right\}$ both be $\tau$-dependent with coefficients $\tau(r)=O(r^{-6-\epsilon})$, but independent of each other. Further, let $n_{x}=\rho_{x}n$ and $n_{y}=\rho_{y}n$ where $n=n_{x}+n_{y}$. Then, under the null hypothesis $P_{x}=P_{y}$, $\varphi\left(\rho_{x}\rho_{y}n\widehat{\text{MMD}}_{k},\rho_{x}\rho_{y}n\widehat{\text{MMD}}_{k,b}\right)\to 0$ in probability as $n\to\infty$, where $\varphi$ is the Prokhorov metric and $\widehat{\text{MMD}}_{k}$ is the MMD between empirical measures.

2 Wild Bootstrap For HSIC

Using HSIC in the context of random processes is not new in the machine learning literature. For a 1-approximating functional of an absolutely regular process , convergence in probability of the empirical HSIC to its population value was shown in . No asymptotic distributions were obtained, however, nor was a statistical test constructed. The asymptotics of a normalized $V$-statistic were obtained in for absolutely regular and $\phi$-mixing processes . Due to the intractability of the null distribution for the test statistic, the authors propose a procedure to approximate its null distribution using circular shifts of the observations leading to tests of instantaneous independence, i.e., of X_{t}\rotatebox[origin={c}]{90.0}{\models}Y_{t}, $\forall t$. This was shown to be consistent under the null (i.e., leading to the correct Type I error), however consistency of the shift procedure under the alternative is a challenging open question (see [8, Section A.2] for further discussion). In contrast, the wild bootstrap guarantees test consistency under both hypotheses: null and alternative, which is a major advantage. In addition, the wild bootstrap can be used in constructing a test for the harder problem of determining independence across multiple lags simultaneously, similar to the one in .

Following symmetrisation, it is shown in that the empirical HSIC can be written as a degree four $V$-statistic with core given by

where we denote by $S_{n}$ the group of permutations over $n$ elements. One-degeneracy of the core under the null hypothesis was stated in [15, Theorem 2], [15, Section A.2, following eq. (11)] shows that $h_{2}$ is a kernel; $h_{0}\geq 0$ follows from the fact that HSIC is a distance. Using Theorems 1,3,2 we can construct an independence test using $h$. Drawback of this test, when implemented in the most straightforward way, is its quadruple computational complexity. To achieve quadratic time complexity, that matches time complexity of HSIC test for i.i.d. data, we modify our bootstrapped statistic.

[15, section A.2, following eq. (11)] The second component of $h$ is $h_{2}(z_{1},z_{2})=\frac{1}{6}\bar{k}(x_{1},x_{2})\bar{l}(y_{1},y_{2}).$

Squared norm of $T_{n}$ is equal to $6B(h_{2})$.

Next we relate $S_{n}$ to $T_{n}$ – we show that the difference between them is asymptotically negligible. We start with a technical lemma.

If $(\bar{k}\times\bar{k},Z_{i})$ is of type $\varDelta$ of order $O(r^{-4})$ (see Definition 2), then

Since $(\bar{k}\times\bar{k},X_{i})$ is of type $\varDelta$, by Lemma 4, the expected value is of order $O(1)$. ∎

If $(\bar{k}\times\bar{k},Z_{i})$, $(\bar{l}\times\bar{l},Z_{i})$ are of type $\varDelta$ of order $O(r^{-4})$, then, under the null, $\|S_{n}\|^{2}-\|T_{n}\|^{2}$ converges to zero in mean square. Under the alternative $\frac{1}{n}(\|S_{n}\|^{2}-\|T_{n}\|^{2})$ converges to zero in mean square.

We first show that $E\|S_{n}-T_{n}\|^{2}\to 0$ both under the null and the alternative. Then, using the fact that $\|T_{n}\|^{2}<\infty$ under the null and $\frac{1}{n}\|T_{n}\|^{2}<\infty$ under alternative we will conclude the proof. The difference $S_{n}-T_{n}$ is

We examine differences separately – it is sufficient to show that each difference converges to zero in mean square.

The expected norm of the first difference is

We used $\|v\otimes u\|=\|v\|\|u\|$ and Cauchy-Schwarz inequality. By Lemma 2 the first term is $O(1)$. The second term is equal to

The expected value converges to zero in mean square by Lemma 4 (the assumption $\sup_{i,j}k(X_{i},X_{j})<\infty$ is satisfied). Using similar reasoning, the second term

also converges to zero. The final term is

$\frac{1}{n}\sum_{i\in N}Q_{i}$ converges in mean square to zero (Lemmas 6, 7). We rewrite the second term

To show that the above expression converges to zero it is sufficient to show that $E\|T_{n}\|^{2}<\infty$ and $E\|S_{n}\|^{2}<\infty$. Under the null hypothesis, by Lemma 4, expected value of $E\|T_{n}\|^{2}=nB_{n}(h_{2})$ is finite. Since $E\|T_{n}-S_{n}\|^{2}\to 0$ we also have $E\|T_{n}-S_{n}\|\to 0$. Therefore we have

it is sufficient to show that $n^{-1}E\|T_{n}\|^{2}<\infty$ and $n^{-1}E\|S_{n}\|^{2}<\infty$. By Theorem 3, $n^{-1}E\|T_{n}\|^{2}<\infty$ is finite and, using the reasoning similar to the one above, we have that $n^{-1}E\|S_{n}\|^{2}<\infty$. ∎

This shows that we can use squared norm of $S_{n}$ as a bootstrapped test statistic. For HSIC we redefine $B_{n}$

Finally, notice that both statistics 7 and 8 can be calculated in quadratic time.

Let $Z_{t}=\left(X_{t},Y_{t}\right)$ be a stationary process that is $\tau$-dependent such that $\sum_{r=1}^{\infty}r^{2}\sqrt{\tau(r)}<\infty$. Under the null hypothesis $B_{n}^{*}$ (7) and $nV_{n}(h)$ (8)converge weakly to the same distribution. Under the alternative hypothesis $B_{n}^{*}$ converges to zero in probability, while $V_{n}(h)$ converges to a positive constant.

By Lemma 1, $6nB_{n}(h_{2})=\|T_{n}\|^{2}$. By definition (7), $B_{n}^{*}=\|S_{n}\|^{2}$ . By Lemma 3, $6nB_{n}(h_{2})-B_{n}^{*}$ converges to zero in mean square. We check assumptions; since process $Z_{t}$ is $\tau$-mixing (of order $o(r^{-4})$ ) and both $\bar{k}$, $\bar{l}$ are canonical, Lemma 11 guarantees that $(\bar{k},Z_{i})$, $(\bar{l},Z_{i})$ are of type $\varDelta$ of order $O(r^{-4})$.

Under the null hypothesis, by Theorem 1, $nV_{n}(h)-6nB_{n}(h_{2})$ converges to zero. We check assumptions; by Lemma 1, $h_{2}$ is a symmetric, one-degenerate, bounded kernel, assumptions concerning $\tau$-mixing are satisfied.

Under the alternative, by Theorem 3 and Lemma 3 respectively, $6B_{n}(h_{2})$ and $\frac{1}{n}B_{n}^{*}-6B_{n}(h_{2})$ converge to zero in mean square. By Theorem 3, $V_{n}(h)$ converges to a positive constant. ∎

We consider two types of tests: instantaneous independence and independence at multiple time lags.

Test of instantaneous independence

Here, the null hypothesis $\mathbf{H_{0}}$ is that $X_{t}$ and $Y_{t}$ are independent at all times $t$, and the alternative hypothesis $\mathbf{H_{1}}$ is that they are dependent. We use Proposition2 directly to bootstrap an appropriate quantile and compare it with a test statistic.

Lag-HSIC

Proposition 2 allows us to construct a test of time series independence that is similar to one designed by . Here, we will be testing against a broader null hypothesis: $X_{t}$ and $Y_{t^{\prime}}$ are independent for $|t-t^{\prime}|<M$ for an arbitrary large but fixed $M$.

Since the time series $Z_{t}=(X_{t},Y_{t})$ is stationary, it suffices to check whether there exists a dependency between $X_{t}$ and $Y_{t+m}$ for $-M\leq m\leq M$. Since each lag corresponds to an individual hypothesis, we will require a Bonferroni correction to attain a desired test level $\alpha$. We therefore define $q=1-\frac{\alpha}{2M+1}$. The shifted time series will be denoted $Z_{t}^{m}=(X_{t},Y_{t+m})$. Let $S_{m,n}=nV_{n}(h,Z^{m})$ denote the value of the normalized HSIC statistic calculated on the shifted process $Z_{t}^{m}$. Let $F_{b,n}$ denote the empirical cumulative distribution function obtained by the bootstrap procedure using $B_{n}^{*}$ (7). The test will then reject the null hypothesis if the event $\mathcal{A}_{n}=\left\{\max_{-M\leq m\leq M}S_{m,n}>F^{-1}_{b,n}(q)\right\}$ occurs. By a simple application of the union bound, it is clear that the asymptotic probability of the Type I error will be $\lim_{n\to\infty}P_{\,\mathbf{H_{0}}}\left(\mathcal{A}_{n}\right)\leq\alpha$. On the other hand, if the alternative holds, there exists some $m$ with $|m|\leq M$ for which $V_{n}(h,Z^{m})=n^{-1}S_{m,n}$ converges to a non-zero constant. In this case

as long as $n^{-1}F^{-1}_{b,n}(q)\to 0$, which follows from the convergence of $B_{n}^{*}$ (7) to zero in probability shown in Proposition 2. Therefore, the Type II error of the multiple lag test is guaranteed to converge to zero as the sample size increases. Our experiments in the next Section demonstrate that while this procedure is defined over a finite range of lags, it results in tests more powerful than the procedure for an infinite number of lags proposed in . We note that a procedure that works for an infinite number of lags, although possible to construct, does not add much practical value under the present assumptions. Indeed, since the $\tau$-mixing assumption applies to the joint sequence $Z_{t}=(X_{t},Y_{t})$, dependence between $X_{t}$ and $Y_{t+m}$ is bound to disappear at a rate of $o(m^{-6})$, i.e., the variables both within and across the two series are assumed to become gradually independent at large lags.

Experiments

We employ MMD in order to diagnose how far an MCMC chain is from its stationary distribution [23, Section 5], by comparing the MCMC sample to a benchmark sample. A hypothesis test of whether the sampler has converged based on the standard permutation-based bootstrap leads to too many rejections of the null hypothesis, due to dependence within the chain. Thus, one would require heavily thinned chains, which is wasteful of samples and computationally burdensome. Our experiments indicate that the wild bootstrap approach allows consistent tests directly on the chains, as it attains a desired number of false positives. To assess performance of the wild bootstrap in determining MCMC convergence, we consider the situation where samples $\{X_{i}\}$ and $\{Y_{i}\}$ are bivariate, and both have the identical marginal distribution given by an elongated normal P=\mathcal{N}\left(\left[\begin{array}[]{cc}0&0\end{array}\right],\left[\begin{array}[]{cc}15.5&14.5\\ 14.5&15.5\end{array}\right]\right). However, they could have arisen either as independent samples, or as outputs of the Gibbs sampler with stationary distribution $P$. Table 1 shows the rejection rates under the significance level $\alpha=0.05$. It is clear that in the case where at least one of the samples is a Gibbs chain, the permutation-based test has a Type I error much larger than $\alpha$. The wild bootstrap using $V_{b1}$ (without artificial degeneration) yields the correct Type I error control in these cases. Consistent with findings in [19, Section 5], $V_{b1}$ mimics the null distribution better than $V_{b2}$. The bootstrapped statistic $\widehat{\text{MMD}}_{k,b}$ in (4) which also relies on the artificially degenerated bootstrap processes, behaves similarly to $V_{b2}$. In the alternative scenario where $\{Y_{i}\}$ was taken from a distribution with the same covariance structure but with the mean set to \mu=\left[\begin{array}[]{cc}2.5&0\end{array}\right], the Type II error for all tests was zero.

Pitch-evoking sounds

Instantaneous independence

To examine instantaneous independence test performance, we compare it with the Shift-HSIC procedure on the ’Extinct Gaussian’ autoregressive process proposed in the [8, Section 4.1]. Using exactly the same setting we compute type I error as a function of the temporal dependence and type II error as a function of extinction rate. Figure 1 shows that all three tests (Shift-HSIC and tests based on $V_{b1}$ and $V_{b2}$) perform similarly.

Lag-HSIC

The KCSD is, to our knowledge, the only test procedure to reject the null hypothesis if there exist $t$,$t^{\prime}$ such that $Z_{t}$ and $Z_{t^{\prime}}$ are dependent. In the experiments, we compare lag-HSIC with KCSD on two kinds of processes: one inspired by econometrics and one from . In lag-HSIC, the number of lags under examination was equal to $\max\{10,\log n\}$, where $n$ is the sample size. We used Gaussian kernels with widths estimated by the median heuristic. The cumulative distribution of the $V$-statistics was approximated by samples from $nV_{b2}$. To model the tail of this distribution, we have fitted the generalized Pareto distribution to the bootstrapped samples ( shows that for a large class of underlying distribution functions such an approximation is valid). The first process is a pair of two time series which share a common variance,

The above set of equations is an instance of the VEC dynamics used in econometrics to model market volatility. The left panel of the Figure 2 presents the Type II error rate: for KCSD it remains at 90% while for lag-HSIC it gradually drops to zero. The Type I error, which we calculated by sampling two independent copies $(X^{(1)}_{t},Y^{(1)}_{t})$ and $(X^{(2)}_{t},Y^{(2)}_{t})$ of the process and performing the tests on the pair $(X^{(1)}_{t},Y^{(2)}_{t})$, was around 5% for both of the tests. Our next experiment is a process sampled according to the dynamics proposed by ,

with parameters $C=.4$, $f_{1}=4Hz$,$f_{2}=20Hz$, and frequency $\frac{1}{T_{s}}=100Hz$. We compared performance of the KCSD algorithm, with parameters set to vales recommended in , and the lag-HSIC algorithm. The Type II error of lag-HSIC, presented in the right panel of the Figure 2, is substantially lower than that of KCSD. The Type I error ($C=0$) is equal or lower than 5% for both procedures. Most oddly, KCSD error seems to converge to zero in steps. This may be due to the method relying on a spectral decomposition of the signals across a fixed set of bands. As the number of samples increases, the quality of the spectrogram will improve, and dependence will become apparent in bands where it was undetectable at shorter signal lengths.

References

Appendix A An Introduction to the Wild Bootstrap

Bootstrap methods aim to evaluate the accuracy of the sample estimates - they are particularly useful when dealing with complicated distributions, or when the assumptions of a parametric procedure are in doubt. Bootstrap methods randomize the dataset used for the sample estimate calculation, so that a new dataset with a similar statistical properties is obtained, e.g. one popular method is resampling. In the wild bootstrap method the observations in the dataset are multiplied by appropriate random numbers. To present the idea behind the wild bootstrap we will discuss a toy example similar to that in , and then relate it to the wild bootstrap method used in this article.

Consider a stationary autoregressive-moving-average random process $\{X_{i}\}_{i\in\mathbf{Z}}$ with zero mean. The normalized sample mean of the process $X_{t}$ has normal distribution

where $\sigma_{\infty}^{2}=\sum_{j=-\infty}^{j=\infty}cov(X_{0},X_{j})$. The variance $\sigma_{\infty}^{2}$ is not easy to estimate (the naive approach of approximating different covariances separately and summing them has several drawbacks, e.g. how many empirical covariances should be calculated?). Using the wild bootstrap method we will construct processes that mimic behaviour of the $X_{t}$ process and then use them to approximate the distribution of the normalized sample mean, $\frac{\sum_{i=1}^{N}X_{i}}{\sqrt{n}}$. The bootstrap process used to to randomize the sample meets the following criteria:

$\{W_{t,n}\}_{1\leq t\leq n}$ is a row-wise strictly stationary triangular array independent of all $X_{t}$, such that $EW_{t,n}=0$ and $\sup_{n}E|W_{t,n}^{2+\sigma}|<\infty$ for some $\sigma>0$.

The autocovariance of the process is given by $EW_{s,n}W_{t,n}=\rho(|s-t|/l_{n})$ for some function $\rho$, such that $\lim_{u\to 0}\rho(u)=1$.

The sequence $\left\{l_{n}\right\}$ is taken such that $\lim_{n\to\infty}l_{n}=\infty$.

A process that fulfils those criteria, given also in the main article, is

We need to show that the distribution of the normalized sample mean of the process $Y_{t}^{n}=W_{t}^{n}X_{t}$, where $|t|\leq n$, mimics the distribution $N(0,\sigma_{\infty}^{2})$. It suffices to calculate the expected value and correlations:

The asymptotic auto-covariance structure of the process $Y_{t}$ is the same as the auto-covariance structure of the process $X_{t}$. Therefore

This mechanism is used in . Recall that, under some assumptions, a normalized V-statistic can be written as

where $\lambda_{k}$ are eigenvalues and $\phi_{k}$ are eigenfunction of the kernel $h$, respectively. Since $E\phi_{k}(X_{i})=0$ (degeneracy condition) one may replace

and conclude, as in the toy example, that the limiting distribution of the single component of the sum $\sum_{k}\lambda_{k}...$ remains the same. One of the main contributions of is in showing that the distribution of the whole sum $\sum_{k}\lambda_{k}\left(\frac{\sum_{i=1}^{n}W_{t}^{n}\phi_{k}(X_{i})}{\sqrt{n}}\right)^{2}$ with the components bootstrapped converges in a particular sense (in probability in Prokhorov metric) to the distribution of the normalized V-statistic, $\frac{1}{n}\sum_{1\leq i,j\leq n}h(X_{i},X_{j})$.

Appendix B Relation between β𝛽\beta,ϕitalic-ϕ\phi and τ𝜏\tau mixing

Strong mixing is historically the most studied type of temporal dependence – a lot of models, example being Markov Chains, are proved to be strongly mixing, therefore it’s useful to relate weak mixing to strong mixing. For a random variable $X$ on a probability space $(\Omega,\mathcal{F},P_{X})$ and $\mathcal{M}\subset\mathcal{F}$ we define

A process is called $\beta$-mixing or absolutely regular if

By we have $\beta(\mathcal{M},\sigma(X))\leq\phi(\mathcal{M},\sigma(X))$ .

[10, Equation 7.6] relates $\tau$-mixing and $\beta$-mixing , as follows: if $Q_{x}$ is the generalized inverse of the tail function

While this definition can be hard to interpret, it can be simplified in the case $E|X|^{p}=M$ for some $p>1$, since via Markov’s inequality $P(|X|>t)\leq\frac{M}{t^{p}}$, and thus $\frac{M}{t^{p}}\leq u$ implies $P(|X|>t)\leq u$. Therefore $Q^{\prime}(u)=\frac{M}{\sqrt[p]{u}}\geq Q_{x}(u)$. As a result, we have the following inequality

provides examples of systems that are tau-mixing. In particular, given that certain assumptions are satisfied causal functions of stationary sequences, iterated random functions, Markov chains, expanding maps are all $\tau$-mixing.

Of particular interest to this work are Markov chains. The assumptions provided by , under which Markov chains are tau-mixing are somehow difficult to check but we can use classical theorems about the $\beta$-mixing). In particular [7, Corollary 3.6] states that a Harris recurrent (chain returns to a fixed set of the state space an infinite number of times) and aperiodic Markov chain satisfies absolute regularity. [7, Theorem 3.7] states that geometric ergodicity $\forall_{x}\|P^{n}(x,\cdot)-\pi\|_{TV}\leq Cq^{n},0<q<1$ implies geometric decay of the $\beta$ coefficient. Interestingly [7, Theorem 3.3] describes situations in which a non-stationary chain $\beta$-mixes exponentially.

Using inequality 18 between $\tau$-mixing coefficient and strong mixing coefficients one can use those classical theorems show that e.g for $p=2$ we have

Appendix C Proofs

In this section we prove the main theorems. As for the notation, $n$ denotes number of observations, $N=\{1,\cdots,n\}$, if $h$ is function then $h\times h$ denotes a product of $h$ with itself, $\lim_{n\to\infty}X_{n}\overset{L_{2}}{=}X$ denotes convergence in mean square

Hoeffding decomposition reduces any $V$-statistic to a sum of canonical $V$-statistics with canonical cores $h_{c}$, which are easier to study in context of non-iid data. As an illustration, consider a canonical core $h$ of $m$ arguments and fix some indexes $i_{1}\leq\cdots\leq i_{m-1}\ll i_{m}$, for a sake of example we may assume that indexes represent time. If observations $Z_{i_{1}},\cdots,Z_{i_{m-1}}$ are independent of the observation $Z_{i_{m}}$, then the expected value of $h(Z_{i_{1}},\cdots,Z_{i_{m}})$, by degeneracy, is equal to zero. If it is reasonable to assume that $Z_{i_{m}}$ is almost independent of $Z_{i_{1}},\cdots,Z_{i_{m-1}}$, maybe because it is so distant in time, then it is also reasonable to expect that for a canonical core $h$ (which is not too complicated )

which follows from the following approximate calculation

Associate with any set of indexes $i_{1},\cdots,i_{m}$ its nearest neighbor within the set. Suppose $i_{r}$ is is an index with the most distant nearest neighbor. We will call $i_{r}$ the most isolated index, and we will refer to its distance to the nearest neighbor as an isolation distance.

Consider a following example, for the set $\{1,5,7\}$, $1$ is the most isolated index and the isolation distance is $4$.

Given a sequence of random variables $Z_{t}$ and a function $h$, if for all sets of indexes $i_{1},\cdots,i_{m}$, with the isolation distance equal to $r$

for some some function $\varDelta$, then we say that the pair $(h,Z_{t})$ is of type $\varDelta$.

The next theorem shows a growth rate of a canonical $V$-statistic when a pair $h,Z_{t}$ is of type $\varDelta$.

Let $(Z_{t},h)$, where $h$ is a function of $m>1$ arguments, be a of type $\varDelta$, with $\varDelta(h,r)=o(r^{-k})$ for some $k$, then

The proof uses a technique similar to [1, Lemma 3]. We will focus on ordered $m$-tuples $1\leq i_{1}\leq\ldots\leq i_{m}\leq n$, and by considering all possible permutations of their indices, we obtain an upper bound

where (strict) inequality stems from the fact that the $m$-tuples with some coinciding entries appear multiple times on the right.

Since $(h,Z_{t})$ is a of type $\varDelta$

where $w(i)$ is an isolating distance of the index set $i=i_{1},\cdots i_{m}$. We need to estimate order of the sum

Let us upper bound the number of ordered $m$-tuples $i$ with $w(i)=w$. Denote $s=\left\lfloor\frac{m}{2}\right\rfloor+1$. $i_{1}$ can take $n$ different values, but since $i_{2}\leq i_{1}+w$, $i_{2}$ can take at most $w+1$ different values. For $2\leq l\leq s-1$, since $\min\left\{i_{2l}-i_{2l-1},i_{2l-1}-i_{2l-2}\right\}\leq w$, we can either let $i_{2l-1}$ take up to $n$ different values and let $i_{2l}$ take up to $w+1$ different values (if $i_{2l}-i_{2l-1}\leq i_{2l-1}-i_{2l-2}$) or let $i_{2l-1}$ take up to $w+1$ different values and let $i_{2l}$ take up to $n$ different values (if $i_{2l}-i_{2l-1}>i_{2l-1}-i_{2l-2}$), upper bounding the total number of choices for $\left[i_{2l-1},i_{2l}\right]$ by $2n(w+1)$. Finally, the last term $i_{m}$ can always have at most $w+1$ different values. This brings the total number of $m$-tuples with $w(i)=w$ to at most $2^{s-2}n^{s-1}(w+1)^{s}$. Thus, the number of $m$-tuples with $w(i)=0$ is $O(n^{s-1})$ and since $Eh\left(Z_{i_{1}},\ldots,Z_{i_{m}}\right)<\infty$, we have

which proves the claim. We have used $\varDelta(h,w)=o(w^{-k})$. ∎

The previous theorem states sufficient conditions for a $V$-statistic or a bootstrapped $V$-statistic to converge to zero.

Let $h$ be a function of $m>1$ arguments and let $(\{Z_{t}\}_{t\in N},h\times h)$ be a of type $\varDelta$, with $\varDelta(h\times h,r)=o(r^{-4})$. If $\{G_{i}\}_{i\in N}$ is a random process, independent of $Z_{t}$, such that $\sup_{i}EG_{i}^{4}<\infty$, with notation $T_{n}=\frac{1}{n^{m-1}}\sum_{i\in N^{m}}G_{i_{1}}G_{i_{2}}h(Z_{i})$,

is uniformly bounded. We get the bound by applying Cauchy-Schwarz iteratively and using assumption $\sup_{i}EG_{i}^{4}<\infty$.

We check that the second non-central moment converges to zero,

Supremum over $n$ is needed since $EG_{i_{1}}G_{i_{2}}G_{j_{1}}G_{j_{2}}$ might change with $n$. Lemma 4, by the assumption that $(h(\cdots)\times h(\cdots),Z_{t})$ is of type $\varDelta$, the growth of the inner sum $\sum_{i,j\in N^{m}}|Eh(Z_{i})h(Z_{j})|$ is at most of order

Since $\varDelta(h\times h,r)=o(r^{-4})$, the growth rate is

For $m=2$ we have assumed existence of an extra term $o(1)$, which concludes the proof. ∎

We next prove that the asymptotic distribution of a $V$-statistic depends on number of terms in the Hoeffding decomposition that are equal to zero.

Let $h$ be a core with $m$ arguments. If $h_{0}=h_{1}=0$, and for all $c>2$ component $(h_{c}\times h_{c},Z_{t})$ is of type $\varDelta$, with $\varDelta(h_{c}\times h_{c},r)=o(r^{-4})$ then

Using Hoeffding decomposition we write the core $h$ as a sum of the components $h_{c}$ ,

$h_{0}=0$ and $h_{1}=0$. By Lemma 4, for $c\geq 3$, $nV_{n}(h_{c})$ converges to zero in mean squared. To see that it suffices to put $Q=1$ and verify that $(h_{c}\times h_{c},Z_{t})$ is of $\varDelta$ type, which is explicitly assumed. ∎

Before we study the asymptotic distribution of a bootstrapped statistic $B_{n}$ we need to sate three simple lemmas that will be frequently used.

By the definition of $W_{i}$, $E(\sum_{i=1}^{n}W_{i})^{2}\leq n2\sum_{r=1}^{n}Cov(W_{0},W_{r})=nO(l_{n})$, where $\sum_{r=1}^{n}Cov(W_{0},W_{r})=O(l_{n})$ follows from bootstrap assumption. Also, by the Bootstrap assumptions we have $\lim_{n\to\infty}\frac{l_{n}^{3}}{n^{2}}=0$. Therefore $\frac{1}{n}\sum_{i=1}^{n}W_{i}$ converges to zero in mean squared. ∎

If $\{W_{i}\}$ is a bootstrap process then

Let $f$ be a function and let $j=\{j_{1},\ldots,j_{q}\}$ be a subset of $\{1,\ldots,m\}$. Then

Each element $f(Z_{i_{j_{1}}},...,Z_{i_{j_{q}}})$ is repeated exactly $n^{m-q}$ times. ∎

We now prove an analogue of the Lemma 5 for bootstrapped statistics $B$.

and $(h_{c},Z_{t})$ for $c>2$ are of type $\varDelta$, with $\varDelta(h_{c}\times h_{c},r)=o(r^{-4})$ then

Using Hoeffding decomposition we write core $h$ as a sum of components $h_{c}$ (the ones with $h_{0},h_{1}$ are equal to zero and therefore omitted)

We will show that for almost all fixed $j_{1}<\cdots<j_{c}$ the sum 19 converges to zero.

Suppose $j_{1}>2$. The sum 19 can be written

By Lemma 4, for $c\geq 3$, $\frac{n}{l_{n}}V_{n}(h_{c})$ converges to zero in mean squared. Indeed, it is sufficient to put $G_{i}=1$ and $T_{n}=nV_{n}(h_{c})$ and notice that $\frac{n}{l_{n}}V_{n}(h_{c})=\frac{1}{l_{n}}=o(1)T_{n}$, since $l_{n}\to\infty$. Consequently, since $(\frac{1}{n}\sum_{i=1}^{n}Q_{i})^{2}$ converges to zero in mean square 6, the product, converges to zero in mean square i.e.

Suppose $j_{1}=2$. The sum 19 can be written

The latter expression $\frac{l_{n}}{n}\sum_{i=1}^{n}Q_{i}$ converges to zero in mean square. The former expression can be further decomposed

We use Lemma 4 for $\frac{1}{l_{n}}T_{+}$ and $\frac{1}{l_{n}}T_{-}$, to show that they converge to zero. We need to check that

Suppose $j_{1}=1$ and $j_{2}>2$. This case is identical to the previous case, up to swapping $i_{1},i_{2}$ in the equation 20.

Finally, suppose $j_{1}=1$ and $j_{2}=2$ and $c>2$. The sum 19 can be written

We again use Lemma 4 to see that this sum converges to zero in mean squared (we checked the assumptions above). We have proved that

So far we avoided expressing results in terms of $\tau$-mixing and degeneracy of a core, now we relate $\varDelta$ formalism to those concepts. We start with a technical lemma.

If $h$ is a Lipschitz continuous core then its components are also Lipschitz continuous.

The auxiliary function used in the Hoeffding decomposition

is Lipschitz, since $h$ is Lipschitz continuous.

$h_{0}$ is obviously Lipschitz continuous. If $h_{k}$ for $k<c$ are Lipschitz continuous then, since $g_{c}$ is Lipschitz continuous, $h_{c}$ is also Lipschitz continuous as a sum of Lipschitz continuous functions. ∎

Let $\left\{Z_{t}\right\}$ be a $\tau$-dependent stationary process and $h$ be a Lipschitz core of $m$ arguments, If for all $c>0$ $(h_{c}\times h_{c},Z_{t})$ and $(h,Z_{t})$ are of type $\varDelta$ with the rate $O(\tau(d))$ then

Let $f=h_{c}\times h_{c}$ or $f=h$. $f$ is canonical and Lipschitz continuous (if $f=h_{c}\times h_{c}$ it follows from Lemma 10). Suppose $i_{r}$ is the isolating index. Further suppose there are $k$ indexes $a_{1},\cdots,a_{k}$ smaller than $i_{r}$ and $m-k-1$ indexes greater than $i_{r}$, namely $a_{k+2},\cdots,a_{m}$. In this notation $a_{k+1}=i_{r}$.

Let us partition the vector $\left(Z_{i_{1}},\ldots,Z_{i_{m}}\right)$ into three parts:

where $a_{k+1}$ is the isolating index. If $k=0$, $A$ is empty and if $k=m-1$, $C$ is empty but this does not change our arguments below. Using Lemma [9, Lemma 5.3], we will construct $B^{*}$ and $C^{**}$ that are independent of $A$ and independent of each other and

where $w$ is an isolating distance [9, Lemma 5.3] assumes that there exists a random variable $\delta$ independent of the vector $(A,B,C)$. This assumption is important only if CDF of the vector is not continuous, we can assume that our space is endowed with such $\delta$.. Let $D=(B,C)$ The [9, Lemma 5.3] guarantees that there exist $D^{*}$ independent of $A$, such that

where $d$ is the $L_{1}$ distance on Euclidean space (non-negativity justifies dropping absolute value). By definition of $\tau$-mixing, $\tau(w)\geq\tau(\sigma(A),D)$. Since $D^{*}=(B^{*},C^{*})$ has the same distribution as $D$ (in particular it has the same $\tau$ dependence structure) we use the lemma again to construct $C^{**}$, independent of $A$ and $B^{*}$, such that

By the triangle inequality we obtain equation 21.

Since $B^{*}$ is a singleton, independent of both $A$ and $C^{**}$, by degeneracy of $f$

Note that $f(A,B^{*},C^{**})$ is just a shorthand, random variables $A,B^{*},C^{**}$ are inserted in the right order. Thus, we have that

In the proof we are going to use [Theorems 2.1, 3.1], which characterise asymptotic properties of $nV_{n}(h_{2})$ and $nB(h_{2})$. Both theorems use similar set of assumptions which we verify upfront. Assumption A2.

(i) $h_{2}$ is one-degenerate and symmetric - this follows from the Hoeffding decomposition;

(ii) $h_{2}$ is a kernel - is one of the assumptions of this theorem;

(iii) $Eh_{2}(Z_{1},Z_{1})<\infty$ – follows from $\sup_{i\in N^{6}}|Eh(Z_{i})|<\infty$ ;

(iv) $h_{2}$ is Lipschitz continuous - follows from the Lemma 10.

Assumption B1, A1. Assumption $B1$, $\sum_{r=1}^{n}r^{2}\sqrt{\tau(r)}<\infty$, is the same as ours, assumption $A1$, $\sum_{r=1}^{n}\sqrt{\tau(r)}<\infty$ is implied. Assumption B2. This assumption about the bootstrap process $W_{t}$ is the same as our Bootstrap assumptions.

Denote by $V$ the weak limit of $nV_{n}(h_{2})$, which exits by the [Theorem 2.1], and let $\mathcal{F}=\sigma(Z_{1},\cdots,Z_{n})$. By [19, Theorem 3.1], since the distribution of $V$ is continuous, we have

in probability. We show that $nB_{n}(h_{2})$ converges to $V$ weakly, by showing pointwise convergence of CDF

To change the order of limit and expectation we have dominated convergence Theorem, justified since $P(nB_{n}(h)<x|\mathcal{F})$ are bounded by 1. The difference $n(B_{n}(h)-V_{n}(h))$ is

By Lemma 9 and Lemma 5 respectively, both

converge to zero in mean square. We check assumptions: since $Z_{t}$ is tau mixing and $h$ is Lipschitz continuous, by Lemma 11 all self products of components and $Z_{t}$, $(h_{c}\times h_{c},Z_{t})$ for $c>0$, are $\varDelta$ type of order $\tau(r)$, of order at least $o(r^{-4})$ (since $\sum_{r=1}^{n}r^{2}\sqrt{\tau(r)}<\infty$). Since $h$ is one degenerate, first and zero component $h_{0},h_{1}$ are equal to zero (and so are $B(h_{0}),B(h_{1})$).

This shows that $nB_{n}(h_{2})$ converges weakly to $V$. ∎

C.2 Proof of Theorem 2

Using Hoeffding decomposition we write the core $h$ as a sum of the components $h_{c}$ ,

By the Lemma 4, for $c\geq 1$, $V_{n}(h_{c})$ converges to zero in probability. The sum associated with $h_{1}$ is

By Lemma 11 $(h_{1}\times h_{1},Z_{t})$ is $\varDelta$ type of order $o(r^{-4})$. Using Lemma 4 we get the growth rate of $E(V_{n}(h_{1}))^{2}=O(\frac{1}{n})$, thus $V_{n}(h_{1})$ converges in mean square to zero. ∎

C.3 Proof of Theorem 3

We show that the second non central moment of $B_{1}$ converges to . The second non central moment is

The inequality in the third line follows from the fact that correlations of the bootstrap process $W_{i}$ are positive (Bootstrap assumption) and

and therefore $EC\left(\frac{1}{n}\sum_{i=1}^{n}W_{i}\right)^{4}\to 0$.

We now prove that $o(n)B_{2}(h)$ converges to zero. Using Hoeffding decomposition we write core $h$ as a sum of components $h_{c}$ and $h_{0}$