A Kernel Independence Test for Random Processes

Introduction

Measures of statistical dependence between pairs of random variables $(X,Y)$ are well established, and have been applied in a wide variety of areas, including fitting causal networks (Pearl, 2000), discovering features which have significant dependence on a label set (Song et al., 2012), and independent component analysis (Hyvärinen et al., 2004). Where pairs of observations are independent and identically distributed, a number of non-parametric tests of independence have been developed (Feuerverger, 1993; Gretton et al., 2007; Székely et al., 2009; Gretton & Györfi, 2010), which determine whether the dependence measure value is statistically significant. These non-parametric tests are consistent against any fixed alternative - they make no assumptions as to the nature of the dependence.

For many data analysis tasks, however, the observations being tested are drawn from a time series: each observation is dependent on its past values. Examples include audio signals, financial data, and brain activity. Given two such random processes, we propose a hypothesis test of instantaneous dependence, of whether the two signals are dependent at a particular time $t$. Our test satisfies two important properties: it is consistent against any fixed alternatives, and it is non-parametric - we do not assume the dependence takes a particular form (such as linear correlation), nor do we require parametric models of the time series. We further avoid making use of a density estimate as an intermediate step, so as to avoid the assumption that the distributions have densities (for instance, when dealing with text or other structured data).

We use as our test statistic the Hilbert-Schmidt Independence Criterion (HSIC) (Gretton et al., 2005, 2007), which can be interpreted as the distance between embeddings of the joint distribution and the product of the marginals in a reproducing kernel Hilbert space (RKHS) (Gretton et al., 2012, Section 7). When characteristic RKHSs are used, the HSIC is zero iff the variables are independent (Sriperumbudur et al., 2010). Under the null hypothesis of independence, $P_{XY}=P_{X}P_{Y}$, the minimum variance estimate of HSIC is a degenerate U-statistic. The distribution of the empirical HSIC under the null is an infinite sum of independent $\chi^{2}$ variables (Gretton et al., 2007), which follows directly from e.g. (Serfling, 2002, Ch. 5). In practice, given a sample $(x_{i},y_{i})_{i=1}^{n}$ of pairs of variables drawn from $P_{XY}$, the null distribution is approximated by a bootstrap procedure, where a histogram is obtained by computing the test statistic on many different permutations $\{x_{i},y_{\pi(i)}\}_{i=1}^{n}$, to decouple $X$ and $Y$.

In the case where the samples $Z_{t}=(X_{t},Y_{t})$ are drawn from a random process, the analysis of the asymptotic behaviour of HSIC requires substantially more effort than in the i.i.d. case. As our main contribution, we obtain both the null and alternative distributions of HSIC for random processes, where the null distribution is defined as $X_{t}$ being independent of $Y_{t}$ at time $t$. Such a test may be used for rejecting causal effects (i.e., whether one signal is not dependent on the values of another signal at a particular delay) or instant coupling (see our first experiment in Section 4.2).We distinguish our case from the problem of ensuring time series are independent simultaneously across all time lags, e.g the null will hold even if $X_{t}=Y_{t-1}$ where $Y_{t}$ is white noise. The null distribution is again an infinite weighted sum of $\chi^{2}$ variables, however these are now correlated, rather than independent. Under the alternative hypothesis, the statistic has an asymptotically normal distribution.

For the test to be used in practice, we require an empirical estimate of the null distribution, which gives the correct test threshold when $Z_{t}=(X_{t},Y_{t})$ is a random process. Evidently, the bootstrap procedure used in the i.i.d. case is incorrect, as the temporal dependence structure within the $Y_{t}$ will be removed. This turns out to cause severe problems in practice, since the permutation procedure will give an increasing rate of false positives as the temporal dependence of the $Y_{t}$ increases (i.e., dependence will be detected between $X_{t}$ and $Y_{t}$, even though none exists, this is also known as a Type I error). Instead, our null estimate is obtained by making shifts of one signal relative to the other, so as to retain the dependence structure within each signal. Consequently, we are able to keep the Type I error at the designed level $\alpha=0.05$. In our experiments, we address three examples: one artificial case consisting of two signals which are dependent but have no correlation, and two real-world examples on forex data. HSIC for random processes reveals dependencies that classical approaches fail to detect. Moreover, our new approach gives the correct Type I error rate, whereas a bootstrap-based approach designed for i.i.d. signals returns too many false positives.

Prior work on testing independence in time series may be categorized in two branches: testing serial dependence within a single time series, and testing dependence between one time series and another. The case of serial dependence turns out to be relatively straightforward, as under the null hypothesis, the samples become independent: thus, the analysis reduces to the i.i.d. case. Pinkse (1998); Diks & Panchenko (2005) provide a quadratic forms function-based serial dependence test which employs the same statistic as HSIC. Due to the simple form of the null hypothesis, the analysis of (Serfling, 2002, Ch. 5) applies. Further work in the context of the serial dependency testing includes simple approaches based on rank statistics e.g. Spearman’s correlation or Kendall’s tau, correlation integrals e.g. (Broock et al., 1996); criteria based on integrated squared distance between densities e.g (Rosenblatt & Wahlen, 1992); KL-divergence based criteria e.g. (Robinson, 1991; Hong & White, 2005); and generalizations of KL-divergence to so called $q$-class entropies e.g. (Granger et al., 2004; Racine & Maasoumi, 2007).

