Kernel Two-Sample and Independence Tests for Non-Stationary Random Processes

Two-sample and independence tests with the kernel-based MMD and HSIC have shown remarkable results on i.i.d. data and stationary random processes. However, these statistics are not directly applicable to non-stationary random processes, a prevalent form of data in many scientific disciplines. In this work, we extend the application of MMD and HSIC to non-stationary settings by assuming access to independent realisations of the underlying random process. These realisations - in the form of non-stationary time-series measured on the same temporal grid - can then be viewed as i.i.d. samples from a multivariate probability distribution, to which MMD and HSIC can be applied. We further show how to choose suitable kernels over these high-dimensional spaces by maximising the estimated test power with respect to the kernel hyper-parameters. In experiments on synthetic data, we demonstrate superior performance of our proposed approaches in terms of test power when compared to current state-of-the-art functional or multivariate two-sample and independence tests. Finally, we employ our methods on a real socio-economic dataset as an example application.


Introduction
Non-stationary processes are the rule rather than the exception in many scientific disciplines such as epidemiology, biology, sociology, economics, or finance. In recent years, there has been a surge of interest in the analysis of problems described by large sets of interrelated variables with few observations over time, often involving complex non-linear and non-stationary behaviours. Examples of such problems include the longitudinal spread of obesity in social networks Christakis and Fowler [2007], disease modelling from time-varying inter-and intra-cellular relationships Barabási et al. [2011], behavioural responses to losses of loved ones within social groups Bond [2017], and the linkage between climate change and the global financial system Battiston et al. [2017]. All such analyses rely on the statistical assessment of the similarity between, or the relationship amongst, noisy time series that exhibit temporal memory. Therefore, the ability to test the statistical significance of homogeneity and dependence between random processes that cannot be assumed to be independent and identically distributed (i.i.d.) is of fundamental importance in many fields.
Kernel-based methods provide a popular framework for homogeneity and independence tests by embedding probability distributions in reproducing kernel Hilbert spaces (RKHSs) [Muandet et al., 2017, Section 2.2]. Of particular interest are the kernel-based two-sample statistic maximum mean discrepancy (MMD) Gretton et al. [2007], which is used to assess whether two samples were drawn from the same distribution, hence testing for homogeneity; and the related Hilbert-Schmidt independence criterion (HSIC) Gretton et al. [2008], which is used to assess dependence between two random variables, thus testing for independence. These methods are non-parametric, i.e., they do not make any assumptions on the underlying distribution or the type of dependence. However, in their original form, both MMD and HSIC assume access to a sample of i.i.d. observations-an assumption that is often violated for temporally-dependent data such as random processes.
Extensions of MMD and HSIC to random processes have been proposed Besserve et al. [2013], Chwialkowski et al. [2014]. Yet, these methods require the random process to be stationary, meaning that its distribution does not change over time. Whilst it is sometimes possible to approximately achieve stationarity with pre-processing techniques such as (seasonal) differencing or square root and power transformations, such approaches become cumbersome and notoriously difficult, particularly with large sets of variables. The stationarity assumption can therefore pose severe limitations in many application areas where multiple non-stationary processes must be taken into consideration. When studying the relationships of climate change to the global financial system, for example, factors such as greenhouse gas emissions, stock market indices, government spending, and corporate profits would have to be transformed or assumed to be stationary over time.
In this paper, we show how the kernel-based statistics MMD and HSIC can be applied to non-stationary random processes. At the heart of our proposed approach is the simple, yet effective idea that realisations of a random process in the form of temporally-dependent measurements (i.e., the observed time series) can be viewed as independent samples from a multivariate probability distribution, provided that they are observed at the same points in time, i.e., over the same temporal grid. Then, MMD and HSIC can be applied on these distributions to test for homogeneity and independence, respectively.
The remainder of this paper is structured as follows. After discussing related work in section 2, we introduce our applications of two-sample and independence testing with MMD and HSIC to non-stationary random processes in section 3. We then carry out experiments on multiple synthetic datasets in section 4 and demonstrate that the proposed tests have higher power compared with current functional or multivariate two-sample and independence tests under the same conditions. We provide an example application of our proposed methods to a socio-economic dataset in section 5 and conclude the paper with a brief discussion in section 6.

