New HSIC-based tests for independence between two stationary multivariate time series

New HSIC-based tests for independence between two stationary multivariate time series

Introduction

Before applying any sophisticated method to describe relationships between two time series, it is important to check whether they are independent or not. If they are dependent, causal analysis techniques, such as copula and multivariate modeling, can be used to investigate the relationship between them, and this may lead to interesting insights or effective predictive models; otherwise, one should analyze them using two independent parsimonious models; see, e.g., Pierce (1977), Schwert (1979), Hong (2001a), Lee and Long (2009), Shao (2009), and Tchahou and Duchesne (2013) for many empirical examples in this context.

Most of the existing methods for testing the independence between two multivariate time series models use a measure based on cross-correlations. Specifically, they aim to check whether the sample cross-correlations of model residuals, up to either certain fixed lag or all valid lags, are significantly different from zeros. The former includes the portmanteau tests (Cheung and Ng, 1996; El Himdi and Roy, 1997; Pham et al. 2003; Hallin and Saidi, 2005 and 2007; Robbins and Fisher, 2015), and the latter with the aid of kernel smooth technique falls in the category of spectral tests (Hong, 2001a and 2001b; Bouhaddioui and Roy, 2006). It must be noted that the idea of using the cross-correlations is a natural extension of the pioneered studies in Haugh (1976) and Hong (1996) for univariate time series models, but in many circumstances it only suffices to convey evidence of uncorrelatedness rather than independence.

Generally speaking, all of the aforementioned tests are designed for investigating the linear dependence (i.e., the cross-correlation in the mean, variance or higher moments) between two model residuals, and hence they could exhibit a lack of power in detecting the non-linear dependence structure. A significant body of research so far has documented the non-linear dependence relationship among a myriad of economic fundamentals; see, e.g., Hiemstra and Jones (1994), Wang et al. (2013), Choudhry et al. (2016), and Diks and Wolski (2016) to name a few. However, less attempts have been made in the literature to account for both linear and nonlinear dependence structure, which shall be two parallel important characteristics to be tested.

To examine the general dependence structure, a direct measure on independence is expected for testing purpose. In the last decade, the Hilbert-Schmidt independence criterion (HSIC) in Gretton et al. (2005) has been extensively used in many fields. Some inspiring works in one- or two-sample independence tests via HSIC include Gretton et al. (2008) and Gretton and Gy\"{o}rfi (2010) for observable i.i.d. data, and Zhang et al. (2009), Zhou (2012) and Fokianos and Pitsillou (2017) for observable dependent or time series data. The last two instead applied the distance covariance (DC) in Sz\'{e}kely et al. (2007), while Sejdinovic et al. (2013) showed that HSIC and DC are equivalent. When the data are un-observable and derived from a fitted statistical model (e.g., the estimated model innovations), the estimation effect has to be taken into account. The original procedure based on HSIC or DC will no longer be valid, and a modification of the above procedure has to be derived for testing purpose. By now, very little work has been done in this context. Two exceptions are Sen and Sen (2014) and Davis et al. (2016) for one-sample independence tests; the former focused on the regression model with independent covariates, and the latter considered the vector AR models but without providing a rigorous way to obtain the critical values of the related test.

This paper proposes some novel one-sided tests for the independence between two stationary multivariate time series. These new tests apply the HSIC to examine the independence between the un-observable innovation vectors of both time series. Among them, the single HSIC-based test is tailored to detect the general dependence between these two innovation vectors at a specific lag $m$, and the joint HSIC-based test is designed for this purpose up to certain lag $M$. Under regular conditions, the limiting null distributions of our HSIC-based tests are established. Next, our HSIC-based tests are shown to be consistent. Moreover, a residual bootstrap method is used to obtain the critical values for our HSIC-based tests, and its validity is justified. Our methodologies are applicable for the general specifications of the time series models driven by i.i.d. innovations. By choosing different lags, our new tests can give investigators more complete information on the general (including both linear and non-linear) dependence relationship between two time series. Finally, the importance of our HSIC-based tests is illustrated by some simulation results and a real example.

This paper is organized as follows. Section 2 introduces our HSIC-based test statistics and some technical assumptions. Section 3 studies the asymptotic properties of our HSIC-based tests. A residual bootstrap method is provided in Section 4. Simulation results are reported in Section 5. One real example is presented in Section 6. Concluding remarks are offered in Section 7. The proofs are provided in the Appendix.

Throughout the paper, $\mathcal{R}=(-\infty,\infty)$, $C$ is a generic constant, $I_{s}$ is the $s\times s$ identity matrix, $1_{s}$ is the $s\times1$ vector of ones, $\otimes$ is the Kronecker product, $A^{T}$ is the transpose of matrix $A$, $\|A\|$ is the Euclidean norm of matrix $A$, $vec(A)$ is the vectorization of $A$, $vech(A)$ is the half vectorization of $A$, $D(A)$ is the diagonal matrix whose main diagonal is the main diagonal of matrix $A$, $\partial_{x}h$ denotes the partial derivative with respect to $x$ for any function $h(x,y,\cdots)$, $o_{p}(1)(O_{p}(1))$ denotes a sequence of random numbers converging to zero (bounded) in probability, “$\to_{d}$” denotes convergence in distribution, and “$\to_{p}$” denotes convergence in probability.

The HSIC-based test statistics

Review of the Hilbert-Schmidt Independence Criterion

In this subsection, we briefly review the Hilbert-Schmidt independence criterion (HSIC) for testing the independence of two random vectors; see, e.g., Gretton et al. (2005) and Gretton et al. (2008) for more details.

Let $\mathcal{U}$ be a metric space, and $k: \mathcal{U}\times \mathcal{U}\to \mathcal{R}$ be a symmetric and positive definite (i.e., $\sum_{i,j}c_{i}c_{j}k(x_{i},x_{j})\geq0$ for all $c_{i}\in\mathcal{R}$) kernel function. There exists a Hilbert space $\mathcal{H}$ (called {\it Reproducing Kernel Hilbert Space} (RKHS)) of functions $f: \mathcal{U}\to \mathcal{R}$ with inner product $\langle \cdot,\cdot \rangle$ such that