In most of the tests of independence of two time series, specific conditions have been enforced, e.g that processes follow a moving average specification or the dependence is linear. Prior work in the context of dependency tests of two time series includes cross covariance based tests e.g. (Haugh, 1976; Hong, 1996; Shao, 2009); and a Generalized Association Measure based criterion (Fadlallah et al., 2012). Some work has been undertaken in the non-parametric case, however. A non-parametric measure of independence for time series, based on the Hilbert-Schmidt Independence Criterion, was proposed by Zhang et al. (2008). While this work established the convergence in probability of the statistic to its population value, no asymptotic distributions were obtained, and the statistic was not used in hypothesis testing. To our knowledge, the only non-parametric independence test for pairs of time series is due to Besserve et al. (2013), which addresses the harder problem of testing independence across all time lags simultaneously. Let $X_{t}$ follow a MA(2) model and put $Y_{t}=X_{t-20}$. This is a case addressed by Besserve et al. (2013), who will reject their null hypothesis, whereas our null is accepted The procedure is to compute the Hilbert-Schmidt norm of a cross-spectral density operator (the Fourier transform of the covariance operator at each time lag). The resulting statistic is a function of frequency, and must be zero at all frequencies for independence, so a correction for multiple hypothesis testing is required. It is not clear how the asymptotic analysis used in the present work would apply to this statistic, and this remains an interesting topic of future study.

The remaining material is organized as follows. In Section 2 we provide a brief introduction to random processes and various mixing conditions, and an expression for our independence statistic, HSIC. In Section 3, we characterize the asymptotic behaviour of HSIC for random variables with temporal dependence, under the null and alternative hypotheses, and establish the test consistency. We propose an empirical procedure for constructing a statistical test, and demonstrate that the earlier bootstrap approach will not work for our case. Section 4 provides experiments on synthetic and real data.

Background

In this section we introduce necessary definitions referring to random processes. We then go on to define a V-statistic estimate of the Hilbert-Schmidt Independence Criterion, which applies in the i.i.d. case.

Next, we formalize a concept of memory of a process. A process is called absolutely regular ($\beta$-mixing) if $\beta(m)\rightarrow 0$, where

The second supremum in the $\beta(m)$ definition is taken over all pairs of finite partitions $\{A_{1},\cdots,A_{I}\}$ and $\{B_{1},\cdots,B_{J}\}$ of the sample space such that $A_{i}\in\mathcal{A}_{1}^{n}$ and $B_{j}\in\mathcal{A}_{n+m}^{\infty}$, and $\mathcal{A}_{b}^{c}$ is a sigma field spanned by a subsequence, $\mathcal{A}_{b}^{c}=\sigma(Z_{b},Z_{b+1},...,Z_{c})$. A process is called uniform mixing ($\phi$-mixing) if $\phi(m)\rightarrow 0$, where

Uniform mixing implies absolute regularity, i.e. $\beta(m)\leq\phi(m)$ (Bradley et al., 2005). Under technical assumptions, Autoregressive Moving Average processes — or more generally Markov Chains — are absolutely regular or uniformly mixing (Doukhan, 1994).

Hilbert-Schmidt Independence Criterion

Let $k$, $l$ be positive definite kernels associated with respective reproducing kernel Hilbert spaces $\mathcal{H}_{\mathcal{X}}$ on $\mathcal{X}$, and $\mathcal{H}_{\mathcal{Y}}$ on $\mathcal{Y}$. We assume that $k$ and $l$ are bounded and continuous. We associate to the random variable $X$ a mean embedding $\mu_{X}(x):=\mathcal{E}_{X}k(X,x)$, such that $\forall f\in\mathcal{H}_{\mathcal{X}}$, $\langle f,\mu_{X}\rangle_{\mathcal{H}_{\mathcal{X}}}=\mathcal{E}_{X}(f(X))$ (Berlinet & Thomas-Agnan, 2004; Smola et al., 2007). We assume $k$, $l$ are characteristic kernels, meaning the mappings $\mu_{X}$ and $\mu_{Y}(y):=\mathcal{E}_{Y}l(Y,y)$ are injective embeddings of the probability measures to the corresponding RKHSs; i.e., distributions have unique embeddings (Fukumizu et al., 2008; Sriperumbudur et al., 2010).

We next recall a measure of statistical dependence, the Hilbert-Schmidt Independence Criterion (HSIC), which can be expressed in terms of expectations of RKHS kernels (Gretton et al., 2005, 2007). Denote a group of permutations over 4 elements by $S_{4}$, with $\pi$ one of its elements, i.e., a permutation of four elements. We define

Let $\gamma$ be an expected value of the function $h$, $\gamma=\mathcal{E}h(Z_{1}^{*},Z_{2}^{*},Z_{3}^{*},Z_{4}^{*})$. This expectation corresponds to HSIC, computed using a function symmetric in its arguments. For $k$ and $l$ characteristic, continuous, translation invariant, and vanishing at infinity, $\gamma$ is equal to zero if and only if the null hypothesis holds (see (Lyons, 2013, Lemma 3.8), applying (Sriperumbudur et al., 2011, Proposition 2), and the note at the end of Section 5).

The value of $\gamma$ corresponds to a distance between embeddings of $(X_{1}^{*},Y_{2}^{*})$ and $(X_{1}^{*},Y_{1}^{*})$ to an RKHS with the product kernel $\kappa=k\cdot l$ (Gretton et al., 2012, Section 7). A biased empirical estimate of the Hilbert-Schmidt Independence Criterion can be expressed as a $V$-statistic (the unbiased estimate is a U-statistic, however the difference will be accounted for when constructing a hypothesis test, through an appropriate null distribution).

V𝑉V statistics.

A $V$-statistic of a $k$-argument, symmetric function $f$ is written