Related work
Two-sample and independence tests on stochastic processes have been widely studied in recent years. Under the stationarity assumption, Besserve et al. [2013] investigate how the kernel cross-spectral density operator may be used to test for independence, and Chwialkowski et al. [2014] formulate a wild bootstrap-based approach for both two-sample and independence tests, which outperforms Besserve et al. [2013] in various experiments. The wild bootstrap in Chwialkowski et al. [2014] approximates the null hypothesis H 0 by assuming there exists a time lag τ such that a pair of measurements at any point in time t, (x i , y i ) t , is independent of (x i , y i ) t±s for s ≥ τ . This method is applicable to test for instantaneous homogeneity and independence in stationary processes, but requires further assumptions to investigate non-instantaneous cases: a maximum lag M ≤ τ must be defined as the largest absolute lag for the test. This results in multiple hypothesis testing requiring adjustment by a Bonferroni correction. Further, Davis et al. [2018] have applied distance correlation Székely et al. [2007], a HSIC-related statistic, to independence testing on stationary random processes.
Beyond the stationarity assumption, two-sample testing in the functional data analysis literature has mostly focused on differences of mean  or covariance structures Fremdt et al. [2012], Panaretos et al. [2010]. However, Pomann et al. [2016] have developed a two-sample test for distributions based on generalisations of a finite-dimensional test by utilising functional principal component analysis, and Wynne and Duncan [2020] have derived kernels over functions to be used with MMD for the two-sample test. Independence testing for functional data using kernels was recently proposed in Górecki et al. [2018], but assumes the samples lie on a finite-dimensional subspace of the function space-an assumption not required in our work. Moreover, Zhang et al. [2018] have developed computationally efficient methods to test for independence on high-dimensional distributions and large sample sizes by using eigenvalues of centred kernel matrices to approximate the distribution under the null hypothesis H 0 instead of simulating a large number of permutations.
3 MMD and HSIC for non-stationary random processes

Notation and assumptions
Let {X t } and {Y t } denote two non-stationary stochastic processes with probability laws P X and P Y , respectively. We assume that we observe m independent realisations of {X t } and n independent realisations of {Y t } in the form of time series measured at T X and T Y time points, respectively. Said differently, the data samples ∼ P X are a set of non-stationary time series, x i = {x i,1 , . . . , x i,T X }, arriving over the same temporal grid, and similarly for ∼ P Y with y i = {y i,1 , . . . , y i,T Y }. Note that the measurements x i,t and y i,t are not independent across time. 1 We may view the realisations x i and y i as samples of multivariate probability distributions of dimension T X and T Y , respectively, which are independent at any given point in time, i.e., x i,t ⊥ ⊥ x j,t and y i,t ⊥ ⊥ y j,t ∀t and ∀i = j. Consequently, we can represent these distributions by their mean embeddings µ X and µ Y in RKHSs, and use these to conduct kernel-based two-sample and independence tests. Given a characteristic kernel k, i.e., the mean embedding µ captures all information of a distribution P Sriperumbudur et al. [2010], the dependence between measurements in time is captured by the ordering of the variables, and the fact that any characteristic kernel k is injective, thus guaranteeing a unique mapping of any probability distribution into a RKHS Sriperumbudur et al. [2011].
For homogeneity testing (P X ? = P Y ), we use the kernel-based MMD statistic and require equal number of measurements T = T X = T Y , but allow different sample sizes, m = n. For independence testing (P XY ? = P X P Y ), we employ the related HSIC, and in this case number of measurements can differ, but we require the same number of realisations, m = n. We now describe how two-sample and independence tests can be performed under these assumptions.