For any Borel probability measure $P$ defined on $\mathcal{U}$, its {\it mean element} $\mu[P]\in\mathcal{H}$ is defined as follows:

where the random variable $U\sim P$. From ((ref))-((ref)), we have $\mu[P](u)=\langle k(\cdot,u),\mu[P]\rangle=E[k(U,u)]$. Furthermore, we say that $\mathcal{H}$ is {\it characteristic} if and only if the map $ P\to \mu[P] $ is injective on the space $\mathcal{P}:=\{P:\int_{\mathcal{U}}k(u,u)dP(u)<\infty\}$.

Likewise, let $\mathcal{G}$ be a second RKHS on a metric space $\mathcal{V}$ with kernel $l$. Let $P_{uv}$ be a Borel probability measure defined on $\mathcal{U}\times\mathcal{V}$, and let $P_{u}$ and $P_{v}$ denote the respective marginal distributions on $\mathcal{U}$ and $\mathcal{V}$, respectively. Assume that

where the random variable $(U,V)\sim P_{uv}$. The HSIC of $P_{uv}$ is defined as

where $(U',V')$ is an i.i.d. copy of $(U,V)$, and $E_{\xi,\zeta}$ (or $E_{\xi}$) denotes the expectation over $(\xi,\zeta)$ (or $\xi$). Following Sejdinovic et al. (2013), if ((ref)) holds and both $\mathcal{H}$ and $\mathcal{G}$ are characteristic, then $$\Pi(U,V)=0\,\,\,\,\mbox{ if and only if }\,\,\,\, P_{uv}=P_{u}\times P_{v}.$$ Therefore, we can test the independence of $U$ and $V$ by examining whether $\Pi(U,V)$ is significantly different from zero.

Suppose the samples $\{(U_{i},V_{i})\}_{i=1}^{n}$ are from $P_{uv}$. Following Gretton et al. (2005), the empirical estimator of $\Pi(U,V)$ is

where $k_{ij}=k(U_{i},U_{j})$, $l_{ij}=l(V_{i},V_{j})$, $K=(k_{ij})$ and $L=(l_{ij})$ are $n\times n$ matrices with entries $k_{ij}$ and $l_{ij}$, respectively, and $H=I_{n}-(1_{n}1_{n}^{T})/n$. Here, each index of the summation $\sum$ is taken from 1 to $n$. If $\{(U_{i},V_{i})\}_{i=1}^{n}$ are i.i.d. samples, Gretton et al. (2005) showed that $\Pi_{n}$ is a consistent estimator of $\Pi(U,V)$.

In order to compute $\Pi_{n}$, we need to choose the kernel functions $k$ and $l$. In the sequel, we assume $\mathcal{U}=\mathcal{R}^{\kappa_{1}}$ and $\mathcal{V}= \mathcal{R}^{\kappa_{2}}$ for two positive integers $\kappa_{1}$ and $\kappa_{2}$. Then, some well known choices (see Peters, 2008; Zhang et al. 2017) for $k$ (or $l$) are given below:

We shall highlight that the HSIC is easy-to-implement in multivariate cases, since the computation cost of $\Pi_{n}$ is $O(n^{2})$ regardless of the dimensions of $U$ and $V$, and many softwares can calculate ((ref)) very fast.

Test statistics

Consider two multivariate time series $Y_{1t}$ and $Y_{2t}$, where $Y_{1t}\in\mathcal{R}^{d_{1}}$ and $Y_{2t}\in\mathcal{R}^{d_{2}}$. Assume that each $Y_{st}$ ($s=1$ or $2$ hereafter) admits the following specification:

where $I_{st}=(Y_{st}^{T},Y_{st-1}^{T},\cdots)^{T}\in\mathcal{R}^{\infty}$ is the information set at time $t$, $\theta_{s0}\in\mathcal{R}^{p_{s}}$ is the true but unknown parameter value of model ((ref)), $\eta_{st}\in\mathcal{R}^{d_{s}}$ is a sequence of i.i.d. innovations such that $\eta_{st}$ and $\mathcal{F}_{st-1}$ are independent, $\mathcal{F}_{st}:=\sigma(I_{st})$ is a sigma-field, and $f_{s}: \mathcal{R}^{\infty}\times\mathcal{R}^{p_{s}}\times\mathcal{R}^{d_{s}}\to \mathcal{R}^{d_{s}}$ is a known measurable function. Model ((ref)) is rich enough to cover many often used models, e.g., the vector AR model in Sim (1980), the BEKK model in Engle and Kroner (1995), the dynamic correlation model in Tse and Tsui (2002), and the vector ARMA-GARCH model in Ling and McAleer (2003) to name a few; see also L\"{u}tkepohl (2005), Bauwens et al. (2006), Silvennoinen and Ter\"{a}svirta (2008), Francq and Zako\"{i}an (2010), and Tsay (2014) for surveys.

Model ((ref)) ensures that each $Y_{st}$ admits a dynamical system generated by the innovation sequence $\{\eta_{st}\}$. A practical question is whether either one of the dynamical systems should include the information from the other one, and this is equivalent to testing the null hypothesis:

If $H_{0}$ is accepted, we can separately study these two systems; otherwise, we may use the information of one system to get a better prediction of the other system. Let $m$ be a given integer. Most of the conventional testing methods for $H_{0}$ in ((ref)) aim to detect the linear dependence between $\eta_{1t}$ and $\eta_{2t+m}$ (or their higher moments) via their cross-correlations. Below, we apply HSIC to examine the general dependence between $\eta_{1t}$ and $\eta_{2t+m}$.