If $j=k-1$ we say that the function is canonical. We refer to a normalized $V$ statistic as a $V$-statistic multiplied by the sample size, $n\cdot V$.

HSIC for random processes

In this section we construct the Hilbert-Schmidt Independence Criterion for random processes, and define its asymptotic behaviour. We then introduce an independence testing procedure for time series.

We introduce two hypotheses: the null hypothesis $\mathbf{H_{0}}$ that $X_{t}$ and $Y_{t}$ are independent, and the alternative hypothesis $\mathbf{H_{1}}$ that they are dependent. To build a statistical test based on $n\cdot V(h,Z)$ we need two main results. First, if null hypothesis holds, we show $n\cdot V(h,Z)$ converges to a random variable. Second, if the null hypothesis does not hold, the $n\cdot V(h,Z)$ estimator diverges to infinity. Following these results, the Type I error (the probability of mistakenly rejecting the null hypothesis) will stabilize at the design parameter $\alpha$, and the Type II error (the probability of mistakenly accepting the null hypothesis when the variables are dependent) will drop to zero, as the sample size increases.

We begin by introducing an auxiliary kernel function $s$, and characterize the normalized $V$-statistic distribution of $s$ using a CLT introduced by (Borisov & Volodko, 2008). We then show that the normalized $V$-statistic associated with the function $s$ has the same asymptotic distribution as the $n\cdot V(h,Z)$ distribution.

By Steinwart & Scovel (2012) Corollary 3.5, the bounded, continuous kernel $s$ has a representationA bounded kernel is compactly embedded into $L^{2}(\mathbf{Z},\mathcal{B}(\mathbf{Z}),P_{\mathbf{Z}})$ (Steinwart & Scovel, 2012).

Let the process $Z_{t}$ have a mixing coefficient smaller than $m^{-3}$ $(\beta(m),\phi(m)\leq m^{-3})$ and satisfy either of the following conditions:

$Z_{t}$ is $\beta$-mixing. For some $\epsilon>0$ and for an even number $c\geq 2$, the following holds

$\sup_{i}\mathcal{E}|e_{i}(X_{1})|^{2+\epsilon}\leq\infty$, where $e_{i}$ is the basis introduced in the Statement 1 and $|\cdot|$ denotes an absolute value.

$\sum_{m=1}^{\infty}\beta^{\epsilon/(2+\epsilon)}(m)<\infty$.

If the null hypothesis holds, then $s$ is a canonical function and a kernel. What is more,

where $\tau_{j}$ is a centred Gaussian sequence with the covariance matrix

We now characterize the asymptotics of $V(h,Z)$.

Under assumptions of Lemma 2, if $\mathbf{H_{0}}$ holds, then the asymptotic distribution of the empirical HSIC (with scaling $n$) is the same as that of $n\cdot V(s,Z)$,

Under assumptions of the Lemma 2, if $\mathbf{H_{1}}$ holds, then $\gamma>0$ and $\sqrt{n}(V(h,Z)-\gamma)$ has asymptotically normal distribution with mean zero and finite variance.

Consequently, if the null hypothesis does not hold then $P(n\cdot V(h,Z)>C)=P(V(h,Z)>\frac{C}{n})\to 1$ for any fixed $C$. Finally, we show that the $\gamma$ estimator is easy to compute. According to Gretton et al. (2007, equation 4), $V(h,Z)=n^{-2}trHKHL,$ where $K_{ab}=k(X_{a},X_{b})$, $L_{ab}=l(Y_{a},Y_{b})$ ,$H_{ij}=\delta_{ij}-n^{-1}$ and $n$ is a sample size.

We begin by showing that the $H_{0}$ distribution of the $\gamma$ estimator obtained via the bootstrap approach of (Diks & Panchenko, 2005; Gretton et al., 2007) gives an incorrect p-value estimate when used with independent random processes. In fact, the null hypothesis obtained by permutation corresponds to the processes being both i.i.d. and independent from each other. Recall the covariance structure of the $\gamma$ estimator from Theorem 1,

We can represent $e_{a}$ and $e_{b}$ as $e_{a}(z)=e^{X}_{u}(x)e^{Y}_{o}(y)$, $e_{b}(z)=e^{X}_{i}(x)e^{Y}_{p}(y)$, as a decomposition of the $\mathbf{Z}$ basis into bases of $\mathbf{X}$ and $\mathbf{Y}$, respectively. Consider a partial sum $T_{n}$ of series from the above equation (3), with $X_{t}$ replaced with its permutation $X_{\pi(t)}$,

Using covariance inequalities from (Doukhan, 1994, Section 1.2.2) we conclude that $\mathcal{E}e^{Y}_{o}(Y_{1})e^{Y}_{p}(Y_{j+1})=O(\Lambda(j)^{\frac{1}{2}})$ and $\mathcal{E}e^{X}_{u}(X_{\pi(1)})e^{X}_{i}(X_{\pi({j+1})})=O(\Lambda(|\pi(j)-\pi(1)|)^{\frac{1}{2}})$ where $\Lambda$ is an appropriate mixing coefficient ($\beta$ or $\phi$). Recall that $0<\Lambda(j)<Cj^{-3}$.

We can therefore reduce the problem to the convergence of a random variable

where $\pi$ is a random permutation drawn from the uniform distribution over the set of $n$-element permutations. In the supplementary material we show that this sum converges in probability to zero.