MMD for non-stationary random processes
Let k : R T × R T → R be a characteristic kernel, such as the Gaussian kernel k(x, y) = exp (− x − y 2 /σ 2 ), which uniquely maps P X and P Y to its associated RKHS H k via the mean embeddings µ X := k(x, ·) dP X (x) and µ Y := k(y, ·) dP Y (y) [Muandet et al., 2017, Section 2.1]. The MMD between P X and P Y in H k is defined as Gretton et al. [2007]: Given samples X and Y, MMD 2 (H k , P X , P Y ) can then be approximated by the following unbiased estimator Gretton et al. [2007]: Henceforth, we drop the implied H k for ease of notation.
Using MMD 2 u (X, Y) as a test statistic, one can construct a statistical two-sample test for the null hypothesis H 0 : P X = P Y against the alternative hypothesis H 1 : P X = P Y Gretton et al. [2012]. Let α be the significance level of the test, i.e., the maximum allowable probability of falsely rejecting H 0 and hence an upper bound on the type-I error. Given α, the threshold c α for the test statistic can be approximated with a permutation test as follows. We first generate P randomly permuted partitions of the set of all realisations X ∪ Y with sizes commensurate with (X, Y), denoted (X p , Y p ), p = 1, . . . , P . We then compute MMD 2 u (X p , Y p ), ∀p, and sort the results in descending order. Finally, we select the statistic at position (1 − α) × P as our empirical thresholdĉ α . The null hypothesis H 0 is then rejected if MMD 2 u (X, Y) >ĉ α . For a computationally less expensive (but generally less accurate) option, the inverse cumulative density function of the Gamma distribution can be computed to approximate the null distribution Gretton et al. [2009].

HSIC for non-stationary random processes
Let P XY denote the joint distribution of {X t } and {Y t }, and let H k and G l be separable RKHSs with characteristic kernels k : respectively. HSIC is then defined as the MMD between P XY and P X P Y Gretton et al. [2008]: Here, ⊗ denotes the tensor product. Recall that we assume an equal number of realisations m for both processes, and let K, L ∈ R m×m be the kernel matrices with entries k ij = k(x i , x j ) and l ij = l(y i , y j ), respectively. Given i.i.d.
samples (X, Y), an unbiased empirical estimator of HSIC(H k , G l , P XY ) is given by [Song et al., 2012, Theorem 2]: where K = K − diag(K) and L = L − diag(L), and 1 is the m × 1 vector of ones. To ease our notation, we henceforth omit the implied H k and G l .
To test HSIC u (XY) for statistical significance, we define the null hypothesis H 0 : P XY = P X P Y and the alternative H 1 : P XY = P X P Y . We broadly repeat the procedure outlined in section 3.2 by bootstrapping the distribution under H 0 via permutations, with the distinction that we only permute the samples Gretton et al. [2008]. HSIC u (XY) is then computed for each permutation (X, Y p ) and the empirical thresholdĉ α is taken as the statistic at position

Maximising the test power
The power of both MMD-based two-sample and HSIC-based independence tests is prone to decay in high dimensional spaces , , as in our setting where each measurement point in time is treated as a separate dimension. Hence, we describe here how a kernel k can be chosen to maximise the test power, i.e., the probability of correctly rejecting H 0 given that it is false. First, note that under H 1 both MMD 2 u (X, Y) [Gretton et al., 2012, Corollary 16] and HSIC u (XY) [Gretton et al., 2008, Theorem 1] are asymptotically Gaussian: where V MMD m (P X , P Y ) and V HSIC m (P XY ) denote the asymptotic variance of MMD , the test power is defined in terms of P 1 , the distributions under H 1 , with equal sample sizes m = n as: where Φ is the cumulative density function of the standard Gaussian distribution, and whereĉ α → c α with increasing sample size. To maximise the test power, we maximise the argument of Φ, which we approximate by maximising We perform this optimisation by splitting our samples (X, Y) into training and testing sets, of which we take the former to learn the kernel hyper-parameters and the latter to conduct the final hypothesis test with the learnt kernel.