To introduce our HSIC-based tests, we need some more notations. Let $\theta_{s}=(\theta_{s1},\theta_{s2},$ $\cdots,\theta_{sp_{s}})\in\Theta_{s}\subset\mathcal{R}^{p_{s}}$ be the unknown parameter of model ((ref)), where $\Theta_{s}$ is a compact parametric space. Assume that $\theta_{s0}$ is an interior point of $\Theta_{s}$, and $Y_{st}$ admits a causal representation, i.e.,

where $g_{s}: \mathcal{R}^{d_{s}}\times\mathcal{R}^{\infty}\times\mathcal{R}^{p_{s}}\to\mathcal{R}^{d_{s}}$ is a measurable function. Moreover, based on the observations $\{Y_{st}\}_{t=1}^{n}$ and (possibly) some assumed initial values, we let

be the residual of model ((ref)), where $\widehat{\theta}_{sn}$ is an estimator of $\theta_{s0}$, and $\widehat{I}_{st}$ is the observed information set up to time $t$.

As for ((ref))-((ref)), our single HSIC-based test statistic on $\widehat{\eta}_{1t}$ and $\widehat{\eta}_{2t+m}$ is

for $m\geq0$, where $\widehat{k}_{ij}=k(\widehat{\eta}_{1i},\widehat{\eta}_{1j})$, $\widehat{l}_{ij}=l(\widehat{\eta}_{2i+m},\widehat{\eta}_{2j+m})$, and $\widehat{K}=(\widehat{k}_{ij})$ and $\widehat{L}=(\widehat{l}_{ij})$ are $N\times N$ matrices with entries $\widehat{k}_{ij}$ and $\widehat{l}_{ij}$, respectively. Here, the effective sample size $N=n-m$, and each index of the summation is taken from $1$ to $N$. Likewise, our single HSIC-based test statistic on $\widehat{\eta}_{1t+m}$ and $\widehat{\eta}_{2t}$ is

for $m\geq0$. Clearly, $S_{1n}(0)=S_{2n}(0)$.

With the help of the single HSIC-based test statistics, we can further define the joint HSIC-based test statistics as follows:

for some specified integer $M\geq0$. The joint test statistic $J_{1n}(M)$ or $J_{2n}(M)$ can detect the general dependence structure of two innovations up to certain lag $M$, while the single test statistic $S_{1n}(m)$ or $S_{2n}(m)$ is used to examine the general dependence structure of two innovations at a specific lag $m$.

Asymptotic theory

This section studies the asymptotics of our HSIC-based test statistics $S_{1n}(m)$ and $J_{1n}(M)$. The asymptotics of $S_{2n}(m)$ and $J_{2n}(M)$ can be derived similarly, and hence the details are omitted for simplicity.

Technical conditions

To derive our asymptotic theory, the following assumptions are needed.

We offer some remarks on the above assumptions. Assumption (ref) is standard for time series models. Assumption (ref)(i) requires some technical moment conditions for the partial derivatives of $g_{st}$. Assumption (ref)(ii) presents some temporal dependence condition on the joint sequence $\{(\eta_{1t}^{T},\eta_{2t}^{T})^{T}\}$. Assumption (ref) is satisfied under mild conditions for most estimators, such as (quasi) maximum likelihood estimator (MLE), least squares estimator (LSE), nonlinear least squares estimator (NLSE) and their robust modifications; see, e.g., Comte and Lieberman (2003), L\"{u}tkepohl (2005), and Hafner and Preminger (2009) for more details. Assumption (ref) is a condition on the truncation of the information set $\widehat{I}_{st-1}$ and is similar to Assumption A5 in Escanciano (2006). Assumption (ref) gives some restrictive conditions for kernel functions $k$ and $l$; these conditions may exclude some kernel functions such as the fractional Brownian motion kernel, but they are usually satisfied by the often used Gaussian kernel, Laplace kernel and inverse multi-quadratics kernel. The conditions in Assumptions (ref)-(ref) may be further relaxed, but they are convenient for presenting our proofs in a simple way.

Some lemmas

This subsection gives some useful lemmas, which are key to study the asymptotics of our test statistics.

Before introducing these lemmas, we present some notations. Let

for $i, j, q, r\in\{1, 2, \cdots, N\}$. With these notations, define

for $a\in\{1,2\}$ and $b\in\{1,\cdots,a+1\}$, where $k_{ij}^{(11)}=\overline{k}_{ij}$, $l_{ij}^{(11)}=l_{ij}$, $k_{ij}^{(12)}=k_{ij}$, $l_{ij}^{(12)}=\overline{l}_{ij}$, $k_{ij}^{(21)}=\widecheck{k}_{ij}$, $l_{ij}^{(21)}=l_{ij}$, $k_{ij}^{(22)}=k_{ij}$, $l_{ij}^{(22)}=\widecheck{l}_{ij}$, $k_{ij}^{(23)}=\overline{k}_{ij}$, and $l_{ij}^{(23)}=\overline{l}_{ij}^{T}$. Then, $S_{1n}^{(0)}(m)$ can be expressed as the $V$-statistic of the form (see Gretton et al. 2005):

for some symmetric kernel $h^{(0)}_{m}$ given by

where the sum is taken over all $4!$ permutations of $(i,j,q,r)$, and $\eta_{t}^{(m)}=(\eta_{1t},\eta_{2t+m})\in\mathcal{R}^{d_{1}}\times\mathcal{R}^{d_{2}}$. Likewise, all $S_{1n}^{(ab)}(m)$ can be expressed as the $V$-statistics for the symmetric kernel $h^{(ab)}_{m}$ given by

where the sum is taken over all $4!$ permutations of $(i,j,q,r)$, and $$\varsigma_{t}^{(m)}=\left(\eta_{1t},\frac{\partial g_{1t}(\theta_{10})}{\partial\theta_{1}},\eta_{2t+m},\frac{\partial g_{2t+m}(\theta_{20})}{\partial\theta_{2}}\right) \in \mathcal{R}^{d_{1}}\times \mathcal{R}^{p_{1}\times d_{1}}\times \mathcal{R}^{d_{2}}\times \mathcal{R}^{p_{2}\times d_{2}}.$$