Since $S_{n}>T_{n}>0$, then $T_{n}$ converges to zero in probability, and consequently the covariance matrix entry $\mathcal{E}\tau_{a}\tau_{b}$ converges to unity for $a=b$, and to zero otherwise. Indeed, the expected value $\mathcal{E}e_{a}((X_{\pi(1)},Y_{1}))e_{b}((X_{\pi(1)},Y_{1}))=0$ if $a\neq b$ and is equal to one otherwise. Note that this is the covariance matrix described by Gretton et al. (2007).

A correct approach to approximating the asymptotic null distribution of $n\cdot V(h,Z)$ under $\mathbf{H}_{0}$ is by shifting of one time series relative to the other. Define the shifted process $S^{c}_{t}=Y_{t+c\text{ mod }n}$ for an integer $c$, $0\leq c\leq n$ and $0\leq t\leq n$. If we let $c$ vary over $0\leq A\leq B\leq n$ for $A$ such that the dependence between $Y_{t+A}$ and $X_{t}$ is negligible, then we can approximate the null distribution with an empirical distribution calculated on points $(V(h,Z^{k}))_{A\leq k\leq B}$, where $Z^{k}_{t}=(X_{t},S^{k}_{t})$. This is due to the fact that the shifted process $S^{c}_{t}$ retains most of the dependence, since it does not scramble the time index.As a illustration, consider $W_{t}=Y_{t-10}$. If $Y_{t}$ is stationary then the dependence structure of $(W_{t_{1}},W_{t_{2}})$ and $(Y_{t_{1}},Y_{t_{2}})$ is the same. If we set $W_{t}=$ $Y_{\pi(t)}$ this property does not hold. We call this method Shift HSIC. In the supplementary material we show that Shift HSIC samples from the correct null distribution.

Experiments

In the experiments we compare Shift HSIC with the Bootstrap HSIC of Gretton et al. (2007). We investigate three cases: an artificial dataset, where two time series are coupled non-linearly; and two forex datasets, where in one case we seek residual dependence after one time series has been used to linearly predict another, and in the other case, we reveal strong dependencies between signals that are not seen via linear correlation.

We investigate two dependent, autoregressive random processes $X_{t}$,$Y_{t}$, specified by

with an autoregressive component $a$. The coupling of the processes is a result of the dependence in the innovations $\epsilon_{t},\eta_{t}$. These $\epsilon_{t},\eta_{t}$ are drawn from an Extinct Gaussian distribution, defined in Algorithm 1. The parameter $p$ (called extinction rate) controls how often a point drawn form a ball $B(0,r)$ dies off. According to Algorithm 1, the probability of seeing a point inside the ball $B(0,r)$ is different than for a two dimensional Gaussian $N(\mathbf{0},Id)$. On the other hand, as $p$ goes to zero, the Extinct Gaussian converges in distribution to $N(\mathbf{0},Id)$. Figure 1 illustrates the joint distribution of $X_{t},Y_{t}$. The left scatter plot in Figure 1 presents $X_{t}$ and $Y_{t}$ generated with an extinction rate of $50\%$, while the right hand plot is generated with an extinction rate of $99.87\%$. Processes used in this experiment had an autoregressive component of $0.2$, and the radius of the innovation process was $1$.

Figure 2 compares the power of the Shift HSIC test and the correlation test. The $X$ axis represents an extinction rate, while the $Y$ axis shows the true positive rate. Shift HSIC is capable of detecting non-linear dependence between $X_{t}$ and $Y_{t}$, which is missed by linear correlation. The red star depicts performance of the KCSD algorithm developed by Besserve et al. (2013), with parameters tuned by its authors: note that this result required using four times as many data points as HSIC.

False positive rates.

We next investigate the rate of false positives for Shift HSIC and Bootstrap HSIC on independent copies of the $AR(1)$ processes used in the previous experiment. To generate independent processes, we first sampled two pairs $(X_{t},Y_{t})$, $(X_{t}^{\prime},Y_{t}^{\prime})$ of time series using (6), and then constructed $Z$ by taking $X$ from the first pair and $Y$ from the second, i.e., $Z_{t}=(X_{t},Y_{t}^{\prime})$. We set an extinction rate to $50\%$. As a reviewer pointed out, the example for the FP rates can be simplified, however we decided to be consistent with the marginal distribution of $X_{t}$,$Y_{t}$ across the experiments. The AR component $a$ in the model (6) controls the memory of a processes - the larger this component, the longer the memory. We performed the Shift HSIC and the Bootstrap HSIC tests on $Z_{t}$ generated under $\mathbf{H_{0}}$ with different AR components. Figure 3 illustrates the results of this experiment. The $X$ axis is indexed by the AR component and $Y$ axis shows the FP rate. As the temporal dependence increases, the Bootstrap HSIC incorrectly gives an increasing number of false positives: thus, it cannot be relied on to detect dependence in time series. The Shift HSIC false positive rate remains at the targeted $5\%$ p-value level.

2 Forex data

We use Foreign Exchange Market quotes to evaluate Shift HSIC performance on the real life data. Practitioners point out that forex time series are noisy and hard to handle, especially at low granulations (smaller then 15 minutes). We decided to work with forex time series to show that Shift HSIC can detect dependence even on such a difficult dataset. The forex time series were granulated to obtain two minute sampling (the granulation function returned the last price in the two minute window). Using the test of Diks & Panchenko (2005), we checked that serial dependence of the differentiated time series decays fast enough to satisfy the assumed mixing conditions (by a differentiated time series, we refer to $(X_{t}-X_{t-1})_{t\in\mathbf{N}}$). The choice of the pairs and trading day (21st January 2013) were arbitrary.