Experimental results on synthetic data
To evaluate our proposed tests empirically, we first apply our homogeneity and independence tests to various nonstationary synthetic datasets. We report test performance usingμ, the percentage of rejection of the null hypothesis H 0 , which becomes the test power once H 0 is false, by repeating the experiments on 200 trials (i.e., 200 independently generated synthetic datasets). We provide 95% confidence intervals computed asμ ± 1.96 μ(1 −μ)/100.

Homogeneity tests with MMD
Setup. We evaluate our MMD-based homogeneity test against shifts in mean and variance of two non-stationary stochastic processes {X t } and {Y t } by establishing if they are correctly accepted or rejected under the null hypothesis H 0 : P X = P Y . For ease of comparison, we adopt the experimental protocol of Pomann et al. [2016] and consider two stochastic processes based on a linear mixed effects model. We generate independent samples X where we set K = 2 with Fourier basis functions φ 1 (t) = √ 2 sin(2πt) and φ 2 (t) = √ 2 cos(2πt). The coefficients ξ Xi,k and ξ Y i,k and the additive noises Xi,t , Y i,t are all independent Gaussian-distributed random variables with means and variances specified below.
The coefficients δ µ and δ σ for mean and variance shifts, respectively, determine the departure from the null hypothesis. Setting δ µ , δ σ = 0 means H 0 is true, whereas δ µ , δ σ > 0 means H 0 is false. Although this is not a necessity, we set the number of independent samples of {X t } and {Y t } to be equal, m = n. To test for statistical significance, we follow the procedure described in section 3.2 and perform permutation tests of P = 5000 partitions for varying values of δ µ and δ σ and different sample sizes m = 100, 200, 300, 500.
Baseline results without test power optimisation. Our baseline results are obtained with a Gaussian kernel k(x, y) = exp (− x − y 2 /σ 2 ) with bandwidth σ equal to the median distance between observations of the aggregated samples. Figure 1 shows how our method (solid lines) compares to Pomann et al. [2016] (dashed lines) for T = 100 discrete time points. For all sample sizes, the type-I error rate lies at or below the allowable probability of false rejection α, and our method significantly outperforms Pomann et al. [2016] for nearly all levels of mean and variance shifts. Both shifts become easier to detect for larger sample sizes. Particularly strong improvements are achieved for mean shifts: our method makes no type-II errors for δ µ ≥ 3 on m = 100 samples, whereas Pomann et al. [2016] only reach such performance with m = 500 samples and δ µ ≥ 4.5. We obtain similar test power results (see Appendix A.1) for coarser realisations with T = 5, 10, 25, 50 over the same interval I = [0, 1].
Results of the optimised test. Next, we apply the method described in section 3.4 to maximise the test power. Specifically, we search for the Gaussian kernel bandwidth σ (over spaces defined in Table 1 in Appendix A.2), that maximises the argument of Φ in our approximations of (7) on our training samples. For demonstrative purposes, we choose to split our dataset equally into training and testing sets although other ratios may lead to higher test power. We find that the test power is significantly improved by our optimisation for the detection of mean shifts. For instance, test power rises fourfold for δ µ = 1 and m = 200 compared to our baseline method. Furthermore, we have no type-II errors once δ µ ≥ 2 for m = 100, as compared to δ µ ≥ 3 for our baseline test and δ µ ≥ 6.5 for Pomann et al. [2016]. In its current form, however, our optimisation does not yield higher test power for the detection of variance shifts, a fact that we discuss in section 6.