Now, we are ready to introduce these three lemmas. The first lemma below gives an important expansion of $S_{1n}(m)$.

The second lemma below is crucial in studying the asymptotics of $S_{1n}^{(0)}(m)$ and $S_{1n}^{(ab)}(m)$ under $H_{0}$.

By standard arguments for V-statistics (see, e.g., Lee (1990)), we have $N[S_{1n}^{(0)}(m)]=N[V_{1n}^{(0)}(m)]+o_{p}(1)$, where

is the $V$-statistic with the kernel function

for $x_{1}, x_{2}\in\mathcal{R}^{d_{1}} \times\mathcal{R}^{d_{2}}$. Under $H_{0}$, $\{\eta_{t}^{(m)}\}$ is a sequence of i.i.d. random variables, and hence Lemma (ref)(i) implies that $V_{1n}^{(0)}(m)$ is a degenerate $V$-statistic of order 1, from which $h_{2m}^{(0)}$ can be expressed as

where $\{\Phi_{jm}(\cdot)\}$ is an orthonormal function in $L_{2}$ norm, and $\lambda_{jm}$ is the eigenvalue corresponding to the eigenfunction $\Phi_{jm}(\cdot)$. That is, $\{\lambda_{jm}\}$ is a finite enumeration of the nonzero eigenvalues of the equation $$E[h_{2m}^{(0)}(x_{1},\eta_{1}^{(m)})\Phi_{jm}(\eta_{1}^{(m)})]=\lambda_{jm}\Phi_{jm}(x_{1}),$$ where $E\Phi_{jm}(\eta_{1}^{(m)})=0$ for all $j\geq1$, and $$E[\Phi_{jm}(\eta_{1}^{(m)})\Phi_{j'm}(\eta_{1}^{(m)})]= \left\{

\right.$$ (see, e.g., Dunford and Schwartz (1963, p.1087)). With (\ref{3.7}) and (\ref{3.8}), we can obtain that under $H_{0}$,

Moreover, we consider $S_{1n}^{(ab)}(m)$, which results from the estimation effect. Under $H_{0}$, $S_{1n}^{(ab)}(m)$ (for $a, b=1, 2$) is a degenerate $V$-statistic of order 1 by Lemma (ref)(ii), and hence $N[S_{1n}^{(ab)}(m)]=O_{p}(1)$, and then its related estimation effect is negligible in view of that $\zeta_{sn}^{T}N[S_{1n}^{(ab)}(m)]=o_{p}(1)$. However, under $H_{0}$, the estimation effect related to $S_{1n}^{(23)}(m)$ is negligible only when $\Upsilon=0$. This is because when $\Upsilon\not=0$, $S_{1n}^{(23)}(m)=O_{p}(1)$ by the law of large numbers for V-statistics, and its related estimation effect is not negligible based on the ground that $N[\zeta_{1n}^{T}S_{1n}^{(23)}(m)\zeta_{2n}]=O_{p}(1)$.

Our third lemma below provides a useful central limit theorem.

Asymptotics of test statistics

Based on Lemmas (ref)-(ref), this subsection studies the asymptotics of our test statistics. Let

First, we give the limiting null distributions of $S_{1n}(m)$ and $J_{1n}(M)$ as follows.

Theorem (ref) shows that $S_{1n}(m)$ and $J_{1n}(M)$ have convergence rate $n^{-1}$ under $H_{0}$. Based on this theorem, we reject $H_{0}$ at the significance level $\alpha$, if $$n[S_{1n}(m)]>c_{m\alpha}\,\,\,\mbox{ or }\,\,\,n[J_{1n}(M)]>c_{\alpha},$$ where $c_{m\alpha}$ and $c_{\alpha}$ are the $\alpha$-th upper percentiles of $\chi_{m}$ and $\sum_{m=0}^{M}\chi_{m}$, respectively. Since the distribution of $\chi_{m}$ depends on $\{Y_{st}\}$ and $\{\pi_{st}\}$, a residual bootstrap method is proposed in Section 4 to obtain the values of $c_{m\alpha}$ and $c_{\alpha}$.

Second, we study the behavior of $S_{1n}(m)$ under the following fixed alternative:

Under $H_{1}^{(m)}$, $h_{2m}^{(0)}$ is not a degenerate kernel of order 1. Hence, the V-statistic $S_{1n}^{(0)}(m)$ can not have the convergence rate $n^{-1}$ as suggested by Lemma (ref)(i), leading to the consistency of $S_{1n}(m)$ in detecting $H_{1}^{(m)}$. Similarly, we can show the consistency of $J_{1n}(M)$ to detect the fixed alternative below:

In the end, we highlight that similar results as in Theorems (ref)-(ref) hold for $S_{2n}(m)$ and $J_{2n}(M)$, which can be implemented in the similar way as $S_{1n}(m)$ and $J_{1n}(M)$, respectively.

Residual bootstrap approximations

In this section, we introduce a residual bootstrap method to approximate the limiting null distributions in Theorem (ref). The residual bootstrap method has been well used in the time series literature; see, e.g., Berkowitz and Kilian (2000), Paparoditis and Politis (2003), Politis (2003), and many others. Our residual bootstrap procedure to obtain the approximation of the critical values $c_{m\alpha}$ and $c_{\alpha}$ is as follows:

{\it Step 1}. Estimate the original model ((ref)) and obtain the residuals $\{\widehat{\eta}_{st}\}_{t=1}^{n}$.

{\it Step 2}. Generate bootstrap innovations $\{\widehat{\eta}_{st}^{*}\}_{t=1}^{n}$ (after standardization) by resampling with replacement from the empirical residuals $\{\widehat{\eta}_{st}\}_{t=1}^{n}$.

{\it Step 3}. Given $\widehat{\theta}_{sn}$ and $\{\widehat{\eta}_{st}^{*}\}_{t=1}^{n}$, generate bootstrap data set $\{Y_{st}^{*}\}_{t=1}^{n}$ according to $$Y_{st}^{*}=f_{s}(\widehat{I}_{st-1}^{*},\widehat{\theta}_{sn},\widehat{\eta}_{st}^{*}),$$ where $\widehat{I}_{st}^{*}$ is the bootstrap observable information set up to time $t$, conditional on some assumed initial values.

{\it Step 4}. Based on $\{Y_{st}^{*}\}_{t=1}^{n}$, compute $\widehat{\theta}_{sn}^{*}$ in the same way as for $\widehat{\theta}_{sn}$, and then calculate the corresponding bootstrap residuals $\{\widehat{\eta}_{st}^{**}\}_{t=1}^{n}$ with $\widehat{\eta}_{st}^{**}:=g_{s}(Y_{st}^{*},\widehat{I}_{st-1}^{*},\widehat{\theta}_{sn}^{*})$.

{\it Step 5}. Calculate the bootstrap test statistic $S_{1n}^{**}(m)$ and $J_{1n}^{**}(M)$ in the same way as for ((ref)) and ((ref)), respectively, with $\widehat{\eta}_{st}^{**}$ replacing $\widehat{\eta}_{st}$.

{\it Step 6}. Repeat steps 1-5 $B$ times to obtain $\{n[S_{1nb}^{**}(m)]; b=1, 2, \cdots, B\}$ and $\{n[J_{1nb}^{**}(M)]; b=1, 2, \cdots, B\}$, then choose their $\alpha$-th upper percentiles, denoted by $c_{\alpha}^{*}$ and $c_{m\alpha}^{*}$, as the approximations of $c_{\alpha}$ and $c_{m\alpha}$, respectively.

In order to prove the validity of the bootstrap procedure in steps 1-6, we need some notations. Let

where $\widehat{\eta}_{t}^{(m*)}=(\widehat{\eta}_{1t}^{*}, \widehat{\eta}_{2t+m}^{*})$ and $ {\varsigma}_{t}^{(m*)}=\big(\widehat{\eta}_{1t}^{*},\frac{\partial g_{1t}(\widehat{\theta}_{1n})}{\partial\theta_{1}},\widehat{\eta}_{2t+m}^{*},\frac{\partial g_{2t+m}(\widehat{\theta}_{2n})}{\partial\theta_{2}}\big). $ Also, let $\zeta_{sn}^*=\widehat{\theta}_{sn}^{*}-\widehat{\theta}_{sn}$, and $\varpi_{n}:=\{Y_{11},Y_{12},\cdots,Y_{1n},$ $Y_{21}, Y_{22},\cdots,Y_{2n}\}$ be the given sample. Denote by $E^{*}$ the expectation conditional on $\varpi_{n}$; by $o_{p}^{*}(1) (O_{p}^{*}(1))$ a sequence of random variables converging to zero (bounded) in probability conditional on $\varpi_{n}$.

Since $\{\widehat{\eta}_{st}^*\}_{t=1}^N$ is an i.i.d sequence conditional on $\varpi_{n}$, a similar argument as for Lemma (ref) implies that

where $S_{1n}^{(0*)}(m)$, $S_{1n}^{(ab*)}(m)$ and $R_{1n}^*(m)$ are defined in the same way as $S_{1n}^{(0)}(m)$, $S_{1n}^{(ab)}(m)$ and $R_{1n}(m)$, respectively, with $\eta_{t}^{(m)}$ and $\varsigma_{t}^{(m)}$ being replaced by $\widehat{\eta}_{t}^{(m*)}$ and $\widehat{\varsigma}_{t}^{(m*)}$, respectively. Moreover, by a similar argument as for Lemma (ref)(i), we can obtain

where $E^{*}\Phi_{jm}^{*}(\widehat{\eta}_{1}^{(m*)})=0$ for all $j\geq1$, and $E^{*}[\Phi_{jm}^{*}(\widehat{\eta}_{1}^{(m*)})\Phi_{j'm}^{*}(\widehat{\eta}_{1}^{(m*)})]=1$ if $j=j'$, and 0 if $j\not=j'$.

Next, we give two technical assumptions.

Assumptions (ref) and (ref) are standard to prove the validity of the bootstrap procedure, and they are similar to those in Assumption A7 of Escanciano (2006). For the (quasi) MLE, LSE and NLSE or, more generally, estimators resulting from a martingale estimating equation (see Heyde, 1997), the function $\pi_s(\cdot)$ required in Assumption (ref) could be expressed as $\pi_s(Y_{st},{I}_{st-1},{\theta}_{s})=\varrho_1(\eta_{st}(\theta_{s}))\times \varrho_2({I}_{st-1},\theta_{s})$ for some functions $\varrho_1(\cdot)$ and $\varrho_2(\cdot)$ with $E(\varrho_1(\eta_{st}(\theta_{s0})))=0$. Then, in those cases, Assumptions (ref) and (ref) are satisfied under some mild conditions on the function $\varrho_2(\cdot)$. Note that the calculation of the bootstrap estimator $\widehat{\theta}_{sn}^{*}$ in step 4 may be time-consuming for some times series models (e.g, multivariate ARCH-type models) when $n$ is large. In view of Assumption (ref), we suggest to generate $\widehat{\theta}_{sn}^{*}$ as $$\widehat{\theta}_{sn}^{*}=\widehat{\theta}_{sn}+\frac{1}{n}\sum_{t}\pi_{s}(Y_{st}^{*},\widehat{I}_{st-1}^{*},\widehat{\theta}_{sn}).$$ This results in saving a lot of compute time. In Section 5, we will apply this method to the conditional variance models, and find that it can generate very precise critical values $c_{m\alpha}$ and $c_{\alpha}$ for the proposed HSIC-based tests.

The following theorem guarantees that when $B$ is large, our bootstrapped critical values $c_{m\alpha}$ and $c_{\alpha}$ from steps 1-6 are valid under the null or the alternative hypothesis.

Simulation studies

In this section, we compare the performance of our HSIC-based tests $S_{sn}(m)$ and $J_{sn}(M)$ $(s=1, 2 \mbox{ hereafter})$ with some well-known existing tests in finite samples.

Conditional mean models

We generate 1000 replications of sample size $n$ from the following two conditional mean models:

where $\{\eta_{1t}\}$ and $\{\eta_{2t}\}$ are two sequences of i.i.d. random vectors. To generate $\{\eta_{1t}\}$ and $\{\eta_{2t}\}$, we need an auxiliary sequence of i.i.d. multivariate normal random vectors $\{u_{t}\}$ with mean zero, where $u_{t}=(u_{1t}, u_{2t}, u_{3t}', u_{4t}')'$ with $u_{1t}, u_{2t}\in\mathcal{R}$ and $u_{3t}, u_{4t}\in\mathcal{R}^{2\times1}$, and its covariance matrix is given by

with

Here, we set $\rho_{2}=0.5$ and $\rho_{3}=0.75$ as in El Himdl and Roy (1997), which have also considered model ((ref)) in their simulations.

Based on $\{u_{t}\}$, we consider six different error generating processes (EGPs):

Clearly, each entry of $\eta_{1t}$ or $\eta_{2t}$ has mean zero and variance one. Let $\rho_{\eta_{1},\eta_{2}}(d)$ be the cross-correlation matrix between $\eta_{1t}$ and $\eta_{2t+d}$. EGP 1 is designed for the null hypothesis, since $\rho_{\eta_{1},\eta_{2}}(d)=0_{2\times2}$ for all $d$ in this case. EGPs 2-6 are set for the alternative hypotheses, since they pose a linear or non-linear dependence structure between $\eta_{1t}$ and $\eta_{2t}$. Specifically, a linear dependence structure between $\eta_{1t}$ and $\eta_{2t}$ exists in EGP 2, with $\rho_{\eta_{1},\eta_{2}}(d)=0.3I_{2}$ for $d=0$, and 0 otherwise; a non-linear dependence structure between $\eta_{1t}$ and $\eta_{2t}$ is induced by the co-factor $u_{1t}$ in EGP 3, the lagged co-factors $u_{1t}$ and $u_{1t+3}$ in EGP 4, and two correlated co-factors $u_{1t}$ and $u_{2t}$ in EGPs 5 and 6. In EGPs 2-6, $\eta_{1t}$ and $\eta_{2t}$ are dependent but un-correlated.

Now, we fit each replication by using the least squares estimation method for model ((ref)). Denote by $\{\widehat{\eta}_{1t}\}$ and $\{\widehat{\eta}_{2t}\}$ the residuals from the fitted models. Based on $\{\widehat{\eta}_{1t}\}$ and $\{\widehat{\eta}_{2t}\}$, we compute $S_{sn}(m)$ and $J_{sn}(M)$ ($S_{sn}$ and $J_{sn}$ in short), with $k$ and $l$ being the Gaussian kernels and $\sigma=1$. The critical values of all HSIC-based tests are obtained by the residual bootstrap method with $B=1000$ in Section 4.

Meanwhile, we also compute the test statistics $G_{sn}(M)$ ($G_{sn}$ in short) in El Himdl and Roy (1997) and the test statistics $W_{sn}(h)$ ($W_{sn}$ in short) in Bouhaddioui and Roy (2006), where

Here, $\widehat{Z}_{n}(m)=n[vec(R_{12}(m))]^{T}[R_{22}^{-1}(0)\otimes R_{11}^{-1}(0)][vec(R_{12}(m))]$, $R_{ij}(m)=D[(\widehat{r}_{ii}(0))^{-1/2}]$ $\widehat{r}_{ij}(m)D[(\widehat{r}_{jj}(0))^{-1/2}]$, $\widehat{r}_{ij}(m)$ is the sample cross-covariance matrix between $\{\widehat{\eta}_{it}\}$ and $\{\widehat{\eta}_{jt+m}\}$, $\widetilde{Z}_{n}(m)$ is defined in the same way as $\widehat{Z}_{n}(m)$ with $\widehat{\eta}_{st}$ being replaced by $\widetilde{\eta}_{st}$, $\widetilde{\eta}_{st}$ is the residual from a fitted VAR($p$) model for $Y_{st}$, $\overline{K}(\cdot)$ is a kernel function, $h$ stands for the bandwidth, $A_{1}=\int_{-\infty}^{\infty}[\overline{K}(z)]^{2}dz$, $B_{1}=\int_{-\infty}^{\infty}[\overline{K}(z)]^{4}dz$, and

Note that $G_{1n}$ is for testing the cross-correlation between $\eta_{1t}$ and $\eta_{2t}$, and $G_{2n}$ is its modified version for small $n$; $W_{1n}$ is towards the same goal as $G_{1n}$ but with ability to detect the cross-correlation beyond lag $M$, and $W_{2n}$ is the modified version of $W_{1n}$. Under certain conditions, the limiting null distribution of $G_{1n}$ or $G_{2n}$ is $\chi^{2}_{(2M+1)d_{1}d_{2}}$, and that of $W_{1n}$ or $W_{2n}$ is $N(0,1)$.

In all simulation studies, we set $m=0$ and $3$ for the single HSIC-based tests $S_{sn}(m)$, and set $M=3$ and $6$ for the joint HSIC-based test $J_{sn}(M)$. Because $S_{1n}(0)=S_{2n}(0)$, the results of $S_{2n}(0)$ are absent. For $G_{sn}(M)$, we choose $M=3, 6$ and $9$. For $W_{sn}(h)$, we follow Hong (1996) to choose $p=3$ (or 6) when $n=100$ (or 200), and use the kernel function $\overline{K}(z)=\sin(\pi z)/(\pi z)$ (Daniel kernel) with the bandwidth $h=h_{1}, h_{2}$ or $h_{3}$, where $h_{1}=[\log(n)], h_{2}=[3n^{0.2}]$, and $h_{3}=[3n^{0.3}]$. The significance level $\alpha$ is set to be $1\%,5\%$ and 10%.

Table (ref) reports the power of all tests for model ((ref)), and the sizes of all tests are corresponding to those in EGP 1. From this table, our findings are as follows:

(i) The sizes of all single HSIC-based tests $S_{sn}$ are close to their nominal ones in most cases, while the sizes of other tests are a little unsatisfactory. For instance, $J_{sn}$ are slightly oversized especially at $\alpha=5\%$ and $10\%$, while $W_{1n}$ (or $W_{2n}$) is slightly oversized (or undersized) when $n=200$ (or 100) at all levels. The size performance of $G_{sn}$ depends on $M$: a larger value of $M$ leads to a more undersized behavior especially at $\alpha=10\%$, although $G_{2n}$ in general has a better performance than $G_{1n}$.

(ii) In all examined cases, the single HSIC-based test $S_{1n}(0)$ is much more powerful than other tests in EGPs 2-3 and 5-6, and the single HSIC-based test $S_{2n}(3)$ has a significant power advantage in EGP 4. These results are expected, since $S_{1n}(0)$ and $S_{2n}(3)$ are tailored to examine the dependence at specific lags $0$ and $3$, respectively, which are the set-ups of our EGPs.

(iii) For the linear dependence case (i.e., EGP 2), the joint HSIC-based tests $J_{sn}$ have a comparable power performance as $G_{sn}$, and they are much less powerful than $W_{1n}(h_{1})$ but much more powerful than $W_{2n}(h_{3})$ when $n=100$. For the non-linear dependence case (i.e., EGPs 3-6), the joint HSIC-based tests $J_{sn}$ in general are much more powerful than the tests $G_{sn}$ and $W_{sn}$ especially when $n=200$. The only exception is $J_{1n}$ in EGP 4, since $J_{1n}$ can not detect the dependence between $\eta_{1t+m}$ and $\eta_{2t}$ at lag $m=3$. In contrast, $J_{2n}$ performs very well here.

(iv) In all examined cases, the power of $J_{sn}$ and $G_{sn}$ decreases as the value of $M$ increase, while this tendency is vague for $W_{sn}$.

Overall, our single HSIC-based tests are very powerful in detecting dependence at specific lags, and our joint HSIC-based tests exhibit a significant power advantage in detecting non-linear dependence, which can not be easily examined by other tests.

Conditional variance models

We generate 1000 replications of sample size $n$ from the following two conditional variance models:

where $\{\eta_{1t}\}$ and $\{\eta_{2t}\}$ are two sequences of i.i.d. random vectors generated as for model ((ref)). Model ((ref)) contains two CC-MGARCH models studied in Tse (2002). For each replication, we fit the above models by using the Gaussian-QMLE method. Denote by $\{\widehat{\eta}_{1t}\}$ and $\{\widehat{\eta}_{2t}\}$ the residuals from the fitted models. Based on $\{\widehat{\eta}_{1t}\}$ and $\{\widehat{\eta}_{2t}\}$, we compute $S_{sn}(m)$ and $J_{sn}(M)$, and their critical values as for model ((ref)).

At the same time, we also compute the test statistics $L_{sn}(M)$ and $T_{sn}(M)$ ($L_{sn}$ and $T_{sn}$ in short) in Tchahou and Duchesne (2013), where

Here, $\rho_{\widehat{q}_{1t},\widehat{q}_{2t}}(m)$ is the sample cross-correlation between $\{\widehat{q}_{1t}\}$ and $\{\widehat{q}_{2t+m}\}$, $C_{ij}(m)$ is the sample cross-covariance matrix between $\{\widehat{\varphi}_{it}\}$ and $\{\widehat{\varphi}_{jt+m}\}$, $\widehat{q}_{st}=\widehat{\eta}_{st}^{T}\widehat{\eta}_{st}$, and $\widehat{\varphi}_{st}=vech(\widehat{\eta}_{st}\widehat{\eta}_{st}^{T})$. It is worth noting that $L_{1n}$ (or $T_{1n}$) is for testing the cross-correlation between two transformed (or original) residuals, and $L_{2n}$ (or $T_{2n}$) is its modified version for small $n$. Under certain conditions, the limiting null distribution of $L_{1n}$ or $L_{2n}$ is $\chi^{2}_{(2M+1)}$, and that of $T_{1n}$ or $T_{2n}$ is $\chi^{2}_{(2M+1)d_{1}^{*}d_{2}^{*}}$, where $d_{s}^{*}=d_{s}(d_{s}+1)/2$ for $s=1, 2$.

In all simulation studies, we choose the values of $m$ and $M$ as for model ((ref)). The significance level $\alpha$ is set to be $1\%,5\%$ and 10%. Table (ref) summarizes the power results of all tests for model ((ref)), and the sizes of all tests are corresponding to those in EGP 1. From this table, our findings are as follows:

(i) The sizes of all tests are close to their nominal ones, although most of $T_{sn}$ are slightly oversized.

(ii) Similar to the results in model ((ref)), the single HSIC-based test $S_{1n}(0)$ or $S_{1n}(3)$ as expected is the most powerful one among all tests.

(iii) For the linear dependence case (i.e., EGP 2), all joint HSIC-based tests $J_{sn}$ are much more powerful than $L_{sn}$ and $T_{sn}$. For the non-linear dependence case (i.e., EGP 3-6), all $J_{sn}$ still have larger power than $L_{sn}$ and $T_{sn}$ in most cases, but this advantage is small especially for $J_{sn}(6)$. There are two exceptions that some $J_{sn}$ exhibit low power: first, $J_{1n}(3)$ and $J_{1n}(6)$ as argued for model ((ref)) have no power in EGP 4; second, $J_{2n}(6)$ is less powerful than most of $L_{sn}$ and $T_{sn}$ especially for $n=200$. Since the cross-correlation between $\eta_{1t}^2$ and $\eta_{2t}^2$ is high in EGPs 2-6, the relative good power performance of $L_{sn}$ and $T_{sn}$ in some cases is not out of our expectation.

(iv) For the tests $J_{sn}$, $L_{sn}$ and $T_{sn}$, their power decreases as the value of $M$ increases in all examined cases.

Overall, our single HSIC-based tests as usual have good power in detecting dependence at specific lags, and our joint HSIC-based tests could be more powerful than other tests in detecting either linear or non-linear dependence.

A real example

In this section, we study two bivariate time series. The first bivariate time series consist of two index series from the Russian market and the Indian market: the Russia Trading System Index (RTSI) and the Bombay Stock Exchange Sensitive Index (BSESI). The second bivariate time series include two Chinese indexes: the ShangHai Securities Composite index (SHSCI) and the ShenZhen Index (SZI). The data are observed on a daily basis (from Monday to Friday), beginning on 8 October 2014, and ending on 29 September 2017. In all there were 1088 days, missing data due to holidays are removed before the analysis, and hence the final data set include $n=672$ daily observations. The resulting four time series are denoted by \{RTSI$_t$; $t=1,\ldots,n\}$, \{BSESI$_t$; $t=1,\ldots,n\}$, \{SHSCI$_t$; $t=1,\ldots,n\}$ and \{SZI$_t$; $t=1,\ldots,n\}$, respectively.

As usual, we consider the log-return of each data set:

An investigation on the ACF and PACF of $Y_{1t,1}, Y_{1t,2},Y_{2t,1},Y_{2t,2}$ and their squares indicates that they do not have a conditional mean structure but a conditional variance structure. Motivated by this, we use the following BEKK model with Gaussian-QMLE method to fit $Y_{1t}$ and $Y_{2t}$:

for $s=1, 2$, where $A_s=C_{s0}^{T}C_{s0}$ with $C_{s0}$ being a triangular $2\times 2$ matrix, and $B_{s1},\cdots,B_{sp}, C_{s1},\ldots,C_{sq}$ are all $2\times 2$ diagonal matrixes. Table (ref) reports the estimates for both fitted models. The p-values of portmanteau tests $Q(3)$, $Q(6)$ and $Q(9)$ in Ling and Li (1997) are $0.7698,0.5179, 0.5967$ for $Y_{1t}$ and $0.5048, 0.7328, 0.8746$ for $Y_{2t}$. This implies that both fitted BEKK models are adequate.

Next, we apply our joint HSIC-based tests $J_{sn}(M)$ to check whether $Y_{1t}$ and $Y_{2t}$ behave independently of each other. As a comparison, we also consider the tests $L_{sn}(M)$ and $T_{sn}(M)$ for the testing purpose. Table (ref) reports the $p$-value for all six tests. From Table (ref), we find that except for $J_{2n}(M)$ with $M\geq7$, all examined joint HSIC-based tests $J_{sn}(M)$ convey strong evidence that $Y_{1t}$ and $Y_{2t}$ are not independent. However, neither $L_{sn}(M)$ nor $T_{sn}(M)$ is able to do this for $M\geq 2$.

To get more information, we further plot the values of the single version of $J_{sn}$, $L_{1n}$ and $T_{1n}$ in Fig 1. That is, Fig 1 plots the values of $S_{sn}(m)$, $L_{1n, s}(m)$, and $T_{1n, s}(m)$ for $m\geq 0$, where

and all notations are inherited from Section 5.2. The limiting null distribution of $L_{1n,s}(m)$ is $\chi^2_1$, and that of $T_{1n,s}(m)$ is $\chi_9^{2}$. Similar to $S_{sn}(m)$, $L_{1n,s}(m)$ and $T_{1n,s}(m)$ capture the linear dependence between $\eta_{1t}$ and $\eta_{1t+m}$ at the specific lag $m$. The corresponding single version results for $L_{2n}$ and $T_{2n}$ are similar to those for $L_{1n}$ and $T_{1n}$, and hence they are not displayed here.

From Fig 1, we first find that all single tests indicate a strong contemporaneously causal relationship between the Chinese market and the Russian and Indian (R&I) market. Second, $S_{1n}(1)$ implies that the R&I market has significant influence on the Chinese market one day later, while according to $S_{2n}(3)$ (or $S_{2n}(10)$), the impact of the Chinese market to the R&I market appears after three (or ten) days. These findings demonstrate an asymmetric causal relationship between two markets. Since none of examined $L_{1n,s}(m)$ and $T_{1n,s}(m)$ can detect a causal relationship for $m\geq1$, the contemporaneous causal relationship mainly results in the significance of $L_{sn}(1)$ and $T_{sn}(1)$ in Table 4, and the lagged causal relationship is possible to be non-linear. As the R&I market has a higher degree of globalization and marketization, it could have a quicker impact to other economies. On the contrary, the Chinese market is more localized, and its influence to other economies tends to be slower but can last for a longer term. This long-term effect may be caused by “the Belt and Road Initiatives” program raised by Chinese government since 2015. Hence, the asymmetric phenomenon between two markets seems reasonable, and it may help the government to make more efficient policy and the investors to design more useful investment strategies.

Concluding remarks

In this paper, we apply the HSIC principle to derive some novel one-sided omnibus tests for detecting independence between two multivariate stationary time series. The resulting HSIC-based tests have asymptotical Gaussian representation under the null hypothesis, and they are shown to be consistent. A residual bootstrap method is used to obtain the critical values for our HSIC-based tests, and its validity is justified. Unlike the existing cross-correlation-based tests for linear dependence, our HSIC-based tests look for the general dependence between two un-observable innovation vectors, and hence they can give investigators more complete information on the causal relationship between two time series. The importance of our HSIC-based tests is illustrated by simulation results and real data analysis. Due to the generality of the HSIC method, the methodology developed in this paper may be applied to many other important testing problems such as testing for model adequacy (Davis et al. 2016), testing for independence among multi-dynamic systems (Pfister et al. 2017), or testing for independence in high dimensional systems (Yao et al. 2017). We leave these interesting topics as potential future study.