Having one Australian dollar we may obtain a quantity of Yen in two ways, either by using AUD/JPY exchange rate explicitly or by buying Canadian dollars and then selling them at the CAD/JPY rate. Let $X_{t}$ be a differentiated AUD/JPY exchange rate and $Y_{t}$ be a differentiated product of exchange rates AUD/CAD$\times$CAD/JPY. We will investigate the relation between these two. Common sense dictates that $Y_{t}$ should behave similarly to $X_{t}$. After examining the cross-correlation of $X_{t}$ and $Y_{t}$, we propose a simple regression model to describe the interaction between the signals,

We fit the model and see that $a_{0}=0.97$, and the remaining coefficients are not bigger then $0.06$ in absolute value. This suggest that most of the dependence is explained by an instantaneous coupling. We further investigate the cross-correlation between residuals $R_{t}=Y_{t}-\hat{Y}_{t}$ and $X_{t}$. We observe no significant correlations in the first 30 lags.

Next we investigate dependence of residuals with lagged values of the explanatory variables, i.e., $R_{t}$ with $X_{t-k}$ for $k\in(0,\cdots,30)$. After calculating p-values using the Bootstrap HSIC and the Shift HSIC, we discover dependence only at lags $4$, $5$, $9$, $13$ and $29$, as presented in the Figure 4. Lack of the dependence at lag zero suggests that the linear model for coupling is reasonable. However, both the Bootstrap HSIC and the Shift HSIC support the hypothesis that there is a strong relation at lag $5$, which is not explained well by the linear model.

The questions remains whether test statistics at lags $4$, $9$, $13$ and $29$ indicate further model misspecification. Under $\mathbf{H_{0}}$, at a significance level $94\%$, we expect 1.8 out of 30 statistics to be higher than the 94th percentile. Excluding the statistic at lag $5$, the Shift HSIC test reports two statistics above this percentile, while Bootstrap HSIC reports four. Should the statistics at the different lags be independent from each other, the probabilities of seeing two and four statistics above the percentile are respectively $25\%$ and $6\%$. Shift HSIC indicates that the model fits the data well, while the Bootstrap HSIC suggests that some non-linear dependencies remain unexplained.

Dependence structure.

The data are five currency pairs. A correlation based independence test, and the Shift HSIC test, were performed on each pair of currencies. The dependencies revealed by these tests are depicted in Figure 5 - nodes represent the time series and edges represent dependence. Shift HSIC reveals a strong coupling between EUR/RUB and USD/JPY, HKD/JPY and XAU/USD that was not found by simple correlation. All edges revealed by Shift HSIC have p-values at most at level $0.03$ - clearly, the Shift HSIC managed to find a strong non-linear dependence. Note that the obtained graphs are cliques.

Proofs

A $U$-statistic of a $k$-argument, symmetric function $f$, is written

A decomposition due to Hoeffding allows us to decompose this U-statistic into a sum of $U$-statistics of canonical functions, $U(h,Z)=\sum_{k=1}^{l}{l\choose k}U(h_{k},Z)$, where $h_{k}(z_{1},...,z_{l})$ are components of the decomposition. According to Serfling (2002, section 5.1.5), each of $h_{1}$,$h_{2}$,$h_{3}$,$h_{4}$ is symmetric and canonical. Note that $h_{k}$ is defined using independent samples $Z^{*}$ - this is because the CLT or LLN state that U-statistics or V-statistics of mixing processes converge to their expected value taken with respect to independent copies, i.e., $Z^{*}$. Under $\mathbf{H}_{0}$, $h_{1}$ is equal to zero everywhere and $h_{2}=\frac{1}{6}s$, where these results were obtained by Gretton et al. (2007).The second result is hard to locate - it is in appendix A.2, text between equations 12 and 13 See supplementary material for details.

In order to characterize $U(h,Z)$, we show that under null hypothesis $U(h_{2},Z)$ converges to a random variable, and both $U(h_{3},Z)$,$U(h_{4},Z)$ converge to zero in a probability (the latter proof can be found in the supplementary material). Bellow we characterise $U(h_{2},Z)$ convergence.

First recall that under null hypothesis $h_{2}=\frac{1}{6}s$. We will check the conditions of (Borisov & Volodko, 2008, Theorem 1) (also available in the supplementary). First, from Mercer’s Theorem (Steinwart & Scovel, 2012, Corollary 3.5), we deduce that the $h_{2}$ coefficients in $L_{2}(\mathbf{Z},\mathcal{B}_{\mathbf{Z}},P_{\mathbf{Z}})$ are absolutely summable. In the supplementary material we show that $\mathcal{E}e_{i}(Z_{1}^{*})=0$. Recall the assumptions of Lemma 2. If A holds then $\sum_{k=1}^{\infty}\phi(k)^{\frac{1}{2}}<\infty$ and $\sup_{i}\mathcal{E}|e_{i}(X_{1})|^{2}=1<\infty$ ($e_{i}$ is an orthonormal eigenfunction). Finally, if B holds then the process $Z_{t}$ is $\alpha$-mixing. The remaining assumptions concerning uniform mixing in Borisov & Volodko (2008) are exactly the same as in this lemma. ∎

(Lemma 2) We use the fact that $h_{2}$ is equal to $s$ up to scaling ($6U(h_{2},Z)=U(s,Z)$), and Lemma 3, to see that $nU(s,Z)\stackrel{{\scriptstyle D}}{{\to}}\sum_{i}^{\infty}\lambda_{i}(\tau_{i}^{2}-1)$. Since $\mathcal{E}s(Z_{t},Z_{t})=\mathcal{E}\sum_{i=1}^{\infty}\lambda_{i}e_{i}(Z_{t})^{2}=\sum_{i=1}^{\infty}\lambda_{i}$, then by the LLN for mixing processes,