Independence tests with HSIC
Setup. To test for independence, the null hypothesis is H 0 : P XY = P X P Y . We assume we observe measurements x i,t and y i,t over temporal grids of length T X and T Y in the interval I = [0, 1], respectively. To measure type-I and type-II error rates, we use the following experimental protocols, partly adopted from Zhang et al. [2018] and Gretton et al. [2008Gretton et al. [ , 2005: • Linear dependence: X is generated as in (9) with µ X (t) = t, basis coefficients ξ Xi,1 ∼ N (0, 10), ξ Xi,2 ∼ N (0, 5), and noise Xi,t ∼ N (0, 0.25). The samples of the second process are Y = {x i,1 + i } m i=1 where i ∼ N (0, 1), as in Zhang et al. [2018]. • Dependence through a shared coefficient: X and Y are generated as in (9) with µ X (t) = µ Y (t) = t and independently sampled ξ Xi,1 , ξ Y i,1 , Xi,t , Y i,t as in the mean shift experiments of section 4.1, but where the stochastic processes now share the second basis function coefficient: ξ Xi,2 = ξ Y i,2 .
Statistical significance is computed using P = 5000 permutations of Y whilst X is kept fixed to approximate the distribution under H 0 . Test power is calculated for varying T = [5,10,25,50,100] and different sample sizes m = n.
Baseline results without test power optimisation. Our baseline results are computed using a Gaussian kernel with σ equal to the median distance between measurements in the corresponding sample. Figure 3 (left) shows the results of our test on the linear dependence experiments, which demonstrate, due to T Y = 1, how dependencies between individual points in time and an entire time series can be detected. We compare our method to: (i) a statistic explicitly aimed at linear dependence, SubCorr = 1 where Corr(·, ·) is the Pearson correlation coefficient; . For both of these methods, the distribution under H 0 is also  Figure 3 (right) displays the power of our independence test for the case of dependent samples through a shared coefficient for varying sample sizes m and measurements T . We compare our results to two spectral methods Zhang et al. [2018] that approximate the distribution under H 0 using eigenvalues of the centred kernel matrices of X and Y: spectral HSIC uses the unbiased estimator (4) as the test statistic with the eigenvalue-based null distribution; and spectral random Fourier feature (RFF) uses a test statistic induced by a number of RFFs (set here to 10) that approximate the kernel matrices of X and Y. Our method and spectral HSIC achieve 20 − 50% improvement in test power compared to spectral RFF. For small numbers of samples (m < 15), our method outperforms spectral HSIC, which converges to the performance of our method with increasing sample size, as we would expect it [Gretton et al., 2009, Theorem 1]. Figure 4 shows the rotation dependence experiments, where θ = 0 corresponds to the null hypothesis (independence) and θ > 0 to the alternative. The distribution hyper-parameters for ξ Xi,k and ξ Y i,k are detailed in Appendix A.3, and we set T X = T Y = T , although equality is not required. As expected, dependence is easier to detect with increasing θ. We observe that denser temporal measurements do not result in enhanced test power. Note that the test power is highly dependent on the distribution of the coefficients of the basis functions ξ Xi,k , ξ Y i,k .
Results of the optimised test. The test power maximisation was applied to the rotation dependence experiments by searching for optimal Gaussian kernel bandwidths σ X and σ Y over pre-defined intervals (specified in Appendix A.2). Figure 4 shows that the test power is improved when the basis function coefficients are drawn from uniform distributions. In this case, the percentage of rejected H 0 is 20 − 40% higher for θ between 0.2 and 0.75 × π/4, but it levels off at 95% once θ ≥ 0.75 × π/4, which is the same level achieved by our baseline method for θ ≥ 0.85 × π/4. With our current test-train split, our optimised test does not improve the test power if the basis function coefficients ξ Xi,k and ξ Y i,k are drawn from student-t or exponential distributions.

Application to a socio-economic dataset
As a further illustration, we apply our method to the United Nations' socio-economic Sustainable Development Goals (SDGs) (see Appendix A.4 for details). Specifically, we investigate whether some so-called Targets of the 17 SDGs have been homogeneous over the last 20 years across low-and high-income countries, and whether certain SDGs in African countries exhibit dependence over the same period. In both setting, we assume countries are independent.
For our homogeneity tests, we classify countries into low-and high-income according to World Bank [2020a]. We use temporal data of 76 Targets for which World Bank [2020b] provides data collected over the last T = 20 years for m = 30 low-income countries and n = 55 high-income countries. Applying our baseline method without test power optimisation, we find that, out of the 76 Targets we have data available for, only 38 have had homogeneous trajectories in low-and high-income countries. For instance, whereas the 'death rate due to road traffic injuries' (Target 3.6) has been homogeneous between these two groups, the 'fight the epidemics of AIDS, tuberculosis, malaria and others' (Target 3.3) has not been homogeneous in low-and high-income countries.
For our independence tests, we consider temporal data from m = n = 49 African countries over T = 20 years, and test any two Targets for pairwise independence. Of the total 2850 possible pairwise combinations, the null hypothesis of independence is rejected for 357. As an illustration, we examine the dependencies of 'implementation of national social protection systems' (Target 1.3) with 'economic growth' (Target 8.1) and the 'proportion of informally employed workers' (Target 8.3). Applying our baseline method, we accept the null hypothesis of independence between Target 1.3 and 8.1, i.e., we find that the 'implementation of national social protection systems' has been independent of economic growth. In contrast, we find that Target 1.3 has been dependent on the 'proportion of informally employed workers' (Target 8.3).

Discussion and conclusion
Building on ideas from functional data analysis, we have presented approaches to testing for homogeneity and independence between two non-stationary random processes with the kernel-based statistics MMD and HSIC. We view independent realisations of the underlying processes as samples from multivariate probability distributions to which MMD and HSIC can be applied. Our tests are shown to outperform current state-of-the-art methods in a range of experiments. Furthermore, we optimise the test power over the choice of kernel and achieve improved results in most settings. However, we also observe that our optimisation procedure does not always yield an increase in test power. We leave the investigation of this behaviour open for future research with the possibility of defining search spaces and step sizes over kernel hyper-parameters differently, or of choosing a gradient-based approach for optimisation Sutherland et al. [2016]. Our results show that small sample sizes of less than 40 independent realisations can already achieve high test power, and that denser measurements over the same time period do not necessarily lead to enhanced test power.

A Appendix
A.1 Results for realisations with varying number of time points, T

A.1.1 MMD
We show here the results for mean and variance shifts for m = n = 100, but the results are similar for all tested sample sizes m = n = 100, 200, 300, 500,

A.1.2 HSIC
Experiments for linear dependence and dependence through shared second basis function coefficient for various T . We find that the granularity of measurements over time does not influence the text power significantly.  Table 1). These search spaces resulted from extensive manual explorations for all shifts and sample sizes. We acknowledge that the test power may be further improved with search spaces of finer granularity.

A.2.2 HSIC
We define search intervals of both σ X and σ Y across all angles θ, but different for the student-t, uniform, and exponential distributions. For student-t and exponential distributions, both σ X and σ Y were chosen as 20 evenly spaced numbers on a linear scale between 1 and 20. For uniform distributions, both σ X and σ Y were chosen as 40 evenly spaced numbers on a linear scale between 1 and 40. These search spaces resulted from extensive manual explorations for all angles and distributions. We acknowledge that the test power may be further improved with search spaces of finer granularity.  [Gretton et al., 2005, Table 3], and Z is a proxy for both X and Y.

A.4 SDG dataset
Data of the Indicators measuring the progress of the Targets of the SDGs can be found at World Bank [2020b]. Each of these Indicators measures the progress towards a specific Target. For instance, an Indicator for Target 1.1, 'by 2030, eradicate extreme poverty for all people everywhere, currently measured as people living on less than $1.90 a day', is the 'proportion of population below the international poverty line, by gender, age, employment status and geographical location (urban/rural)'. Each of the Targets belongs to one specific Goal (e.g., Target 1.1 belongs to Goal 1, 'end poverty in all its forms everywhere'). There are 17 such Goals, which are commonly referred to as the Sustainable Development Goals (SDGs). We compute averages over all Indicators belonging to one Target for our analyses in Section 5.
The dataset of World Bank [2020b] has many missing values, especially for the time span 2000-2005. We impute these values using a weighted average across countries (where data is available) with weights inversely proportional to the Euclidean distance between indicators.