We use a relationship between $U$ and $V$ statistics,

(Theorem 1) We operate under the null hypothesis. Recall that $U(h,Z)$ can be decomposed as $U(h,Z)=\sum_{k=1}^{4}{4\choose k}U(h_{k},Z)$. Here $h_{1}\equiv 0$. We show in the supplementary material that $U(h_{3},Z)$ and $U(h_{4},Z)$ tend to zero in probability. From Lemma 3,

We define an auxiliary symmetric function $w$,

It is obvious that $\mathcal{E}w(Z_{1}^{*},Z_{2}^{*},Z_{3}^{*})$ $=$ $6\mathcal{E}h(Z_{1}^{*},Z_{1}^{*},Z_{2}^{*},Z_{3}^{*})$. We consider the difference between the unnormalized $V$ and $U$ statistics,

where $\sum_{i\in C_{m}}$ denotes summation over all $n\choose m$ combinations of $m$ distinct elements $\{i_{1},\cdots,i_{m}\}$ from $\{1,\cdots,n\}$. The difference is equal to the sum over $4$-tuples with at least one pair of equal elements. We can choose such tuples in $\binom{4}{2}=6$ ways. Observe that $w$ covers the choice of all these six tuples. Since for any $z_{1},z_{2}\in\mathbf{Z}$, $h(z_{1},z_{1},z_{1},z_{2})=0$, then $w$ is zero whenever more than two indices are equal. Therefore we can sum $w$ over distinct indices $z_{1},z_{2},z_{3}$,

We see that $S_{n}$ is almost a $U$-statistic ($U(w,Z)$). By the CLT for $U$-statistics from Denker & Keller (1983), Theorem 1(c), we obtain

On the other hand, via the relation $h_{2}=\frac{1}{6}s$ and the $h_{2}$ definition, we get $\mathcal{E}s(Z_{1}^{*},Z_{1}^{*})=6\mathcal{E}h(Z_{1}^{*},Z_{1}^{*},Z_{2}^{*},Z_{3}^{*})$, and therefore

We normalize by $\frac{1}{n(n-1)(n-2)}$, and take the limit in $n$,

We substitute (9) and (8) on the right hand side, and use equation (7) from Lemma 2 to replace $\lim_{n\to\infty}\frac{1}{n}\sum_{i}^{n}s(Z_{i},Z_{i})$ with $\sum_{i=1}^{\infty}\lambda_{i}$, yielding

(Theorem 2) If the null hypothesis does not hold, then $\gamma>0$ (Gretton et al., 2005). In this case $h$ is nondegenerate, and we can use Denker & Keller (1983, Theorem 1(c)) to see that $\frac{\sqrt{n}}{4\ \sqrt{\sigma}}(V(h,Z)-\gamma)\sim N(0,1)$, where $\sigma$ is finite (see the note below Theorem 1 of (Denker & Keller, 1983), stating that in case (c) $\sigma^{2}$ is finite, and the note above Theorem 1 stating that $\sigma^{2}=\lim_{n\to\infty}n^{-1}\sigma_{n}^{2}$ ). ∎

(Lemma 1) We use Lemma 1 and Theorem 4 from Gretton et al. (2005) to show that $\mathcal{E}h(Z_{1}^{*},Z_{2}^{*},Z_{3}^{*},Z_{4}^{*})=0$ iff $(X_{1}^{*},Y_{1}^{*})$ has a product distribution. Since $Z_{1}^{*}\stackrel{{\scriptstyle D}}{{=}}Z_{1}$ and $Z_{t}\stackrel{{\scriptstyle D}}{{=}}Z_{1}$, we infer that $X_{t}$ is independent from $Y_{t}$ iff $\mathcal{E}h(Z_{1}^{*},Z_{2}^{*},Z_{3}^{*},Z_{4}^{*})=0$. ∎

The authors thank the reviewers and colleagues for helpful feedback, especially M. Skomra, D. Toczydlowska, and A. Zaremba.

References

Appendix A A Kernel Independence Test for Random Processes - Supplementary

The sections in the supplementary material are in the same order those in the article. In particular, the $n$-th reference to the supplementary in the article is $n$-th subsection in the supplementary material.

The arXiv version of the report and supplementary may be found at: http://arxiv.org/abs/1402.4501

Before we start, we cite (Yoshihara, 1976, Lemma 1), which will be used below.

Note that if a function $g$ is symmetric, then we can always reorder its arguments if necessary.

Let $\pi$ be a permutation drawn from a uniform distribution over the set of $n$-element permutations. We will prove that the random variable

converges to zero in probability at rate $O(n^{-1})$. Since $0\leq S_{n}\leq Q_{n}$, then $S_{n}$ converges to zero in probability at the same rate.

$\mathcal{E}|\pi(1)-\pi(i)|^{-{\frac{3}{2}}}=O(n^{-1})$.

Let $j$ be a positive integer smaller than $n$. Observe that the sum $\sum_{i}^{n}|j-i|^{-{\frac{3}{2}}}$ is finite,

where $\zeta(\cdot)$ is the Riemann zeta function. Now expand the expected value $\mathcal{E}|\pi(1)-\pi(i)|^{-{\frac{3}{2}}}$ using a conditional expected value,

If $k\neq j$ are positive integers smaller than $n$, then

We use the inequality (10) and properties of a conditional expected value.

$Q_{n}$ converges to zero in probability. The convergence rate is $\frac{1}{n}$.

First, using Lemma 5 , we compute the expected value of $Q_{n}$

Next, using Lemma 6, we compute the second moment

Using the Chebyshev’s inequality we obtain the required result. ∎

A.2 Testing procedure - Shift HSIC samples from the right distribution

We will investigate the value of the $V$-statistic for a shifted process i.e. $nV(h,Z^{k})$.

If the null hypothesis holds, then $X_{t}$ and $Y_{t+k}$ are independent for any $k$. To see this, suppose that there exists $k$ for which $X_{t}$ and $Y_{t+k}$ are dependent (the processes are stationary, so this is true for all $t$). The observation $X_{t}$ depends on its past values: in particular, $X_{t-k}$ is a parent of $X_{t}$. If in addition $X_{t-k}\rightarrow Y_{t}$, then $Y_{t}$ and $X_{t}$ will be dependent, as they share a parent.

We will use this fact to show that the $nV(h,Z^{k})$ has the same distribution as the $nV(h,Z)$. Recall the covariance structure of $nV(h,Z)$ from Theorem 1,

We represent $e_{a}$ and $e_{b}$ as $e_{a}(z)=e^{X}_{u}(x)e^{Y}_{o}(y)$, $e_{b}(z)=e^{X}_{i}(x)e^{Y}_{p}(y)$. This represents a decomposition of the basis of $\mathbf{Z}$ into basis of $\mathbf{X},\mathbf{Y}$, respectively. Consider one of the above infinite sums with $Y_{t}$ replaced with the shifted process $S^{k}_{t}$,

We obtain the following covariance structure for $nV(h,Z^{k})$,

We have used the fact that $Y_{t}$ is stationary, $\mathcal{E}e^{Y}_{o}(Y_{1+k})e^{Y}_{p}(Y_{j+1+k})=e^{Y}_{o}(Y_{1})e^{Y}_{p}(Y_{j+1})$ and that the pairs $(X_{1},X_{1+j})$, $(Y_{1+k},Y_{j+1+k})$ are independent (because $X_{t}$ and $Y_{t+k}$ are independent for all shifts $k$). For the second term,

we have used covariance inequalities from Doukhan (1994, section 1.2.2) and our bounds on mixing coefficients to obtain that when $j\geq n-k$, then $\mathcal{E}|e^{X}_{u}(X_{1})e^{X}_{i}(X_{j+1})|\leq(n-k)^{-\frac{3}{2}}$ (and by e.g. Holders inequality, $\mathcal{E}e^{Y}_{o}(Y_{1+k})e^{Y}_{p}(Y_{1+(n-j)})$ is finite). The first component takes the form

Here $\sum_{j=1}^{n-k}\mathcal{E}e_{a}(Z_{1})e_{b}(Z_{j+1})$ converges to $\sum_{j=1}^{\infty}\mathcal{E}e_{a}(Z_{1})e_{b}(Z_{j+1})$ from equation (14). Since $\mathcal{E}e_{a}(Z^{k}_{1})e_{b}(Z^{k}_{1})=\mathcal{E}e_{a}(Z_{1})e_{b}(Z_{1})$, the covariance structure from equation (14) is recovered.

Null hypothesis does not hold.

In this case, the dependence between $X_{t}$ and $Y_{t+k}$ decreases as $k$ increases, since the mixing coefficients for each of the time series converges to zero. In the limit of large $k$ and $n$, the normalized $V$-statistic will converge to the null distribution, where $X_{t}$ and $Y_{t}$ are independent random processes. The proof of this result under the assumed mixing conditions, with suitable conditions on the increase of $k$ with $n$, is a topic of future work (the next two sections give an outline of the results that would need to be established for the shifted process).

A.3 Proofs - Hoeffding decomposition

The Hoeffding decomposition (e.g. Serfling, 2002) allows us to decompose U-statistics into a sum of simpler $U$-statistics that can be easier to analyse. In the following section we will perform a Hoeffding decomposition of $U(h,Z)$ and investigate some of its properties. In the sequel we assume that $k$ and $l$ are bounded kernels. For the $U$ statistic $U(h,Z)$, we call the function $h$ a core.

Any U-statistic can be written as a sum of V-statistics with degenerate cores. To show this, we define the auxiliary functions

for each $c=1,...,m-1$ and put $g_{m}=h$.

We assume the expected value of the core with respect to starred $\{Z_{t}\}$ is zero, i.e., $\mathcal{E}h(Z_{1}^{*},\cdots,Z_{m}^{*})=0$. The canonical functions that enable the core decomposition are

We call these functions components of a core.

The U-statistic of a core function $h$ can be written as a sum of U-statistics with degenerate cores,

Recall that $\sum_{i\in C_{m}}$ denotes summation over all $n\choose m$ combinations of $m$ distinct elements $\{i_{1},\cdots,i_{m}\}$ from $\{1,\cdots,n\}$.

Under $\mathbf{H_{0}}$, $\forall z\in\mathbf{Z}$ $h_{1}(z)=0$.

We use the shorthand notation $k(a,b)\equiv k(x_{a},x_{b})$, $l(a,b)\equiv l(y_{a},y_{b})$, such that

Let us expand this expression. By using the symmetry of $k$ and $l$, and writing the arguments in lexicographical order, we obtain

Finally we introduce colours to picture grouping of terms that will cancel each other during integration.

We will show that brown terms of equation (16) cancel each other. Recall that $h_{1}(z_{1})=\mathcal{E}h(z_{1},Z_{2}^{*},Z_{3}^{*},Z_{4}^{*})$. Without loss of generality we may assume that we integrate with respect to all variables but $x_{a}$ and $y_{a}$. Observe that

Define $q=\mathcal{E}k(x_{a},X_{b}^{*})$, $p=\mathcal{E}l(y_{a},Y_{b}^{*})$. Therefore, after integration, the brown terms of the equation can be written as

Similar reasoning shows that red, green and violet terms cancel out. ∎

A component of a core function is a canonical core.

We will use induction by components’ index to show that $h_{c}$ is degenerate. The expected value of the first component is zero, indeed $\mathcal{E}h_{1}(Z_{1}^{*})=\mathcal{E}h(Z_{1}^{*},...,Z_{m}^{*})=0$. Suppose that for all $c^{\prime}$ smaller then $c$ degeneracy holds. Using component symmetry it is enough to show that the expected value $\mathcal{E}h_{c}(z_{1},...,Z_{c}^{*})$ is equal to zero. We can write

Now the first sum $\sum_{1\leq i_{1}<...<i_{c^{\prime}}\leq c-1}h_{c^{\prime}}(z_{i_{1}},...,z_{i_{c^{\prime}}})$ does not contain term $z_{c}$ so integration with respect to $Z_{c}^{*}$ does not affect it. On the other hand, by induction assumption $\mathcal{E}\sum_{1\leq i_{1}<...<i_{c^{\prime}-1}<c}h_{c^{\prime}}(z_{i_{1}},...,Z_{c}^{*})=0$. Obviously $\mathcal{E}g_{c}(z_{1},...,Z_{c}^{*})=g_{c-1}(z_{1},...,z_{c-1})$. Using these observations we obtain

Since the set $\{1\leq i_{1}<...<i_{c-1}\leq c-1\}$ contains only one sequence,

For this nice simplification we have used definition of the component $h_{c-1}$. ∎

We use that $h_{2}$ is canonial, and the exact form of $\mathcal{E}h(z_{1},z_{2},Z_{3}^{*},Z_{4}^{*})$ from (Gretton et al., 2007), Section A.2, text between equation 12 and 13. ∎

Under $\mathbf{H_{0}}$, $h_{2}=\frac{1}{6}s$.

Let $N:=\{1,\cdots,n\}$, and let $B$ be a set of all strictly increasing $4$-tuples, $B\subset N^{4}$. A $U$-statistic can be expressed as sum over elements of $B$,

If the variance of this random variable goes to zero,

then using Chebyshev’s inequality we can conclude that it converges to a constant in probability. To show this, we use Lemma 3 from Arcones (1998). We see that the first condition of Theorem 1 from Arcones (1998) is met, since $h_{4}$ is bounded and the mixing coefficient converges to zero. Therefore, by the fact that $h_{4}$ is canonical, we can use Lemma 3 from Arcones (1998), which states that

for some $p>2$ and $M=\parallel h\parallel_{\infty}$ . Take $p$ such that $\frac{3(p-2)}{p}=2.5$ and use inequality $\beta(m)\leq m^{-3}$ to obtain

We now need to show that $\mathcal{E}nU(h_{4},Z)$ converges to zero. We will use Lemma 4 with $\delta=2$, and that $\beta(k)^{\frac{2}{3}}\leq k^{-2}$,

for some constant $M$ as in Lemma 4. Next

We have used the fact that $\sum_{d=a+3}^{n}\frac{1}{(d-a)^{2}}\leq 2\zeta(2)$.

The reasoning for $U(h_{3},Z)$ is similar. ∎

A.5 Proofs - Borisov & Volodko (2008, Theorem 1)

Let $m$ be the number of arguments of a symmetric kernel $f$. Let one of the following two sets of conditions be fulfilled:

The stationary sequence $X_{i}$ satisfies $\theta$-mixing and

$\sum_{k=1}^{\infty}\phi(k)^{\frac{1}{2}}<\infty$,

$\sup_{i}\mathcal{E}|e_{i}(X_{1})|^{2}<\infty$.

The stationary sequence $X_{i}$ satisfies $\alpha$-mixing. For some $\epsilon>0$ and for an even number $c\geq 2$ the following holds:

$\sup_{i}\mathcal{E}|e_{i}(X_{1})|^{2+\epsilon}\leq\infty$,

$\sum_{k=1}^{\infty}k^{c-2}\alpha^{\epsilon/(c+\epsilon)}(k)<\infty$

where $e_{i}(X_{1})$ are a basis of $L_{2}(X,F)$. Then, for any degenerate kernel $f(t_{1},...,t_{m})\in L_{2}(X_{m},F_{m})$, under conditions

$\sum_{i_{1},...,i_{m}}^{\infty}|f_{i_{1},...,i_{m}}|<\infty$, where $f_{i_{1},...,i_{m}}$ are the coefficient of $f$ in $L_{2}(X_{m},F_{m})$,

$e_{0}=1$ or $\mathcal{E}e_{i}(Z_{j})=0$ for all $i$,

where $\tau_{j}$ is a centred Gaussian sequence with the covariance matrix

$\nu_{j}(i_{1},...,i_{m}):=\sum_{r=1}^{m}\delta_{j,i_{r}}$, and $H_{k}(x)$ are the Hermite polynomials,

A.6 Proofs - Expected value of the eigenfunctions

From the eigenvalue equation $\lambda_{i}\mathcal{E}e_{i}(z)=\mathcal{E}h_{2}(z,Z_{2}^{*})e_{i}(Z_{2}^{*})$, $h_{2}$ degeneracy, and the independence of $Z_{1}^{*}$ and $Z_{2}^{*}$, we conclude that