Testing Serial Independence of Object-Valued Time Series

TL;DR

This paper introduces a fully nonparametric, general method for testing serial independence in object-valued time series within metric spaces, utilizing distance covariance and spectral analysis.

Contribution

It extends distance covariance to object-valued time series and develops a new spectral density-based test with bootstrap calibration, applicable beyond Euclidean spaces.

Findings

Effective in detecting nonlinear dependence

Applicable to complex, non-Euclidean data

Validated through simulations and real data

Abstract

We propose a novel method for testing serial independence of object-valued time series in metric spaces, which is more general than Euclidean or Hilbert spaces. The proposed method is fully nonparametric, free of tuning parameters, and can capture all nonlinear pairwise dependence. The key concept used in this paper is the distance covariance in metric spaces, which is extended to auto distance covariance for object-valued time series. Furthermore, we propose a generalized spectral density function to account for pairwise dependence at all lags and construct a Cramer-von Mises type test statistic. New theoretical arguments are developed to establish the asymptotic behavior of the test statistic. A wild bootstrap is also introduced to obtain the critical values of the non-pivotal limiting null distribution. Extensive numerical simulations and two real data applications are conducted to illustrate the effectiveness and versatility of our proposed method.

Tables

Table 1. Table 1 : Empirical rejection rate of functional time series when n = 200 𝑛 200 n=200

DGP		Proposed				GK			HHHR				Zhang
		CvM-P	CvM-B	KS-P	KS-B	$h=1$	$h=3$	$h=5$	$h=5$	$h=10$	$h=30$	$h=50$	$b=1$	$b=5$	$b=10$	$b^{\ast }$
BM	F	0.059	0.057	0.059	0.060	0.048	0.052	0.040	0.068	0.064	0.064	0.071	0.053	0.057	0.062	0.062
	B	0.055	0.058	0.057	0.060	0.047	0.050	0.041	0.068	0.064	0.063	0.071	0.055	0.054	0.060	0.062
BB	F	0.052	0.051	0.050	0.050	0.040	0.044	0.034	0.057	0.054	0.068	0.072	0.043	0.038	0.036	0.043
	B	0.051	0.054	0.047	0.051	0.039	0.042	0.032	0.057	0.054	0.068	0.072	0.046	0.038	0.037	0.044
FARCH(0.5)	F	0.113	0.085	0.119	0.081	0.127	0.111	0.087	0.106	0.099	0.086	0.086	0.052	0.046	0.051	0.056
	B	0.114	0.090	0.114	0.081	0.127	0.111	0.087	0.106	0.100	0.086	0.087	0.047	0.043	0.050	0.053
FARCH(0.9)	F	0.296	0.194	0.282	0.173	0.288	0.279	0.240	0.213	0.164	0.107	0.090	0.045	0.037	0.041	0.045
	B	0.294	0.184	0.275	0.174	0.287	0.278	0.241	0.213	0.163	0.108	0.090	0.043	0.037	0.038	0.046
FNMA-BM	F	0.839	0.626	0.810	0.585	0.373	0.222	0.172	0.180	0.149	0.102	0.084	0.039	0.039	0.045	0.046
	B	0.841	0.624	0.804	0.591	0.372	0.225	0.171	0.178	0.148	0.101	0.084	0.042	0.040	0.045	0.045
FNMA-BB	F	0.774	0.518	0.730	0.496	0.512	0.269	0.188	0.211	0.164	0.113	0.110	0.044	0.028	0.026	0.035
	B	0.774	0.518	0.730	0.494	0.492	0.262	0.185	0.208	0.164	0.110	0.108	0.042	0.027	0.028	0.033
FAR-G-BM	F	0.984	0.984	0.981	0.981	0.964	0.832	0.726	0.921	0.822	0.662	0.559	0.973	0.970	0.958	0.973
	B	0.985	0.984	0.981	0.980	0.961	0.828	0.728	0.921	0.818	0.656	0.557	0.972	0.968	0.957	0.972
FAR-G-BB	F	0.984	0.985	0.978	0.981	0.972	0.827	0.692	0.866	0.741	0.540	0.453	0.957	0.942	0.924	0.944
	B	0.985	0.983	0.979	0.980	0.966	0.818	0.682	0.866	0.741	0.541	0.454	0.960	0.938	0.932	0.943
FAR-W-BM	F	0.973	0.972	0.971	0.970	0.919	0.786	0.678	0.918	0.864	0.718	0.636	0.984	0.979	0.968	0.983
	B	0.974	0.973	0.966	0.971	0.920	0.788	0.680	0.920	0.865	0.719	0.635	0.984	0.974	0.975	0.980
FAR-W-BB	F	0.968	0.969	0.966	0.964	0.933	0.706	0.550	0.851	0.715	0.520	0.441	0.950	0.914	0.905	0.931
	B	0.974	0.971	0.967	0.971	0.936	0.712	0.562	0.857	0.719	0.521	0.444	0.953	0.931	0.913	0.939

Table 2. Table 2 : Empirical rejection rate of covariance matrix time series generated by CAW p ​ ( ρ ) subscript CAW 𝑝 𝜌 \text{CAW}_{p}(\rho) model

Method			$p=2$			$p=5$			$p=8$
			$\rho =0$	$\rho =0.5$	$\rho =0.75$	$\rho =0$	$\rho =0.5$	$\rho =0.75$	$\rho =0$	$\rho =0.5$	$\rho =0.75$
Proposed	Euc	CvM-P	0.043	0.419	0.997	0.067	0.570	1.000	0.042	0.745	1.000
		CvM-B	0.044	0.412	0.999	0.069	0.558	1.000	0.046	0.744	1.000
		KS-P	0.045	0.403	0.997	0.061	0.529	1.000	0.046	0.681	1.000
		KS-B	0.048	0.394	0.997	0.072	0.517	1.000	0.053	0.694	1.000
	log-Euc	CvM-P	0.051	0.347	0.997	0.059	0.309	1.000	0.061	0.256	1.000
		CvM-B	0.049	0.349	0.996	0.058	0.316	1.000	0.061	0.254	1.000
		KS-P	0.048	0.340	0.993	0.060	0.292	1.000	0.063	0.226	1.000
		KS-B	0.051	0.339	0.994	0.065	0.294	1.000	0.068	0.223	1.000
	Chol	CvM-P	0.052	0.387	0.995	0.045	0.340	1.000	0.059	0.348	1.000
		CvM-B	0.052	0.378	0.996	0.045	0.339	1.000	0.065	0.356	1.000
		KS-P	0.046	0.355	0.995	0.049	0.314	1.000	0.061	0.305	1.000
		KS-B	0.050	0.356	0.995	0.051	0.321	1.000	0.065	0.315	1.000
	Riemann	CvM-P	0.053	0.346	0.996	0.062	0.316	1.000	0.058	0.290	1.000
		CvM-B	0.053	0.345	0.995	0.061	0.314	1.000	0.062	0.287	1.000
		KS-P	0.051	0.342	0.995	0.068	0.293	1.000	0.058	0.262	1.000
		KS-B	0.053	0.340	0.994	0.068	0.298	1.000	0.070	0.259	1.000
	Vech-Euc	CvM-P	0.048	0.399	0.994	0.059	0.527	1.000	0.044	0.703	1.000
		CvM-B	0.049	0.389	0.994	0.062	0.517	1.000	0.042	0.697	1.000
		KS-P	0.052	0.376	0.992	0.058	0.500	1.000	0.046	0.635	1.000
		KS-B	0.049	0.374	0.993	0.061	0.491	1.000	0.047	0.639	1.000
mADCV	TC	$\lambda =1/10$	0.034	0.146	0.899	0	0	0.898	0	0	0.666
		$\lambda =2/10$	0.036	0.131	0.871	0	0.001	0.931	0	0	0.895
		$\lambda =3/10$	0.037	0.119	0.822	0	0.003	0.963	0	0	0.983
	DAN	$\lambda =1/10$	0.029	0.253	0.984	0	0.002	0.980	0	0	0.843
		$\lambda =2/10$	0.029	0.209	0.971	0	0.003	0.973	0	0	0.915
		$\lambda =3/10$	0.034	0.179	0.939	0	0.003	0.979	0	0	0.974
	PAR	$\lambda =1/10$	0.026	0.243	0.983	0	0.003	0.977	0	0	0.837
		$\lambda =2/10$	0.026	0.190	0.964	0	0.002	0.969	0	0	0.885
		$\lambda =3/10$	0.037	0.175	0.932	0	0.003	0.977	0	0	0.974
	BAR	$\lambda =1/10$	0.026	0.257	0.985	0	0.005	0.979	0	0	0.744
		$\lambda =2/10$	0.026	0.213	0.976	0	0.004	0.981	0	0	0.890
		$\lambda =3/10$	0.030	0.200	0.958	0	0.003	0.981	0	0	0.953

Table 3. Table 3 : Empirical rejection rate of univariate distributions when n = 200 𝑛 200 n=200

			W1		W2		KS		KL		IS		LS
			Permt	Boot	Permt	Boot	Permt	Boot	Permt	Boot	Permt	Boot	Permt	Boot
ATM(0)		CvM	0.049	0.052	0.049	0.050	0.050	0.050	0.050	0.050	0.047	0.045	0.050	0.049
		KS	0.049	0.052	0.050	0.051	0.050	0.057	0.049	0.050	0.049	0.050	0.053	0.049
ATM(1)	$\rho =0.2$	CvM	0.186	0.183	0.178	0.180	0.178	0.173	0.239	0.234	0.272	0.262	0.203	0.201
		KS	0.172	0.176	0.175	0.171	0.164	0.164	0.224	0.228	0.264	0.253	0.188	0.186
	$\rho =0.5$	CvM	0.624	0.623	0.622	0.624	0.603	0.595	0.745	0.732	0.814	0.789	0.696	0.689
		KS	0.607	0.604	0.599	0.603	0.579	0.575	0.728	0.716	0.800	0.782	0.678	0.676
	$\rho =1.0$	CvM	1.000	1.000	1.000	1.000	1.000	1.000	1.000	1.000	0.970	0.978	1.000	1.000
		KS	1.000	1.000	1.000	1.000	1.000	1.000	1.000	1.000	0.970	0.977	1.000	1.000
ATM(4)	$\rho =0.5$	CvM	0.176	0.180	0.178	0.176	0.173	0.171	0.240	0.245	0.279	0.265	0.186	0.186
		KS	0.184	0.188	0.188	0.186	0.175	0.180	0.253	0.253	0.289	0.280	0.202	0.198
	$\rho =0.8$	CvM	0.730	0.732	0.732	0.737	0.717	0.714	0.846	0.834	0.862	0.824	0.762	0.759
		KS	0.735	0.736	0.742	0.745	0.724	0.718	0.855	0.840	0.870	0.838	0.756	0.750
	$\rho =1.0$	CvM	0.996	0.996	0.997	0.997	0.995	0.995	0.998	0.998	0.996	0.988	0.998	0.997
		KS	0.994	0.996	0.996	0.996	0.995	0.996	0.998	0.999	0.997	0.989	0.998	0.997

Table 4. Table 4 : P-values for white noise testing in CIDRs of IBM and MCD in 2007.

Test	CvM	KS	GK1	GK3	GK5	HHHR5	HHHR10	HHHR30	HHHR50	Zhang1	Zhang5	Zhang10
IBM	0.008	0.010	0.003	0.008	0.004	0.000	0.000	0.002	0.100	0.040	0.010	0.027
MCD	0.858	0.938	0.416	0.044	0.098	0.785	0.207	0.424	0.153	0.552	0.301	0.344

Table 5. Table 5 : Empirical rejection rate of univariate time series when n = 200 𝑛 200 n=200

Method	IID				NMA(2)				ARCH(2)				TAR(2)
Proposed	CvM-P	0.052	CvM-B	0.054	CvM-P	1.000	CvM-B	0.997	CvM-P	0.765	CvM-B	0.582	CvM-P	0.995	CvM-B	0.994
	KS-P	0.055	KS-B	0.056	KS-P	1.000	KS-B	0.995	KS-P	0.734	KS-B	0.564	KS-P	0.993	KS-B	0.991
BP	BP-1	0.045	BP-9	0.045	BP-1	0.416	BP-9	0.272	BP-1	0.336	BP-9	0.358	BP-1	0.051	BP-9	0.053
	BP-3	0.047	BP-12	0.039	BP-3	0.374	BP-12	0.227	BP-3	0.426	BP-12	0.327	BP-3	0.055	BP-12	0.056
	BP-6	0.047	BP-15	0.038	BP-6	0.297	BP-15	0.204	BP-6	0.402	BP-15	0.293	BP-6	0.050	BP-15	0.045
LB	LB-1	0.049	LB-9	0.051	LB-1	0.419	LB-9	0.281	LB-1	0.339	LB-9	0.374	LB-1	0.051	LB-9	0.059
	LB-3	0.050	LB-12	0.050	LB-3	0.386	LB-12	0.246	LB-3	0.435	LB-12	0.343	LB-3	0.062	LB-12	0.072
	LB-6	0.050	LB-15	0.051	LB-6	0.308	LB-15	0.229	LB-6	0.415	LB-15	0.318	LB-6	0.055	LB-15	0.065
H	H-D3	0.068	H-P3	0.067	H-D3	0.463	H-P3	0.475	H-D3	0.411	H-P3	0.438	H-D3	0.074	H-P3	0.078
	H-D6	0.067	H-P6	0.069	H-D6	0.453	H-P6	0.452	H-D6	0.476	H-P6	0.494	H-D6	0.077	H-P6	0.078
	H-D9	0.067	H-P9	0.066	H-D9	0.432	H-P9	0.414	H-D9	0.492	H-P9	0.491	H-D9	0.078	H-P9	0.080
	H-D12	0.066	H-P12	0.069	H-D12	0.391	H-P12	0.374	H-D12	0.475	H-P12	0.465	H-D12	0.074	H-P12	0.074
	H-D15	0.071	H-P15	0.069	H-D15	0.358	H-P15	0.350	H-D15	0.455	H-P15	0.452	H-D15	0.077	H-P15	0.076
	H-Q3	0.067	H-B3	0.068	H-Q3	0.470	H-B3	0.468	H-Q3	0.410	H-B3	0.401	H-Q3	0.076	H-B3	0.077
	H-Q6	0.066	H-B6	0.070	H-Q6	0.456	H-B6	0.467	H-Q6	0.483	H-B6	0.471	H-Q6	0.078	H-B6	0.078
	H-Q9	0.068	H-B9	0.070	H-Q9	0.433	H-B9	0.448	H-Q9	0.493	H-B9	0.498	H-Q9	0.075	H-B9	0.078
	H-Q12	0.069	H-B12	0.069	H-Q12	0.388	H-B12	0.420	H-Q12	0.478	H-B12	0.492	H-Q12	0.076	H-B12	0.079
	H-Q15	0.073	H-B15	0.070	H-Q15	0.364	H-B15	0.389	H-Q15	0.455	H-B15	0.472	H-Q15	0.080	H-B15	0.079
FP	FP-D3	0.052	FP-P3	0.054	FP-D3	1.000	FP-P3	1.000	FP-D3	0.899	FP-P3	0.905	FP-D3	0.998	FP-P3	0.998
	FP-D9	0.042	FP-P9	0.042	FP-D9	0.966	FP-P9	0.970	FP-D9	0.793	FP-P9	0.794	FP-D9	0.938	FP-P9	0.931
	FP-D15	0.035	FP-P15	0.030	FP-D15	0.729	FP-P15	0.790	FP-D15	0.618	FP-P15	0.654	FP-D15	0.698	FP-P15	0.752
	FP-D21	0.036	FP-P21	0.030	FP-D21	0.347	FP-P21	0.521	FP-D21	0.428	FP-P21	0.518	FP-D21	0.422	FP-P21	0.564
	FP-D25	0.034	FP-P25	0.030	FP-D25	0.278	FP-P25	0.390	FP-D25	0.396	FP-P25	0.443	FP-D25	0.355	FP-P25	0.446
	FP-Q3	0.054	FP-B3	0.055	FP-Q3	1.000	FP-B3	1.000	FP-Q3	0.904	FP-B3	0.904	FP-Q3	0.999	FP-B3	0.999
	FP-Q9	0.046	FP-B9	0.046	FP-Q9	0.981	FP-B9	0.991	FP-Q9	0.816	FP-B9	0.839	FP-Q9	0.946	FP-B9	0.972
	FP-Q15	0.037	FP-B15	0.036	FP-Q15	0.850	FP-B15	0.935	FP-Q15	0.690	FP-B15	0.753	FP-Q15	0.809	FP-B15	0.906
	FP-Q21	0.033	FP-B21	0.029	FP-Q21	0.606	FP-B21	0.806	FP-Q21	0.557	FP-B21	0.653	FP-Q21	0.618	FP-B21	0.772
	FP-Q25	0.033	FP-B25	0.032	FP-Q25	0.474	FP-B25	0.690	FP-Q25	0.482	FP-B25	0.578	FP-Q25	0.515	FP-B25	0.682
EL	EL25-1.2	0.163	EL-50-1.2	0.209	EL25-1.2	0.127	EL50-1.2	0.170	EL25-1.2	0.111	EL50-1.2	0.163	EL25-1.2	0.180	EL50-1.2	0.240
	EL25-1.5	0.120	EL50-1.5	0.152	EL25-1.5	0.080	EL50-1.5	0.104	EL25-1.5	0.072	EL50-1.5	0.108	EL25-1.5	0.124	EL50-1.5	0.168
	EL25-1.8	0.090	EL50-1.8	0.113	EL25-1.8	0.070	EL50-1.8	0.078	EL25-1.8	0.054	EL50-1.8	0.066	EL25-1.8	0.093	EL50-1.8	0.119
	EL25-2.1	0.078	EL50-2.1	0.090	EL25-2.1	0.066	EL50-2.1	0.068	EL25-2.1	0.045	EL50-2.1	0.051	EL25-2.1	0.074	EL50-2.1	0.088
	EL25-2.4	0.070	EL50-2.4	0.076	EL25-2.4	0.064	EL50-2.4	0.064	EL25-2.4	0.043	EL50-2.4	0.046	EL25-2.4	0.066	EL50-2.4	0.073
	EL25-2.7	0.069	EL50-2.7	0.071	EL25-2.7	0.063	EL50-2.7	0.063	EL25-2.7	0.041	EL50-2.7	0.042	EL25-2.7	0.063	EL50-2.7	0.067
	EL25-3.0	0.068	EL50-3.0	0.069	EL25-3.0	0.063	EL50-3.0	0.063	EL25-3.0	0.040	EL50-3.0	0.041	EL25-3.0	0.063	EL50-3.0	0.065

Table 6. Table 6 : Empirical rejection rate of bivariate time series when n = 200 𝑛 200 n=200

Method			IID			NMA(2)			ARCH(2)			MAR(2)
			$\rho =0$	$\rho =0.4$	$\rho =0.7$	$\rho =0$	$\rho =0.4$	$\rho =0.7$	$\rho =0$	$\rho =0.4$	$\rho =0.7$	$\rho =0$	$\rho =0.4$	$\rho =0.7$
Proposed	CvM	Permt	0.053	0.052	0.056	0.997	0.999	1.000	0.348	0.346	0.348	1.000	1.000	1.000
		Boot	0.057	0.053	0.056	0.978	0.989	1.000	0.261	0.240	0.242	1.000	1.000	1.000
	KS	Permt	0.051	0.053	0.056	0.996	0.997	1.000	0.342	0.319	0.335	1.000	1.000	1.000
		Boot	0.051	0.058	0.060	0.971	0.984	0.997	0.224	0.230	0.234	1.000	1.000	1.000
mLB		$h=1$	0.115	0.034	0.016	0.522	0.470	0.395	0.139	0.066	0.042	0.157	0.020	0.012
		$h=3$	0.172	0.025	0.004	0.522	0.468	0.338	0.167	0.046	0.035	0.160	0.020	0.005
		$h=6$	0.194	0.009	0.001	0.467	0.417	0.270	0.130	0.031	0.017	0.150	0.008	0.000
		$h=9$	0.213	0.007	0.001	0.442	0.403	0.229	0.106	0.021	0.009	0.146	0.004	0.000
mADCV	TC	$\lambda =1/10$	0.047	0.038	0.038	0.779	0.756	0.682	0.242	0.256	0.259	0.997	1.000	1.000
		$\lambda =2/10$	0.038	0.044	0.054	0.732	0.714	0.637	0.244	0.250	0.255	0.991	0.999	1.000
		$\lambda =3/10$	0.056	0.049	0.043	0.672	0.647	0.594	0.223	0.239	0.236	0.972	0.995	0.999
	DAN	$\lambda =1/10$	0.047	0.045	0.048	0.869	0.846	0.762	0.233	0.237	0.250	1.000	1.000	1.000
		$\lambda =2/10$	0.056	0.041	0.049	0.873	0.833	0.768	0.271	0.268	0.279	1.000	1.000	1.000
		$\lambda =3/10$	0.047	0.046	0.051	0.844	0.811	0.742	0.272	0.285	0.291	0.999	1.000	1.000
	PAR	$\lambda =1/10$	0.051	0.043	0.042	0.866	0.831	0.755	0.239	0.256	0.245	1.000	1.000	1.000
		$\lambda =2/10$	0.050	0.038	0.045	0.855	0.822	0.744	0.269	0.274	0.291	0.999	1.000	1.000
		$\lambda =3/10$	0.055	0.045	0.044	0.825	0.800	0.721	0.269	0.280	0.279	0.999	1.000	1.000
	BAR	$\lambda =1/10$	0.051	0.044	0.049	0.871	0.839	0.763	0.220	0.225	0.226	1.000	1.000	1.000
		$\lambda =2/10$	0.052	0.039	0.037	0.869	0.839	0.758	0.255	0.262	0.274	1.000	1.000	1.000
		$\lambda =3/10$	0.053	0.046	0.049	0.859	0.830	0.751	0.274	0.278	0.288	0.999	1.000	1.000
mADCF	TC	$\lambda =1/10$	0.046	0.043	0.045	0.782	0.746	0.677	0.232	0.260	0.258	0.997	1.000	1.000
		$\lambda =2/10$	0.045	0.042	0.053	0.710	0.694	0.626	0.226	0.256	0.250	0.987	1.000	0.999
		$\lambda =3/10$	0.064	0.046	0.046	0.639	0.638	0.577	0.222	0.235	0.236	0.964	0.993	0.999
	DAN	$\lambda =1/10$	0.053	0.042	0.045	0.873	0.845	0.749	0.226	0.217	0.239	1.000	1.000	1.000
		$\lambda =2/10$	0.055	0.036	0.046	0.863	0.823	0.744	0.265	0.260	0.266	1.000	1.000	1.000
		$\lambda =3/10$	0.058	0.044	0.048	0.832	0.816	0.727	0.264	0.279	0.288	0.998	1.000	1.000
	PAR	$\lambda =1/10$	0.052	0.037	0.044	0.862	0.841	0.760	0.235	0.233	0.250	1.000	1.000	1.000
		$\lambda =2/10$	0.058	0.041	0.043	0.860	0.830	0.750	0.257	0.266	0.282	0.999	1.000	1.000
		$\lambda =3/10$	0.057	0.045	0.053	0.825	0.800	0.715	0.262	0.277	0.286	0.998	1.000	1.000
	BAR	$\lambda =1/10$	0.052	0.046	0.051	0.868	0.839	0.753	0.217	0.211	0.224	1.000	1.000	1.000
		$\lambda =2/10$	0.053	0.042	0.049	0.868	0.837	0.757	0.252	0.249	0.273	1.000	1.000	1.000
		$\lambda =3/10$	0.054	0.043	0.050	0.860	0.829	0.736	0.268	0.272	0.278	0.999	1.000	1.000

Table 7. Table 7 : Empirical rejection rate of multivariate time series with different n 𝑛 n and p 𝑝 p

Method			$n=100$					$n=250$
			$p=2$	$p=5$	$p=10$	$p=20$	$p=30$	$p=2$	$p=5$	$p=10$	$p=20$	$p=30$
Proposed	CvM	Permt	0.488	0.606	0.749	0.899	0.970	0.929	0.994	0.998	1.000	1.000
		Boot	0.485	0.604	0.755	0.898	0.970	0.932	0.992	0.998	1.000	1.000
	KS	Permt	0.465	0.572	0.705	0.863	0.968	0.911	0.988	0.998	0.999	1.000
		Boot	0.465	0.579	0.708	0.862	0.966	0.912	0.990	0.998	0.999	1.000
mLB		$h=1$	0.616	0.729	0.773	0.835	0.872	0.930	0.966	0.977	0.970	0.978
		$h=3$	0.523	0.674	0.726	0.823	0.881	0.858	0.926	0.955	0.953	0.966
		$h=6$	0.460	0.629	0.694	0.828	0.893	0.782	0.875	0.920	0.931	0.940
		$h=9$	0.435	0.601	0.689	0.824	0.906	0.737	0.846	0.890	0.908	0.930
mADCV	TC	$\lambda =1/10$	0.214	0.105	0.010	0.000	0.000	0.620	0.639	0.462	0.088	0.007
		$\lambda =2/10$	0.192	0.105	0.025	0.000	0.000	0.493	0.550	0.428	0.148	0.038
		$\lambda =3/10$	0.191	0.128	0.035	0.002	0.000	0.433	0.497	0.462	0.274	0.137
	DAN	$\lambda =1/10$	0.360	0.190	0.024	0.000	0.000	0.858	0.912	0.809	0.331	0.037
		$\lambda =2/10$	0.305	0.171	0.031	0.000	0.000	0.760	0.814	0.722	0.343	0.095
		$\lambda =3/10$	0.269	0.186	0.047	0.001	0.000	0.666	0.730	0.657	0.358	0.134
	PAR	$\lambda =1/10$	0.346	0.180	0.021	0.000	0.000	0.833	0.893	0.768	0.270	0.024
		$\lambda =2/10$	0.282	0.158	0.027	0.000	0.000	0.739	0.786	0.667	0.242	0.051
		$\lambda =3/10$	0.257	0.172	0.041	0.000	0.000	0.636	0.711	0.629	0.322	0.119
	BAR	$\lambda =1/10$	0.362	0.191	0.023	0.000	0.000	0.865	0.923	0.825	0.336	0.039
		$\lambda =2/10$	0.318	0.176	0.029	0.000	0.000	0.794	0.859	0.766	0.319	0.061
		$\lambda =3/10$	0.296	0.182	0.042	0.000	0.000	0.728	0.794	0.707	0.370	0.118
mADCF	TC	$\lambda =1/10$	0.208	0.106	0.012	0.000	0.000	0.614	0.630	0.439	0.074	0.007
		$\lambda =2/10$	0.190	0.104	0.021	0.000	0.000	0.477	0.533	0.410	0.129	0.033
		$\lambda =3/10$	0.195	0.125	0.034	0.002	0.000	0.422	0.494	0.483	0.294	0.164
	DAN	$\lambda =1/10$	0.362	0.186	0.026	0.000	0.000	0.852	0.915	0.819	0.358	0.054
		$\lambda =2/10$	0.293	0.163	0.030	0.000	0.000	0.763	0.820	0.721	0.337	0.099
		$\lambda =3/10$	0.264	0.155	0.030	0.001	0.000	0.657	0.730	0.633	0.328	0.136
	PAR	$\lambda =1/10$	0.338	0.174	0.021	0.000	0.000	0.832	0.886	0.772	0.269	0.033
		$\lambda =2/10$	0.280	0.152	0.024	0.000	0.000	0.736	0.791	0.672	0.285	0.069
		$\lambda =3/10$	0.252	0.155	0.034	0.001	0.000	0.623	0.705	0.618	0.332	0.132
	BAR	$\lambda =1/10$	0.357	0.177	0.014	0.000	0.000	0.867	0.923	0.815	0.304	0.031
		$\lambda =2/10$	0.317	0.167	0.023	0.000	0.000	0.794	0.847	0.730	0.274	0.054
		$\lambda =3/10$	0.283	0.160	0.028	0.001	0.000	0.724	0.774	0.666	0.326	0.102

Equations

$$\mathrm{dcov}(X,Y)$$

$$=$$

$$\mathrm{dcov}(X,Y)$$

$$=$$

$$M_p(\Omega )=\{\nu \in M(\Omega ):\text{ for some }\omega \in \Omega ,{\int }_{\Omega }{d}^p(\omega ,x)\mathrm{d}\nu (x)<\infty \}.$$

$$\begin{array}{rl}& d^{(1)}(x)={\int }_{\Omega }d(x,x^{\prime })\mathrm{d}\nu (x^{\prime }),D={\int }_{\Omega }{d}^{(1)}(x)\mathrm{d}\nu (x),\\ & d_{\nu }(x,x^{\prime })=d(x,x^{\prime })-d^{(1)}(x)-d^{(1)}(x^{\prime })+D.\end{array}$$

$$\mathrm{dcov}(X,Y)=\mathbb{E}[d_{{\nu }_X}(X,X^{\prime })d_{{\nu }_Y}(Y,Y^{\prime })],$$

$${\iint }_{{\Omega }^2}d(x_1,x_2)\mathrm{d}{\nu }_{-}(x_1)\mathrm{d}{\nu }_{-}(x_2)<0.$$

$$V(k)=\mathbb{E}[d_{\nu }(X_t,X_t^{\prime })d_{\nu }(X_{t-k},X_{t-k}^{\prime })],$$

$$H_0:\{X_t{\}}_{t=1}^n\text{mbox}isi.i.d.\text{mbox}v.s.H_a:V(k)\ne 0\text{mbox}forsomek\ne 0.$$

$${\tilde{a}}_{ij}(k)=\{\begin{array}{ll}{a}_{ij}-\frac{{\sum }_{t=k+1}^n{a}_{it}}{n-k-2}-\frac{{\sum }_{t^{\prime }=k+1}^n{a}_{t^{\prime }j}}{n-k-2}+\frac{{\sum }_{k+1\le t\ne t^{\prime }\le n}{a}_{t{t}^{\prime }}}{(n-k-1)(n-k-2)},& i\ne j\\ 0,& i=j.\end{array}$$

$$V_n(k)=\frac{1}{(n-k)(n-k-3)}\sum _{i,j=k+1}^n{\tilde{a}}_{ij}(k){\tilde{b}}_{ij}(k).$$

$$V_n(k)={\left(\genfrac{}{}{0ex}{}{n-k}{4}\right)}^{-1}\sum _{k+1\le i<j<q<r\le n}h(Z_i^{(k)},Z_j^{(k)},Z_q^{(k)},Z_r^{(k)}),$$

$$h\left(Z_i,Z_j,Z_q,Z_r\right)=\frac{1}{24}\sum _{(i_1,i_2,i_3,i_4)}^{(i,j,q,r)}d(X_{i_1},X_{i_2})\left[d(Y_{i_3},Y_{i_4})+d(Y_{i_1},Y_{i_2})-2d(Y_{i_1},Y_{i_3})\right],$$

$${\bar{a}}_{i.}(k)=\frac{1}{n-k}\sum _{t=1}^n{a}_{it},{\bar{a}}_{.j}(k)=\frac{1}{n-k}\sum _{t=k+1}^n{a}_{tj},{\bar{a}}_{..}(k)=\frac{1}{(n-k{)}^2}\sum _{t,t^{\prime }=k+1}^n{a}_{t,t^{\prime }},$$

$${\left\{(n-k)V_n(k)\right\}}_{k=1}^K{\to }_d{\left\{{\xi }_k\right\}}_{k=1}^K,$$

$${\lambda }_{\ell }{e}_{\ell }(z)=\mathbb{E}[d_{\nu }(x,X)d_{\nu }(y,Y)]e_{\ell }(z),$$

$$f(\zeta )=(2\pi {)}^{-1}\sum _{k=-\infty }^{\infty }V(k)e^{-ik\zeta },\zeta \in [-\pi ,\pi ],$$

$$S_n(\zeta )=\sum _{k=1}^{n-4}(n-k)V_n(k){\Psi }_k(\zeta ),\zeta \in [0,\pi ],$$

$$\text{textsc}{CvM}_n={\int }_0^{\pi }{S}_n^2(\zeta )\mathrm{d}\zeta .$$

$$\text{textsc}{KS}_n=\underset{\zeta \in [0,\pi ]}{\sup }|S_n(\zeta )|.$$

$$\{S_n(\zeta ){\}}_{\zeta \in [0,\pi ]}\Rightarrow \{S(\zeta ){\}}_{\zeta \in [0,\pi ]},\text{in }{L}_2[0,\pi ],$$

$$\text{textsc}{CvM}_n{\to }_d{\int }_0^{\pi }{S}^2(\zeta )\mathrm{d}\zeta .$$

$$V_n^{\ast }(k{)}^{(b)}=\frac{1}{(n-k)(n-k-3)}\sum _{k+1\le i\ne j\le n}{w}_i^{(b)}(k){\tilde{a}}_{ij}{\tilde{b}}_{ij}{w}_j^{(b)}(k),$$

$$S_n^{\ast }(\zeta )=\sum _{k=1}^{n-4}(n-k)V_n^{\ast }(k{)}^{(b)}{\Psi }_k(\zeta ),\zeta \in [0,\pi ],$$

$$T_n^{\ast }{\to }_{d^{\ast }}T,\text{mbox}inprobability.$$

$$\{(n-k)V_n^{\ast }(k){\}}_{k=1}^K{\to }_{d\ast }{\left\{{\xi }_k\right\}}_{k=1}^K\text{mbox}inprobability,$$

$$\text{textsc}{CvM}_n^{\ast }{\to }_{d^{\ast }}{\int }_0^{\pi }{S}^2(\zeta )\mathrm{d}\zeta ,\text{mbox}inprobability,$$

$$\mathbb{P}(\text{textsc}{CvM}_n>\text{textsc}{CvM}_n^{\ast }(1-\alpha ))\to 1.$$

$${\beta}\odot T(x)=\left\{\begin{array}[]{ll}x+\beta(T(x)-x),&0<\beta\leq 1\\ x,&\beta=0\\ x+\beta(x-T^{-1}(x)),&-1\leq\beta<0\end{array}\right.$$

$$X_t(s)=100\{\log [P_t(s)/P_t(0)]\},s\in [0,1],$$

Full text

Keywords: Distance covariance; Non-Euclidean valued data; Random object; Spectral test; White noise testing.

1 Introduction

Random objects in general metric spaces have become increasingly common in modern statistical and econometric research. For example, the intraday return path of a financial asset (Aue et al.,, 2017), the annual composition of energy sources (Zhu and Müller, 2023b, ), social networks (Board and Meyer-ter Vehn,, 2021), EEG scans or MRI fiber tracts of patients (Kurtek et al.,, 2012) can all be viewed as random objects in certain metric spaces, although they are typically given specific names as functional data, compositional data, network data, image data and curve data, among others. The concept of random objects also undoubtedly include classical notions such as vectors, covariance matrices and distributions. Instead of building specific models for each of them, by viewing these objects as random elements in metric spaces, we may be able to simplify the modeling and inference while preserving the ability of extracting meaningful information and patterns in a unified fashion (Petersen and Müller,, 2019; Dubey and Müller,, 2019, 2020; Zhang et al.,, 2022).

For many endeavors in this area, the data they analyzed is collected with a natural ordering, i.e., the data is object-valued time series. However, most existing modeling and inference techniques either presume temporal independence or construct time series models without conducting diagnostic checking to assess the goodness-of-fit. This is mainly due to the unavailability of appropriate tests. As a result, researchers often overlook the serial dependence in their data and fail to account for its impact on their analyses. Consequently, the validity and reliability of their findings could be compromised. This motivates us to develop a new test for serial independence of object-valued time series.

Testing serial independence has a long and rich history in statistics and econometrics, with a vast literature that cannot be exhaustively listed. The early work dates back to Box and Pierce, (1970) and Ljung and Box, (1978). Since then, numerous tests have been proposed for univariate random variables or multivariate vectors in Euclidean space, with much attention devoted to testing for second-order uncorrelatedness, see Li and McLeod, (1981), Deo, (2000), Lobato, (2001), Escanciano and Lobato, (2009), Shao, (2011), to name a few. These tests typically capture linear serial dependence in data, and have no power against nonlinear dependence. Nonlinear serial dependence is indeed prevalent among many real-world time series, and many parametric nonlinear models have been proposed to capture nonlinear dependence. A prominent example is GARCH model, which implies uncorrelatedness but is serially dependent. Note that several tests have been developed to target at higher order dependence (Li and Mak,, 1994; Ling and Li,, 1997) and at general nonlinear dependence (Hong,, 1999; Escanciano and Velasco,, 2006). We also note the recent developments for testing white noise hypothesis in functional time series in Hilbert space, see e.g. Gabrys and Kokoszka, (2007), Horváth et al., (2013) and Zhang, (2016).

Despite many tests available for testing the serial independence/uncorrelatedness in Euclidean and Hilbert spaces, they cannot be directly used for testing serial independence of object-valued time series, because of the lack of classic algebraic operations in general metric space, such as addition, multiplication and taking inner product. To fill this gap, we propose to build a new test based on distance covariance, which was originally proposed by Székely et al., (2007) to measure dependence among two random vectors in Euclidean space using characteristic functions, and later extended by Zhou, (2012) into time series setting by using auto-distance covariance (ADCV). A crucial feature of the (auto-)distance covariance is that it can capture both linear and nonlinear serial dependence in data. Although many researchers have employed this idea into serial dependence testing, they are only valid in Euclidean space (Fokianos and Pitsillou,, 2017, 2018; Davis et al.,, 2018).

Based on the concept of distance covariance in metric space by Lyons, (2013), we are able to extend ADCV in Euclidean space into general metric spaces, which then naturally serves the goal of testing independence at fixed lags. To take into account pairwise (in)dependence at all lags, we then propose a generalized spectral density/distribution function using ADCVs in the same spirit of classical spectral density/distribution function using autocovariances. Following the developments in Shao, (2011), who proposed a spectrum-based test for white noise hypothesis of a univariate time series, we develop a Cramér von–Mises (CvM) type test statistic and study the limiting behavior of the test under both the null and alternatives. Our test significantly enhances classical spectrum-based tests in two fundamental ways. First, it captures both linear and nonlinear serial dependencies, a capability lacking in conventional methods except for Hong, (1999) and Escanciano and Velasco, (2006). Second, its versatility extends to data objects in general metric spaces, which is much broader than the Euclidean or Hilbert space considered in the literature. To the best of our knowledge, this is the first formal attempt at testing temporal independence for object-valued time series in a unified manner. Unlike conventional autocovariance function or spectrum based tests for Euclidean time series, the estimation of ADCVs is based on U-statistics. Therefore, new theoretical arguments are developed to establish the asymptotic theory for the proposed test statistic. Since the limiting null distribution is nonstandard and non-pivotal, we propose a wild bootstrap approach to facilitate practical implementation of our test. The bootstrap consistency is also established.

We now introduce the notation. Let $(\Omega,d)$ be a separable metric space, and let $(\mathbb{S},\mathbb{P},\mathcal{F})$ be the probability space. For a vector $a\in\mathbb{R}^{q}$, we denote its Euclidean norm as $|a|_{q}$. Denote $``\to_{p}"$ and $``\to_{d}"$ the convergence in probability and in distribution, respectively. Denote $L_{2}[0,\pi]$ the Hilbert space $\mathbb{H}$ of all square integrable functions on $[0,\pi]$ (with respect to Lebesgue measure) with inner product $\langle f,g\rangle=\int_{[0,\pi]}f(\zeta)g^{c}(\zeta)\mathrm{d}(\zeta)$ where $g^{c}(\zeta)$ denotes the complex conjugate of $g(\zeta)$, and the norm $\|f\|=\langle f,f\rangle^{1/2}$. We denote $``\Rightarrow"$ as weak convergence in $L_{2}[0,\pi]$.

The rest of the paper is organized as follows. We first provide backgrounds of (auto-) distance covariance in metric spaces in Section 2. Section 3 then introduces our distance covariance based test statistics, and investigates their asymptotic distributions under null and alternatives. Section 4 provides the wild bootstrap algorithm for approximating the limiting null distribution. Extensive numerical experiments are conducted in Section 5 with competing methods for testing serial independence of the functional time series in Hilbert space, the covariance matrix time series, and the univariate distributional time series. Section 6 illustrates the usefulness and versatility of our tests via two meaningful real data applications in financial data and human mortality data. Section 7 concludes. Additional numerical results and all the technical proofs are provided in the supplement. The code is available at https://github.com/hjgao117/JiangGaoShao.

2 Preliminaries

In this section, we provide some background on the concept of distance covariance in metric spaces and its use in quantifying dependence. The extension to auto-distance covariance (ADCV) for the time series setting is also introduced.

2.1 Distance covariance

The concept of distance covariance was first introduced by Székely et al., (2007) as a measure of dependence between two random vectors $X\in\mathbb{R}^{p}$ and $Y\in\mathbb{R}^{q}$. It is defined as the weighted integral of the discrepancy between the joint characteristic function and the product of the marginal characteristic functions of $(X,Y)$.

Definition 2.1** (Distance Covariance in Euclidean Space).**

For $X\in\mathbb{R}^{p}$ and $Y\in\mathbb{R}^{q}$, the distance covariance is given by

[TABLE]

where $c_{d}=\pi^{(1+d)/2}/\Gamma((1+d)/2)$ and $\Gamma$ is the Gamma function.

Note that for notational simplicity, we shall use $\mathrm{dcov}$ instead of $\mathrm{dcov}^{2}$ as used in the original definition of distance covariance. Székely et al., (2007) also managed to derive the following alternative definition when $\mathbb{E}[|X|_{p}^{2}+|Y|_{q}^{2}]<\infty$,

[TABLE]

with $(X^{\prime},Y^{\prime})$ and $(X^{\prime\prime},Y^{\prime\prime})$ being independent copies of $(X,Y)$.

Distance covariance can be used to characterize the dependence between $X$ and $Y$ due to the following crucial property.

Proposition 2.1**.**

$\mathrm{dcov}(X,Y)\geq 0$, and the equality holds iff $X$ is independent of $Y$.

2.1 has proven to be a very powerful tool in Euclidean space, with its use widely appeared in mutual dependence testing (Yao et al.,, 2018), feature screening (Li et al.,, 2012), dimension reduction (Sheng and Yin,, 2016), among many other applications. However, 2.1 is not easily extended to general metric spaces because conventional algebraic manipulation such as addition and multiplication may not be applied. Instead, by viewing $|\cdot|_{q}$ as a metric in Euclidean space, the alternative definition in (1) is extendable.

2.2 Distance covariance in metric space

Let $(\Omega,d)$ be a separable metric space. Let $M(\Omega)$ be the set of probability measures on $\Omega$, we denote the subset of $M(\Omega)$ possessing $p$th moment as

[TABLE]

For $\nu\in M_{1}(\Omega)$, define (Lyons,, 2013, Lemma  2.1)

[TABLE]

Clearly, for two independent random objects $X,X^{\prime}\sim\nu(\cdot)$, we can write $d^{(1)}(X)=\mathbb{E}[d(X,X^{\prime})|X]$, and $D=\mathbb{E}[d(X,X^{\prime})]$.

Note that if $X$ and $Y$ are two random vectors taking values in conventional $q$-dimensional Euclidean space (i.e., $(\Omega,d)=(\mathbb{R}^{q},|\cdot|_{q})$) with marginal distributions being $\nu_{X}$ and $\nu_{Y}$ respectively, then (1) is equivalent to $\mathrm{dcov}(X,Y)=\mathbb{E}[d_{\nu_{X}}(X,X^{\prime})d_{\nu_{Y}}(Y,Y^{\prime})]$. Therefore, Lyons, (2013) proposed the following definition of distance covariance in general metric spaces.

Definition 2.2** (Distance Covariance in Metric Space).**

For $X,Y$ taking values in $(\Omega,d)$ whose marginals are $\nu_{X},\nu_{Y}\in M_{1}(\Omega)$ respectively, the distance covariance between $X$ and $Y$ is given as

[TABLE]

where $(X^{\prime},Y^{\prime})$ is an independent copy of $(X,Y)$.

However, as pointed out by Lyons, (2013), in general metric spaces, the above definition alone is insufficient for 2.1 to hold, and additional topological assumption is required. He then introduced the concept of strong negative type to resolve this issue.

Definition 2.3** (Strong Negative Type).**

We say $(\Omega,d)$ is of strong negative type, if for $\nu_{1},\nu_{2}\in M_{1}(\Omega)$ such that $\nu_{1}\neq\nu_{2}$, and $\nu_{-}=\nu_{1}-\nu_{2}$, we have

[TABLE]

Note that the class of metric spaces of strong negative type is actually quite large, for example, every separable Hilbert space is of strong negative type. However, we also note there are many spaces that do not satisfy 2.3, e.g., $\mathbb{R}^{q}$ with $L_{p}$-metric for $3\leq q\leq\infty$, $2<p\leq\infty$ are not of strong negative type. We refer to Lyons, (2013, 2014, 2020) for more discussions and examples. With the notion of strong negative type, distance covariance then completely characterizes the (in)dependence in metric space, given in the following proposition.

Proposition 2.2**.**

[Theorem 3.11 in Lyons, (2013)] If $(\Omega,d)$ is of strong negative type, for $X,Y$ taking values in $(\Omega,d)$ whose marginals are $\nu_{X},\nu_{Y}\in M_{1}(\Omega)$ respectively, 2.1 continues to hold.

2.3 Auto-distance covariance

The distance covariance, as defined in 2.2, measures the dependence between two random objects and might not be readily applicable to the time series context. Therefore, in order to address this limitation, we introduce the concept of auto-distance covariance (ADCV). This notion was first proposed by Zhou, (2012) for conventional Euclidean valued time series by measuring temporal dependence between $\{X_{t}\}_{t\in\mathbb{Z}}\in\mathbb{R}^{q}$ and its lagged observation $\{X_{t-k}\}_{t\in\mathbb{Z}}$ for a fixed lag order $k\in\mathbb{Z}$. Here we generalize the idea for time series objects in metric space.

Definition 2.4** (ADCV).**

Assume that $\{X_{t}\}_{t\in\mathbb{Z}}$ is a sequence of strict stationary time series taking values in $(\Omega,d)$. For $k\in\mathbb{Z}$, we call

[TABLE]

the auto-distance covariance of $X_{t}$ at lag $k$, where $\{X_{t}^{\prime}\}_{t\in\mathbb{Z}}$ is an independent copy of $\{X_{t}\}_{t\in\mathbb{Z}}$.

It is clear that $V(k)=V(-k)$ for $k<0$, and $V(0)=\mathbb{E}[d_{\nu}^{2}(X_{t},X_{t}^{\prime})]>0$. In addition, by 2.2, $V(k)=0$, $k\neq 0$, iff $X_{t}$ and $X_{t-k}$ are independent of each other. We will exploit this property to build our test statistics given in the next section.

Remark 2.1**.**

In Zhou, (2012), the ADCV for Euclidean valued time series is based on 2.1 by replacing $(X,Y)$ with $(X_{t},X_{t-k})$. However, for general metric spaces, such treatment is no longer valid. We thus call for 2.2. It can be shown that for $(\Omega,d)=(\mathbb{R}^{q},|\cdot|_{q})$, 2.4 is equivalent to the one in Zhou, (2012).

3 Test statistics and asymptotics

Given a sequence of stationary random objects $\{X_{t}\}_{t=1}^{n}$ that reside in a separable metric space $(\Omega,d)$ of strong negative type, we are interested in testing the serial independence of $\{X_{t}\}_{t=1}^{n}$. The hypothesis testing problem is formulated as

[TABLE]

3.1 Test statistics at fixed lags

To illustrate the idea, we first consider a relatively simpler task by forming a test based on ADCV at fixed lags $1\leq k\leq K$. This suggests that we find an empirical estimator for $V(k)$. Intuitively, under $H_{0}$, $(X_{t^{\prime}},X_{t^{\prime}-k})$ naturally forms an independent copy of $(X_{t},X_{t-k})$ if $t^{\prime}\neq t$, which motivates us to estimate $V(k)$ by replacing terms involving (conditional) expectations in $d_{\nu}(X_{t},X_{t-k})$ with their empirical counterparts.

Motivated by the estimator in Euclidean space (Székely et al.,, 2007), we construct the empirical estimator $V_{n}(k)$ by adopting the $\mathcal{U}$-centering approach in Székely and Rizzo, (2014). Specifically, let $a_{ij}(k)=d(X_{i},X_{j})$ for $k+1\leq i,j\leq n$, we denote $\{\widetilde{a}_{ij}(k)\}_{i,j=k+1,\cdots,n}$ as its $\mathcal{U}$-centered version:

[TABLE]

Define $\tilde{b}_{ij}(k)$ similarly for $b_{ij}(k)=d(X_{i-k},X_{j-k})$ with $k+1\leq i,j\leq n$. Intuitively, for large $n$, $\tilde{a}_{ij}(k)$ (or $\tilde{b}_{ij}(k)$) approximates the value of $d_{\nu}(X_{i},X_{j})$ (or $d_{\nu}(X_{i-k},X_{j-k})$).

We then estimate $V(k)$ by

[TABLE]

The $\mathcal{U}$-centering approach was originally proposed by Székely and Rizzo, (2014) to provide an unbiased estimator of distance covariance. It is adopted here for simplifying technical analysis as we can alternatively rewrite $V_{n}(k)$ by the following fourth order U-statistic (Zhang et al.,, 2018),

[TABLE]

where $Z_{i}^{(k)}=(X_{i},X_{i-k})$, and the kernel function is given by

[TABLE]

such that $Z_{\ell}=(X_{\ell},Y_{\ell})$, $\ell\in\{i,j,q,r\}$, and $\sum_{(i_{1},i_{2},i_{3},i_{4})}^{(i,j,q,r)}$ denotes the summation over all 24($=4!$) permutations of the 4-tuple of indices $(i,j,q,r)$.

Remark 3.1**.**

Alternatively, one could use the conventional empirical estimator of $V(k)$ given by $\widehat{V}_{n}(k)=\frac{1}{(n-k)^{2}}\sum_{i,j=k+1}^{n}\hat{a}_{ij}(k)\hat{b}_{ij}(k),$ where $\hat{a}_{ij}(k)=a_{ij}(k)-\bar{a}_{i.}(k)-\bar{a}_{.j}(k)+\bar{a}_{..}(k)$,

[TABLE]

with $\hat{b}_{ij}(k)$ similarly defined. Note that $\widehat{V}_{n}(k)$ is a combination of $V$-statistics, the technical analysis could be more involved than the $\mathcal{U}$-centering approach (Székely and Rizzo,, 2014).

To analyze the asymptotic behaviors of $V_{n}(k)$, we make the following moment assumption.

Assumption 3.1**.**

The marginal distribution $\nu(\cdot)\in M_{4}(\Omega)$.

The following theorem derives the limiting distribution of (5) by exploiting the properties of U-statistics (e.g. Lee, (1990)). We emphasize here that even under $H_{0}$, $V_{n}(k)$ is not a conventional U-statistic with kernels applied to i.i.d. samples. For example, $Z_{t}^{(k)}=(X_{t},X_{t-k})$ and $Z_{t+k}^{(k)}=(X_{t+k},X_{t})$ are not independent of each other. Therefore, more delicate technical treatments are required.

Theorem 3.1**.**

Under $H_{0}$, and suppose Assumption 3.1 holds. Then, fix $K\geq 1$, as $n\to\infty$,

[TABLE]

where $\xi_{k}=_{d}\sum_{\ell=1}^{\infty}\lambda_{\ell}\{[G_{\ell}^{(k)}]^{2}-1\}$. Here $\{G_{\ell}^{(k)}\}_{\ell=1,2,\cdots,\infty;k=1,\cdots,K}$ is a sequence of i.i.d. standard Gaussian random variables, and $\{\lambda_{\ell}\}_{\ell=1}^{\infty}$ and $\{e_{\ell}(\cdot)\}_{\ell=1}^{\infty}$ are sequences of nonzero eigenvalues and orthonormal eigenfunctions corresponding to

[TABLE]

where $z=(x,y)\in\Omega^{2}$ and $Z=(X,Y)\in\Omega^{2}$ with $X$ and $Y$ being independent copies of $X_{t}$.

3.1 improves Zhou, (2012) under $H_{0}$ not only by extending the result in Euclidean space to more general metric spaces but also by providing joint convergence of sample ADCVs. In addition, 3.1 is crucial for proving the asymptotics of spectrum based test below. In practice, note that $\{\xi_{k}\}_{k=1}^{K}$ is a sequence of centered mixture of i.i.d. $\chi^{2}(1)$ random variables, which are non-pivotal. Below we adopt a wild bootstrap method to approximate the limiting null distributions, and details are deferred to Section 4.

3.2 Generalized spectral test

For our testing purpose, it is natural to combine ADCVs at all lags. We propose the following generalized spectral density,

[TABLE]

and generalized spectral distribution function, $F(\zeta)=\int_{0}^{\zeta}f(z)\mathrm{d}z,~{}\zeta\in[0,\pi].$ Clearly, $f(\cdot)$ (or $F(\cdot)$) is motivated by the spectral density (or distribution) function in Euclidean space and Hilbert space, where the classical spectral density function is defined by replacing $V(k)$ in (8) by $\gamma(k)$, and $\gamma(k)=\mathrm{Cov}(X_{t},X_{t-k})$ is the auto-covariance at lag $k$.

Our generalized spectral density (or function) serves as a “generalization” of the classical spectral density (or function) in two significant aspects. First, we aim to test serial independence, which is a stronger notion than serial uncorrelatedness. In the conventional spectral density function, $\gamma(k)$ is the auto-covariance measuring linear correlation between $X_{t}$ and its lag-$k$ observation $X_{t-k}$, which implies that spectrum based tests can only capture linear serial dependence in the second-order structure. By contrast, if we consider $V(k)$, the auto-distance covariance at lag $k$, we can additionally measure the lag $k$ nonlinear dependence. Second, the definition of auto-covariance requires certain algebraic operation such as subtraction and inner product in the conventional Euclidean or Hilbert space, whereas the auto-distance covariance does not, making it applicable to more general metric spaces.

Under $H_{0}$, it holds that $F(\zeta)=V(0)\Psi_{0}(\zeta),$ where $\Psi_{k}(\zeta)=\sin(k\zeta)/(k\pi)\mathbf{1}(k\neq 0)+\zeta/(2\pi)\mathbf{1}(k=0)$, so that we can construct a test by comparing the empirical estimator of $F(\zeta)-V(0)\Psi_{0}(\zeta)$ and zero. In particular, we define the process

[TABLE]

and consider the following Cramér von–Mises (CvM) type statistic

[TABLE]

Similar CvM type statistics have been adopted in Euclidean space and Hilbert space for testing second-order white noise or martingale differences, see, e.g. Deo, (2000), Escanciano and Velasco, (2006), Shao, (2011), Zhang, (2016) among many others.

Remark 3.2**.**

One could also consider KS (Kolmogorov–Smirnov) type statistics of the form

[TABLE]

However, the “$\sup$” functional is no longer a continuous map in $L_{2}$ Hilbert space, and therefore the process convergence result below and the continuous mapping theorem do not lead to the asymptotic null distribution of KS test statistic, and a rigorous theoretical investigation is beyond the scope of this paper. Nevertheless, our numerical studies in Section 5 suggest that $\textsc{KS}_{n}$ also delivers satisfactory performance.

Remark 3.3**.**

*We briefly discuss the difference of our tests with the distance covariance based tests by Fokianos and Pitsillou, (2017, 2018) in Euclidean space. First, their test is designed for Euclidean valued time series, while ours has much more broad applications in other metric spaces. Second, their test is founded on the basis of Hong, (1999) by noticing the connection between the characteristic function based definition of distance covariance (c.f. 2.1) and the dependence metric defined in Hong, (1999) in terms of the weighting function. However their approach is no longer valid in metric space. Third, their test statistic is based on a kernel weighted sum of empirical ADCVs at a set of lags. Therefore, they require an additional tuning parameter to truncate the kernel, which is typically known to affect the convergence rate and hence power. On the contrary, our test is tuning-free, and is more powerful even in the Euclidean setting, as demonstrated in Appendix A.1 in the supplement. *

Our next theorem obtains the weak convergence result for $S_{n}(\zeta)$ in the sense of $L_{2}$ metric.

Theorem 3.2**.**

Under $H_{0}$, suppose Assumption 3.1 holds, then

[TABLE]

where $\{S(\zeta)\}_{\zeta\in[0,\pi]}=_{d}\left\{\sum_{k=1}^{\infty}\xi_{k}\Psi_{k}(\zeta)\right\}_{\zeta\in[0,\pi]}$ with $\{\xi_{k}\}$ defined in 3.1.

Based on 3.2 and note that the integral functional is continuous, the following corollary is a natural consequence of continuous mapping theorem.

Corollary 3.1**.**

Under $H_{0}$, suppose Assumption 3.1 holds, then

[TABLE]

Therefore, given significance level $\alpha\in(0,1)$, we reject $H_{0}$ if $\textsc{CvM}_{n}\geq\textsc{CvM}(1-\alpha),$ where $\textsc{CvM}(1-\alpha)$ denotes the $1-\alpha$ quantile of $\int_{0}^{\pi}S^{2}(\zeta)\mathrm{d}\zeta$. For practical implementation, we need to invoke the bootstrap method in Section 4 to approximate the above critical value.

To study the behaviors of the tests under $H_{a}$, we impose the following $\beta$-mixing (absolutely regularity) conditions on $\{X_{t}\}_{t\in\mathbb{Z}}$. Basically, it requires weak temporal dependence of the process and similar assumptions are adopted in Fokianos and Pitsillou, (2017, 2018) and Davis et al., (2018) for distance covariance based inference in time series. Note that weak temporal dependence is also required in Zhou, (2012) using physical dependence measures (Wu,, 2005).

Assumption 3.2**.**

$\{X_{t}\}_{t=1}^{n}$* is a stationary $\beta$-mixing process, with marginal distribution $\nu(\cdot)\in M_{2+2\delta}(\Omega)$ for some $\delta>0$. Furthermore, the mixing coefficient satisfies that $\beta(n)=O(n^{-r})$ for some $r>0$.*

3.3 below derives the asymptotics of the test statistics under $H_{a}$.

Theorem 3.3**.**

Suppose Assumption 3.2 hold, we have $V_{n}(k)\to_{a.s.}V(k)$ for $k=1,2,\cdots$. Under $H_{a}$, as $n\to\infty$, $n^{-1}\textsc{CvM}_{n}\to_{p}\sum_{k=1}^{\infty}V^{2}(k)\int_{0}^{\pi}\Psi_{k}^{2}(\zeta)\mathrm{d}\zeta.$

Clearly, under $H_{a},$ we have $V(k)>0$ for some $k\neq 0$, which implies that our test has asymptotically power 1 under $H_{a}$.

4 Bootstrap

In this section, we use wild bootstrap method to approximate the limiting null distribution of $\textsc{CvM}_{n}$, and show its asymptotic validity.

Algorithm 4.1**.**

Let $\mathbf{w}=(w_{1},\cdots,w_{n})$ be a sequence of independent random variables following the Rademacher distribution: $\mathbb{P}(w_{i}=1)=\mathbb{P}(w_{i}=-1)=\frac{1}{2}.$

For each $k=1,\cdots,n-4$, define the doubly weighted bootstrapped version of $V_{n}$ by

[TABLE]

*where $\tilde{a}_{ij}$ and $\tilde{b}_{ij}$ are $\mathcal{U}$ centering version of $a_{ij}=d(X_{i},X_{j})$ and $b_{ij}=d(X_{i-k},X_{j-k})$ respectively, and $\mathbf{w}^{(b)}(k)$ is an independent random realization of $\mathbf{w}$ with $w_{i}^{(b)}(k)$ being the *$i$th element of $\mathbf{w}^{(b)}(k)$. 2. 2.

Obtain the bootstrapped process

[TABLE]

and compute $\textsc{CvM}_{n}^{*(b)}=\int_{0}^{\pi}[S_{n}^{*(b)}]^{2}(\zeta)\mathrm{d}\zeta.$ 3. 3.

Repeat steps 1-2 $B$ times and collect $\{\textsc{CvM}_{n}^{*(b)}\}_{b=1}^{B}$. 4. 4.

*Denote the *$\alpha$th upper percentile of $\{\textsc{CvM}_{n}^{*(b)}\}_{b=1}^{B}$ by $\textsc{CvM}_{n}^{*}(1-\alpha)$, and reject $H_{0}$ if $\textsc{CvM}_{n}>\textsc{CvM}_{n}^{*}(1-\alpha)$.

For each fixed $k$, the bootstrapped auto-distance covariance estimator $V_{n}^{*}(k)^{(b)}$ can be viewed as a realization of the independent wild bootstrap for degenerate U-statistics, see Dehling and Mikosch, (1994) and Lee et al., (2020). They are further combined to obtain the bootstrap realization of $\textsc{CvM}_{n}^{*(b)}$. Note here we use independent copies of wild bootstrap weights $\mathbf{w}(k)$ for different $k$’s. This is implied by 3.1 that despite common components in constructing $\{V_{n}(k)\}_{k=1}^{K}$, their asymptotic distributions $\{\xi_{k}\}_{k=1}^{K}$ are independent of each other under $H_{0}$.

Remark 4.1**.**

For our test targeting at the global null of serial independence, independent wild bootstrap is sufficient. If one is interested in making statistical inference for $V_{n}(k)$ alone, potential serial dependence at other lag orders should be accounted for. A possible extension along this direction would be to pinpoint the minimum lag at which dependencies exist by sequentially testing the nullity of $V_{n}(k)$. To achieve this, the above wild bootstrap method should be adjusted accordingly. For example, Leucht and Neumann, (2013) used the dependent wild bootstrap in Shao, (2010) by constructing temporarily dependent weights $\{w_{i}\}_{i=1}^{n}$, as an extension of wild bootstrap in Dehling and Mikosch, (1994) for degenerate U-statistics. See also Zhou, (2012) for a subsampling based procedure.

Remark 4.2**.**

Another approach to approximate the limiting null distribution is by using permutation. Although both bootstrap and permutation methods can be based on pre-calculated pairwise distance metrics $d(X_{i},X_{j})$ for $1\leq i<j\leq n$, permutation is less computationally efficient. This is because in the $\mathcal{U}$-centering step, permutation method needs to re-calculate $\tilde{a}_{ij}$ (or $\tilde{b}_{ij}$) in (4), which involves the summation of $d(X_{\tau(i)},X_{\tau(j)})$ with $\tau(i)$ being the permutation order of $i$th original observation. In contrast, bootstrap method does not require such re-calculation, because $\tilde{a}_{ij}$ (or $\tilde{b}_{ij}$) is kept constant, as shown in (10), across all bootstrap statistics.

We then proceed to analyze the asymptotic behavior of the bootstrapped test statistics. We first introduce the notion of convergence in distribution in probability, see Definition 2 in Li et al., (2003).

Definition 4.1** (convergence in distribution in probability).**

For a sequence of bootstrapped statistics $T_{n}^{*}$ which depends on the random sample $\{X_{t}\}_{t=1}^{n}$, we say that $T_{n}^{*}|(X_{1},X_{2},\cdots)$ converges to $T|(X_{1},X_{2},\cdots)$ in distribution in probability if for any subsequence $T_{n}^{\prime}$, there exists a further subsequence $T_{n}^{\prime\prime}$ such that $T_{n}^{\prime\prime}|(X_{1},X_{2},\cdots)$ converges to $T|(X_{1},X_{2},\cdots)$ in distribution for almost every sequence $(X_{1},X_{2},\cdots)$. The following notation is used,

[TABLE]

Theorem 4.1**.**

Under $H_{0}$, suppose Assumptions 3.1 holds. For $\mathbf{w}$ such that $\mathbb{E}w_{i}=0$, $\mathbb{E}w_{i}^{2}=1$ and $\mathbb{E}w_{i}^{4}<\infty$. Then, for any fixed $K\geq 1$, as $n\to\infty$,

[TABLE]

where $\xi_{k}$ is defined in Theorem 3.1.

Theorem 4.2**.**

Under the same conditions of 4.1,

[TABLE]

where $\{S(\zeta)\}_{\zeta\in[0,\pi]}=_{d}\left\{\sum_{k=1}^{\infty}\xi_{k}\Psi_{k}(\zeta)\right\}_{\zeta\in[0,\pi]}$, with $\{\xi_{k}\}$ defined in 3.1. Furthermore, we have $\mathbb{P}(\textsc{CvM}_{n}>\textsc{CvM}_{n}^{*}(1-\alpha))\to\alpha.$

Therefore, the wild bootstrap in 4.1 provides a consistent approximation of the limiting null distribution of our test statistics. When $H_{0}$ is violated, we show that the bootstrapped tests have consistent powers.

Theorem 4.3**.**

Suppose 3.2 with $\delta=1$, then under $H_{a}$, for $\alpha\in(0,1)$, we have

[TABLE]

5 Simulation studies

In this section, we look into the finite sample performance of the proposed method on several different types of data objects. Specifically, the performance on the functional time series in $L_{2}[0,1]$ is investigated in Section 5.1 whereas that on a sequence of covariance matrices can be found in Section 5.2. Finally, Section 5.3 reports the simulation results on a sequence of univariate distributions. Additional results of time series in the Euclidean space are provided in Appendix A of the supplement.

In addition to the wild bootstrap method (denoted by Boot) introduced in the previous section, we also apply the standard permutation test (denoted by Permt) to simulate the critical value for the proposed CvM test. Furthermore, we outline the results of KS in (9) using the same bootstrapped or permutated samples of $S_{n}^{*}(\cdot)$. For each data structure, both bootstrap and permutation conduct $R=300$ replicates. The empirical rejection rates in Section 5.1 and Section 5.3 are averaged over $M=2500$ Monte Carlo replicates, whereas those in Section 5.2 are based on $M=1000$ replicates due to the computational expenses.

5.1 Functional time series

We consider the scenario where $\Omega=\{f:[0,1]\mapsto\mathbb{R},\int_{0}^{1}f^{2}(\tau)d\tau<\infty\}$ with $L_{2}$ metric, i.e. d(f,g)=\big{(}{\int_{0}^{1}(f(t)-g(t))^{2}dt}\big{)}^{1/2} for any $f,g\in\Omega$. Let $\mathcal{B}:=\{B(\tau):\tau\in[0,1]\}$ denote the standard Brownian motion, and we mimic the setups in Zhang, (2016) to generate $\{Y_{t}(\tau):\tau\in[0,1]\}_{t=1}^{n}$ from the following data generating processes (DGPs):

1. (i)

BM: $Y_{t}(\tau)=\varepsilon_{t}(\tau)$ for $\tau\in[0,1]$, where $\{\varepsilon_{t}(\tau)\}_{t=1}^{n}\stackrel{{\scriptstyle i.i.d.}}{{\sim}}\mathcal{B}$; 2. (ii)

BB: $Y_{t}(\tau)=\varepsilon_{t}(\tau)-\tau\varepsilon_{t}(1)$ for $\tau\in[0,1]$, where $\{\varepsilon_{t}(\tau)\}_{t=1}^{n}\stackrel{{\scriptstyle i.i.d.}}{{\sim}}\mathcal{B}$; 3. (iii)

FARCH: $Y_{t}(\tau_{1})=\epsilon_{t}(\tau_{1})\sqrt{\tau_{1}+\int_{0}^{1}\rho\exp\left({(\tau_{1}^{2}+\tau_{2}^{2})/2}\right)Y_{t-1}^{2}(\tau_{2})d\tau_{2}}$ for $\tau\in[0,1]$, where $\{\varepsilon_{t}(\tau)\}_{t=1}^{n}\stackrel{{\scriptstyle i.i.d.}}{{\sim}}\mathcal{B}$. We consider $\rho\in\{0.3418,0.6153\}$, which corresponds to the cases when the Hilbert–Schmidt norm of the Gaussian kernel equals 0.5 and 0.9, respectively. 4. (iv)

FNMA: $Y_{t}(\tau)=\varepsilon_{t}(\tau)\varepsilon_{t-1}(\tau)$ for $\tau\in[0,1]$, where $\{\varepsilon_{t}(\tau)\}_{t=1}^{n}$ is either BM or BB. 5. (v)

FAR: $Y_{t}(\tau_{1})=\int_{0}^{1}\mathcal{K}(\tau_{1},\tau_{2})Y_{t-1}(\tau_{2})d\tau_{2}+\varepsilon_{t}(\tau_{1})$, where $\{\varepsilon_{t}(\tau)\}_{t=1}^{n}$ is generated by either BM or BB, and we consider the Gaussian kernel $\mathcal{K}_{G}(\tau_{1},\tau_{2})=c_{G}\exp\left({(\tau_{1}^{2}+\tau_{2}^{2})/2}\right)$ and the Wiener kernel $\mathcal{K}_{W}(\tau_{1},\tau_{2})=c_{W}\min\{\tau_{1},\tau_{2}\}$. For each kernel, we choose the coefficient $c$ so that the corresponding Hilbert–Schmidt norm equal to 0.3, i.e. $c_{G}=0.2051$ for $\mathcal{K}_{G}$ and $c_{W}=0.7346$ for $\mathcal{K}_{W}$.

Among the existing tests that target at the serial dependence within a sequence of functional data, we consider those proposed in Gabrys and Kokoszka, (2007), Horváth et al., (2013) and Zhang, (2016) for comparison. In particular,

1. (i)

Gabrys and Kokoszka, (2007) (denoted by GK) proposed the portmanteau test based on the functional principal component analysis (fPCA). We consider the lag truncation number $h\in\{1,3,5\}$ and select the smallest $P$ such that the cumulative variation explained by the first $P$ principal components is above 90%. 2. (ii)

Horváth et al., (2013) (denoted by HHHR) presented an independence test based on the empirical correlation functions. We aggregate the recommendations of Horváth et al., (2013) and Zhang, (2016) to consider the lag truncation number $h\in\{5,10,30,50\}$. Each integral involved in HHHR is approximated by the Riemann sum with $p=100$ points. 3. (iii)

Zhang, (2016) proposed a Cramér–von Mises type test (denoted by Zhang) based on the functional periodogram. The critical value of Zhang is approximated using block bootstraps with $R=300$ replicates, and we consider the block size $b\in\{1,5,10\}$ as well as the optimal size $b^{\ast}$ chosen by the minimal volatility method. All the integrals are approximated by the Riemann sum with $p=50$ points to ease the computational burden.

We fix $n=200$ and generate the data on a grid of $T=1000$ equally spaced points in $[0,1]$ for each functional observation. We use Fourier and B-splines (order 4) with $B=20$ basis functions to obtain the functional data.

Table 1 includes the empirical rejection rate under each DGP, where P and B are shorthands for Permt and Boot respectively, and F and B denote the Fourier basis and B-splines respectively. Overall, the type of basis functions has limited influence on the finite sample performance for all the methods in Table 1. Under the null models BM and BB, the empirical size of CvM and KS are both close to 0.05, and the wild bootstrap has similar size accuracy with the permutation test. In contrast, the size accuracy of GK and HHHR depends on the lag truncation parameter $h$, and moderate size distortion can be observed with a larger number of lags. We also note that Zhang achieves decent size accuracy with a smaller block size $b$, but has some size distortion with a large $b$.

As for the power behaviors, our tests, especially with Permt, demonstrate dominating performance in FNMA and FAR models. While they are slightly inferior in FARCH(0.5), the performances get better when dependence level increases in FARCH(0.9). In comparison, both GK and HHHR seem to be sensitive to truncation parameter selections, leading to a considerable loss in power, especially with a larger lag truncation parameter. Note that both FARCH and FNMA models have serial dependence but no serial correlation, thus corresponding to null models in Zhang, (2016). Consequently, the power of Zhang against these models is close to the nominal level. We also note that GK and Zhang can have comparable power with our tests but only when a smaller bandwidth (or block size) parameter is used. Overall, our tests are among the best across all settings.

5.2 Covariance matrix

Next, we consider the scenario where the observations are a sequence of covariance matrices $\{Y_{t}\}_{t=1}^{n}\in\mathbb{R}^{p\times p}$, which is relatively new in literature.

We generate the time series from a conditional autoregressive Wishart model (Golosnoy et al.,, 2012) with dimension $p\in\{2,5,8\}$. In particular, let $\mathcal{W}_{p}(d,V)$ denote the Wishart distribution with degrees of freedom $d$ and the scale matrix $V$, and we generate $\{\varepsilon_{t}\}_{t=1}^{n}\stackrel{{\scriptstyle i.i.d.}}{{\sim}}\mathcal{W}_{p}(10,I_{p})\in\mathbb{R}^{p\times p}$. Furthermore, we generate $Y_{t}=\frac{1}{10}\mbox{chol}(\Sigma_{t})\varepsilon_{t}\mbox{chol}(\Sigma_{t})$, where $\mbox{chol}(\Sigma_{t})$ denotes the lower-triangular Cholesky factor of $\Sigma_{t}$, and $\Sigma_{t}=CC^{\top}+\rho AY_{t-1}A^{\top}+\rho B\Sigma_{t-1}B^{\top}$ with $A=0.7I_{p},$ $B=0.5I_{p}$, and $C=I_{p}$. Here we use the parameter $\rho\in\{0,0.5,75\}$ to control the temporal dependence and it is trivial that $\rho=0$ corresponds to the null model. We use $\text{CAW}_{p}(\rho)$ to denote the model with dimension $p$ and parameter $\rho$.

Note that in this case, we have $\Omega=\{A:\mbox{symmetric, positive definite matrix}\in\mathbb{R}^{p\times p}\}$. To evaluate the distance between two covariance matrices $A,B\in\Omega$, we consider the following metrics: (i) Euclidean metric: $d(A,B)=\|A-B\|_{F}$; (ii) log-Euclidean metric: $d(A,B)=\|\mbox{logm}(A)-\mbox{logm}(B)\|_{F}$, where $\mbox{logm}(\cdot)$ denotes the logarithm of a matrix; (iii) Cholesky metric: $d(A,B)=\|\mbox{chol}(A)-\mbox{chol}(B)\|_{F}$, where $\mbox{chol}(A)$ denotes the Cholesky decomposition of $A$; (iv) Riemann metric: $d(A,B)=\|\mbox{logm}(A^{-1/2}BA^{-1/2})\|_{F}$, where $A^{1/2}$ denotes the matrix square root of $A$.

In view of the fact that there is little literature regarding this testing problem, we use the half vectorization to convert a $p\times p$ covariance matrix into a $\frac{1}{2}p(p+1)$-dimensional vector and apply the multivariate testing methods to the resulting time series. In addition to our proposed test with the Euclidean metric (denoted by Vech-Euc), we also adopt the test statistic mADCV proposed in Fokianos and Pitsillou, (2018), which can be implemented via the R package dCovTS. In particular, mADCV is based on the auto-distance covariance matrix and requires a kernel function $k(\cdot)$ and a bandwidth parameter $p_{n}$. We follow the setup of Fokianos and Pitsillou, (2018) to consider the truncated kernel (TC), the Daniell kernel (DAN), the Parzen kernel (PAR), and the Bartlett kernel (BAR). Here, we set $p_{n}=\lfloor{3n^{\lambda}}\rfloor$ with $\lambda\in\{1/10,2/10,3/10\}$.

According to Table 2, our test achieves accurate size, regardless of the dimension $p$ and the choice of metric. Under the alternative, our proposed test with the Euclidean metric demonstrates the highest power when applied to either matrix-valued time series or half-vectorized multivariate time series. Interestingly, the power increases as the dimension $p$ grows, indicating its advantage in the moderate-dimensional scenario. On the other hand, the log-Euclidean metric and the Riemann metric exhibit high power when $p$ is small, but they do show a slight power loss against moderate values of $p$. This discrepancy in performance could be attributed to the fact that different metrics excel at capturing distinct dependence structures.

In contrast, mADCV has a conservative size under the null and have noticeable power loss under the alternative, especially when $p$ is moderate and $\rho$ is small. When $p=5$ and $p=8$, the half-vectorized covariance matrix is converted to the multivariate time series in $\mathbb{R}^{15}$ and $\mathbb{R}^{36}$, respectively. Note that mADCV involves the estimation of the long-run covariance matrix, and the estimation tends to be less accurate as the dimension grows, which may explain why mADCV worsens against growing $p$.

5.3 Univariate distribution

Lastly, we consider the case where $\Omega$ is the set of cumulative distribution functions of a random variable that takes value from $[0,1]$. In this case, we mimic the settings in Zhu and Müller, 2023a to generate a sequence of CDFs $\{Y_{t}(x):x\in[0,1]\}_{t=1}^{n}$, with $n=200$, from an autoregressive transport model of order $p$ denoted by ATM(p).

Specifically, let $g(x)$ denote the natural cubic spline passing through points $(0,0)$, $(0.33,0.2)$, $(0.66,0.8)$, $(1,1)$, and we generate a random sample $\{\xi_{t}\}_{t=1}^{n}$ from a uniform distribution over $(-1,1)$. For $x\in\mathcal{S}:=[0,1]$, we define $h(x)=\frac{1}{2}\left({(1-\xi_{t})g(x)+(1+\xi_{t})x}\right)$ and $\varepsilon_{t}(x)=\frac{1}{2}\left({(1+\xi_{t})g(h^{-1}(x))+(1{-}\xi_{t})h^{-1}(x)}\right)$. Then we generate a sequence of quantile functions $\{T_{t}\}_{t=1}^{n}$ by (i) ATM(0): $T_{t}=\varepsilon_{t}$; (ii) ATM(1): $T_{t}=\beta_{1}\odot T_{t-1}\oplus\varepsilon_{t}$; (iii) $T_{t}=\beta_{4}\odot T_{t-4}\oplus\beta_{3}\odot T_{t-3}\oplus\beta_{2}\odot T_{t-2}\oplus\beta_{1}\odot T_{t-1}\oplus\varepsilon_{t}$, where $(T_{1}\oplus T_{2})(x)=T_{2}\circ T_{1}(x)=T_{2}(T_{1}(x))$, and for $\beta\in[-1,1]$,

[TABLE]

For ATM(1), we consider $\beta_{1}=0.5\rho$ with $\rho\in\{0.2,0.5,1\}$, whereas for ATM(4), we consider $\beta=\rho(0.2,-0.5,0.1,-0.3)^{\top}$ with $\rho\in\{0.5,0.8,1\}$. Based on the simulated quantile functions $\{T_{t}(x):x\in\mathcal{S}\}_{t=1}^{n}$, we generate $\{Y_{t}(x):x\in[0,1]\}_{t=1}^{n}$ by $Y_{t}(x)=T_{t}^{-1}(x)$.

For the metric space $\Omega$, we consider several metrics to evaluate the distance between two distributions. In particular, we use $F,G$ to denote two CDFs, use $F^{-1},G^{-1}$ to denote the corresponding quantile functions, and use $f,g$ to denote the corresponding density functions, then the metrics are given by (i) 1-Wasserstein (W1): $d_{W1}(F,G)=\int_{0}^{1}\left|{F^{-1}(t)-G^{-1}(t)}\right|\mathrm{d}t$; (ii) 2-Wasserstein (W2): d_{W2}(F,G)=\big{(}{\int_{0}^{1}\big{(}{F^{-1}(t)-G^{-1}(t)}\big{)}^{2}\mathrm{d}t}\big{)}^{1/2}; (iii) Kolmogorov–Smirnov (KS): $d_{KS}(F,G)=\sup_{x\in\mathcal{S}}\left|{F(x)-G(x)}\right|$; (iv) Kullback–Leibler (KL): d_{KL}(F,G)=\frac{1}{2}\big{(}{\int_{x\in\mathcal{S}}\log\frac{F(\mathrm{d}x)}{G(\mathrm{d}x)}F(\mathrm{d}x)+\int_{x\in\mathcal{S}}\log\frac{G(\mathrm{d}x)}{F(\mathrm{d}x)}G(\mathrm{d}x)}\big{)}; (v) Itakura–Saito (IS): $d_{IS}(F,G)=\frac{1}{2|\mathcal{S}|}\int_{x\in\mathcal{S}}\left({\frac{f(x)}{g(x)}-\log\frac{f(x)}{g(x)}-1}\right)\mathrm{d}x+\frac{1}{2|\mathcal{S}|}\int_{x\in\mathcal{S}}\left({\frac{g(x)}{f(x)}-\log\frac{g(x)}{f(x)}-1}\right)\mathrm{d}x$; (vi) Log-Spectral (LS): $d_{LS}(F,G)=\left({\frac{1}{|\mathcal{S}|}\int_{x\in\mathcal{S}}\left({10\log_{10}\frac{f(x)}{g(x)}}\right)^{2}\mathrm{d}x}\right)^{1/2}$.

In practice, we compute the Wasserstein distances directly using the simulated quantile functions $\{T_{t}(x)\}_{t=1}^{n}$, and compute all the other distances using R package seewave with the probability functions estimated from the simulated CDFs. To the best of our knowledge, there is no existing literature for testing the serial dependence within a sequence of univariate distributions, hence we only present the empirical rejection rates of the proposed method in Table 3.

The overall pattern shown in Table 3 generally matches those observed in Section 5.1 and Section 5.2. Specifically, the type of test statistic and the method used to approximate the critical value have little influence on both the size and the power. Among all the metrics examined in Table 3, we observe no significant disparities in terms of the size accuracy, whereas the Itakura-Saito metric and the Kullback-Leibler metric exhibit the dominating power in most cases.

6 Applications

In this section, we demonstrate the versatility of the proposed tests in two real applications.

6.1 Cumulative intraday returns

Modelling financial returns using functional data approach has attracted much attention in the literature, see e.g. Aue et al., (2017), Shang, (2017), Cerovecki et al., (2019). In the first application, we apply the proposed tests to the cumulative intraday returns (CIDRs) of two stocks of International Business Machines Corporation (IBM) and McDonald’s Corporation (MCD) in 2007. Following Gabrys et al., (2010), we define the CIDRs as

[TABLE]

where $P_{t}(s)$ denote the stock price at rescaled trading period $s$, and $P_{t}(0)$ denote the market opening price. The daily trading data is collected at 1-minute frequency with a total 390 data points per day. The data is then smoothed using 20 Fourier basis functions. We consider the proposed CvM test and KS test, along with other competing methods, including GK at lag $h={1,3,5}$, HHHR with truncation lag $h={5,10,30,50}$, and Zhang with $b={1,5,10}$. The p-values of all tests are summarized in Table 4.

From Table 4, we find that all tests except for HHHR50 reject $H_{0}$ at 5% level for IBM stock, and all except for GK3 fail to reject $H_{0}$ for MCD. This finding is consistent with Zhang, (2016) using 5-minute data and cubic B-splines basis functions. Here, we note that the constructed price process could be contaminated by microstructure noise (Aït-Sahalia and Yu,, 2009) and we conjecture that their presence could compromise power due to diminished signal-to-noise ratio. We leave this for future research.

6.2 Mortality Data

Human Mortality Database (https://www.mortality.org/Home/Index) has provided the scientific researchers a high-quality harmonized mortality and population estimates. For example, Figure 1 below gives a visual demonstration of the female (of age between 20 and 105) mortality distribution time series in U.K. from 1970 to 2014.

The yearly age-at-death distribution for a given country can be naturally viewed as a random element in the space of univariate distributions, which has been analyzed in many studies regarding modelling aspects for non-Euclidean valued random objects, see, e.g. Petersen and Müller, (2019), Dubey and Müller, (2020).

However, many of them implicitly assume temporal independence for the sequence of random distributions, which could be doubtful. We now formally test for this assumption.

We analyze the yearly female mortality data in 27 European developed countries from 1970 to 2014. After applying our proposed CvM test and KS test with Wasserstein-1 and Wasserstein-2 metrics, we find that all tests for all the countries reject $H_{0}$ at 5% level. This provides strong evidence of temporal dependence in mortality distribution time series. Therefore, more caution should be taken when it comes to the modelling of mortality data, as temporal dependence may not be negligible.

7 Conclusion

In this paper, we have proposed a new method of testing serial independence for object-valued time series. It builds upon the distance covariance in metric spaces, and captures all pairwise dependence. It is fully nonparametric, easy to implement with no tuning parameters, and is broadly applicable to many different data types as long as they live in a metric space of strong negative type. Our numerical results demonstrate accurate size for both permutation and bootstrap based tests and very competitive power performance even in conventional Euclidean and Hilbert spaces.

This paper is concerned with serial dependence testing of object-valued time series, which is typically conducted before any specific model is fitted. There is a recent literature on modelling temporal dependence in object-valued time series, e.g. distributional autoregression in Zhang et al., (2022), Zhu and Müller, 2023a , Ghodrati and Panaretos, (2023), and spherical autoregression in Zhu and Müller, 2023b , among others. It would be interesting to investigate how to extend the test in this paper to model diagnostic checking. This requires the researchers to properly define model residuals, which is a challenging task for object-valued data in metric spaces. In addition, we assume the object-valued data are fully observed in the paper whereas in real applications they are usually generated by pre-processing Euclidean-valued data. In the realm of Hilbert space, Zhang and Wang, (2016) have demonstrated that pre-smoothing in functional data analysis can exhibit a phase transition phenomenon, depending on the resolution of discretized observations. It would also be helpful to investigate how pre-processing would impact the behavior of our test. Finally, we assume the distance metric $d(\cdot,\cdot)$ is given a priori but its choice can play a significant role in data analysis; it usually depends on the underlying inference problem and the characteristics of the dataset. How to choose the distance metric adaptively is very important yet also very challenging. We leave these topics for future research.

\appendixpage

The supplementary material is organized as follows. Appendix A contains additional simulation studies of the Euclidean time seires. Appendix B provides some auxiliary lemmas and their proofs, and Appendix C gives all the proofs of theorems in the main text.

Appendix A Additional Numerical Studies

In this section, we carry out some additional simulation studies based on the Euclidean time series. Appendix A.1 includes the simulation results of the univariate time series whereas Appendix A.2 contains those of the bivariate time series. In Appendix A.3, we further investigate the finite sample performance of the proposed test statistics against different sample sizes $n$ and dimension $p$. Throughout this section, all the empirical rejection rates are average over $M=2500$ Monte Carlo replicates following the same computing procedures used in Section 5.

A.1 Univariate time series

Let $\Omega=\mathbb{R}$ and $d(x,y)=|x-y|$ for any $x,y\in\mathbb{R}$. We follow the setups of (Fokianos and Pitsillou,, 2017) to generate a time series $\{Y_{t}\}_{t=1}^{n}\in\mathbb{R}$. Specifically, we set $n=200$ and consider four data generating processes (DGPs): (i) IID: $Y_{t}=\epsilon_{t}$; (ii) NMA(2): $Y_{t}=\epsilon_{t}\epsilon_{t-1}\epsilon_{t-2}$; (iii) ARCH(2): $Y_{t}=\sigma_{t}\epsilon_{t}$, $\sigma_{t}^{2}=0.5+0.8Y_{t-1}^{2}+0.1Y_{t-2}^{2}$; (iv) TAR(1): Y_{t}=\left\{\begin{array}[]{ll}-1.5Y_{t-1}+\epsilon_{t},&Y_{t-1}<0\\ 0.5Y_{t-1}+\epsilon_{t},&Y_{t-1}\geq 0\end{array}\right., where $\{\epsilon_{t}\}\stackrel{{\scriptstyle i.i.d.}}{{\sim}}\mathcal{N}(0,1)$ in each DGP.

Apart from that of the proposed test statistic, we also report the performance of the following methods with a wide range of parameters for comparison:

1. (i)

the classic Box–Pierce test (denoted by BP) and Ljung–Box test (denoted by LB). Both BP and LB can be implemented using R function Box.test and we consider the number of lags $K\in\{1,3,6,9,12,15\}$; 2. (ii)

the consistent testing for serial correlation proposed by Hong, (1996) (denoted by Hong) with a pre-specified kernel $k(\cdot)$ and a bandwidth parameter $p_{n}$. In particular, we use the Daniell kernel (D), the Parzen kernel (P), the Bartlett kernel (B), and the QS kernel (Q). As for the selection of $p_{n}$, it is suggested by Hong, (1996) to select $p_{n}\in\{\lfloor{\log(n)}\rfloor,\lfloor{3*n^{0.2}}\rfloor,\lfloor{3*n^{0.3}}\rfloor\}$, which corresponds to $p_{n}\in\{5,9,15\}$ with $n=200$. Here, we consider a larger set of $p_{n}$ with $p_{n}\in\{3,6,9,12,15\}$. 3. (iii)

the automatic Portmanteau test in Escanciano and Lobato, (2009) (denoted by EL) with an upper bound $d\in\{25,50\}$. In Escanciano and Lobato, (2009), the recommended value of $q$ is 2.4. We consider a larger selection of $q\in\{1.8,2.1,2.4,2.7,3.0\}$. 4. (iv)

the consistent testing for pairwise dependence by Fokianos and Pitsillou, (2017) (denoted by FP) with the same kernel function $k(\cdot)$ and bandwidth parameter $p_{n}$ as those used for Hong. In the original paper, the suggested selection of $p_{n}$ include $\{\lfloor{n^{0.2}}\rfloor,\lfloor{n^{0.4}}\rfloor,\lfloor{n^{0.6}}\rfloor\}$, which corresponds to $\{3,9,25\}$. In our setting, we consider $p_{n}\in\{3,9,15,21,25\}$.

Table 5 summarizes the empirical rejection rates of all the methods under each DGP. Here, H-D3 stands for the test statistic Hong with the Daniell kernel and $p_{n}=3$, EL25-2.4 stands for the test statistic EL with $d=25$ and $q=2.4$. Similarly, the other concatenations of initials represent the corresponding parameter combinations.

Under the null, our method together with LB both achieve accurate size. For BP and FP, the size appears conservative when the truncation number/bandwidth becomes large. We also observe some degree of over-rejection with H and EL for the range of tuning parameters under examination. The power of our test against the alternative depends on the specific DGP. In particular, the proposed test is among the most powerful tests against the TAR(2) model. Note that both NMA(2) and ARCH(2) models have serial dependence but no serial correlation, and our method has very high power against NMA(2) model but exhibits slight power loss against the ARCH(2) model when comparing to the most powerful test in this case (i.e., FP with $p_{n}=3$). This observation is consistent with that in Section 5.1.

As shown from the simulation results, the performance of all the competing methods depends on the choice of tuning parameters. It is still possible that our test is outperformed by one of these methods with carefully selected tuning parameters under specific DGP. For example, FP achieves a very accurate size and the highest power but only when specific parameters are chosen. Although empirical recommendations of tuning parameter selection are provided for each comparison method, there is no theoretical or practical guarantee, making parameter selection difficult in practice. In comparison, the proposed test is tuning-free and is convenient to implement, which is a natural advantage in real-world applications.

A.2 Bivariate time series

Consider $\Omega=\mathbb{R}^{2}$ with the Euclidean distance $d(x,y)=|x-y|_{2}$, and we generate $\{\epsilon_{t}\}\stackrel{{\scriptstyle i.i.d.}}{{\sim}}\mathcal{N}_{2}(0,\Sigma)$ with \Sigma=\left(\begin{array}[]{cc}1&\rho\\ \rho&1\\ \end{array}\right) and $\rho\in\{0,0.4,0.7\}$. We fix $n=200$ and generate $\{Y_{t}\}_{t=1}^{n}$ by: (i) IID: $Y_{t}=\epsilon_{t}$; (ii) NMA(2): $Y_{t,i}=\epsilon_{t,i}\epsilon_{t-1,i}\epsilon_{t-2,i}$; (iii) ARCH(2): $Y_{t,i}=h_{t,i}^{1/2}\epsilon_{t,i}$, where

[TABLE]

and (iv) MAR(2): Y_{t}=\left(\begin{array}[]{cc}0.04&-0.1\\ 0.11&0.5\\ \end{array}\right)Y_{t-1}+\epsilon_{t}.

The multivariate Ljung–Box test statistic mLB proposed in Hosking, (1980) is used for comparison and we select the number of lags $h$ from $\{1,3,6,9\}$. Additionally, we apply the multivariate testings mADCV and mADCF in Fokianos and Pitsillou, (2018) using the same kernel functions and bandwidth parameters as mentioned in Section 5.2.

All the simulation results are reported in Table 6. Similar to the observations in Appendix A.1, KS and CvM have comparable size accuracy and power across the table, and for both types of statistics, Permt yields higher power than Boot on the serially-uncorrelated time series NMA(2) and ARCH(2). Additionally, both test statistics seem robust to the componentwise dependence within the data, implying that the proposed test can handle various dependence structures.

By contrast, mLB exhibits noticable size distortion and substantial power loss. As for mADCV and mADCF, both tests are accurate in size but have slight power loss when compared to our proposed test statistic, especially when using the truncated kernel. The choice of the bandwidth also has some impact on the finite sample performance of mADCV and mADCF. In particular, both test statistics have higher power with a smaller $\lambda$ under NMA(2) whereas with a larger $\lambda$ under ARCH(2). In addition, the strengthening componentwise dependence of the NMA(2) model leads to a more severe power loss of mADCV and mADCF.

A.3 Multivariate time series

Lastly, we mimic Table 1 in Zhou, (2012) to investigate the impact of $n$ and $p$ on the power of our test. To this end, We generate a $p$-dimensinoal time series from a VAR(1) model, i.e.

[TABLE]

where $\{\varepsilon_{t}\}_{t=1}^{n}\stackrel{{\scriptstyle iid}}{{\sim}}\mathcal{N}_{p}(0,I_{p})$ is generated from the standard normal distribution in $\mathbb{R}^{p}$. Here, we consider $n\in\{100,250\}$, $p\in\{2,5,10,20,30\}$ and we fix $\rho=0.2$. The comparison methods are the same as in Appendix A.2.

According to the numerical results in Table 7, the power of all the methods increases significantly as $n$ increases. As the dimension $p$ increases, the rejection rates of the proposed methods and mLB both increase. In particular, when $p$ is small, the proposed test could be outperformed by mLB using carefully selected parameters. However, under the moderate-dimensional setting, our method has the highest power, regardless of the sample size.

As $p$ increases, we observe noticeable power loss of mADCV and mADCF, which is consistent with our observation in Section 5.2. We conjecture this is due to deteriorated long-run covariance estimation when dimension increases. This is a new finding as the simulation studies in Fokianos and Pitsillou, (2018) are limited to the two-dimensional case.

Appendix B Auxiliary Lemmas

Lemma B.1**.**

For the kernel $h$ in (6), let $z_{i}=(x_{i},y_{i})$ and $Z_{i}=(X_{i},Y_{i})$ i.i.d., $i=1,\cdots,4$ and

[TABLE]

where $c=1,2,3,4$, and define $\mathbb{E}h(Z_{1},\cdots,Z_{4})=\theta$. Then,

[TABLE]

Furthermore, let $h^{(1)}(z)=h_{1}(z)-\theta$ and $h^{(c)}=h_{c}(z_{1},\cdots,z_{c})-\sum_{j=1}^{c-1}\sum_{(c,j)}h^{(j)}(z_{i_{1}},\cdots,z_{i_{j}})-\theta$, where ${\sum}_{(c,j)}$ denotes the summation over all $j$-subsets $\{({i_{1}},\cdots,{i_{j}})\}$ of $\{1,\cdots,c\}$. Then,

[TABLE]

where $U_{n}(h^{(c)})$ is the $U$ statistic based on $h^{(c)}$.

Proof.

See Section 9 in Han and Shen, (2021) or Section 1.2 of the supplement of Zhang et al., (2018). (11) follows directly by Hoeffding decomposition. ∎

Lemma B.2**.**

Denote by $\{\lambda_{\ell}\}_{\ell=1}^{\infty}$, and $\{e_{\ell}\}_{\ell=1}^{\infty}$ as the eigenvalues and orthonormal eigenfunctions corresponding to (7), then under Assumptions 3.1,

[TABLE]

where the series converges absolutely and uniformly on $(z,z^{\prime})\in\Omega^{2}\times\Omega^{2}$.

Proof.

When $(\Omega,d)$ is of strong negative type, by Lemma 12 and Proposition 18 in Sejdinovic et al., (2013), we know that $-d_{\nu}(\cdot,\cdot)$ is a symmetric positive definite kernel induced by $d$. Therefore, we obtain the product kernel $\mathcal{K}(\cdot,\cdot)$ (Steinwart and Christmann,, 2008, Lemma  4.6 ). By Cauchy-Schwarz inequality and triangle inequality,

[TABLE]

Then, by Proposition 1 and Theorem 2 in Sun, (2005), the result follows. ∎

Lemma B.3**.**

Define $\mathcal{F}_{t}=\sigma(X_{t},X_{t-1},\cdots)$, and $Z_{t}^{(k)}=(X_{t},X_{t-k})$. Under $H_{0}$, for any fixed $K$ and $M$, we have (i). $\left(\{e_{\ell}(Z_{t}^{(k)})\}_{\ell=1,\cdots,M;k=1,\cdots,K},\mathcal{F}_{t-1}\right)$ forms a sequence of martingale difference vectors; (ii). $\mathbb{E}[e_{\ell_{1}}(Z_{t}^{(k_{1})})e_{\ell_{2}}(Z_{t}^{(k_{2})})]=\mathbf{1}(\ell_{1}=\ell_{2},k_{1}=k_{2})$; (iii). for $i_{1}<j_{1},i_{2}<j_{2}$, $\mathbb{E}\left\{e_{\ell_{1}}(Z_{i_{1}}^{(k_{1})})e_{\ell_{1}}(Z_{j_{1}}^{(k_{1})})e_{\ell_{2}}(Z_{i_{2}}^{(k_{2})})e_{\ell_{2}}(Z_{j_{2}}^{(k_{2})})\right\}=0$ except for the case when the index subset $\{i_{1}-k_{1},i_{1},j_{1}-k_{1},j_{1}\}$ is identical to $\{i_{2}-k_{2},i_{2},j_{2}-k_{2},j_{2}\}$.

Proof.

(i). The M.D.S. claim follows by verifying for each $k$ and $\ell$. We first show that

[TABLE]

where $\mathcal{K}(z,z^{\prime})$ is defined in (12) of B.2.

By direct calculation,

[TABLE]

where the second last equality holds by the independence between $X_{t}$ and $\mathcal{F}_{t-1}$, and the last by the fact that $\mathbb{E}[d_{\nu}(x,X_{t})]=0$.

By (7), we have

[TABLE]

taking conditional expectation w.r.t. $\mathcal{F}_{t-1}$ on both sides of the above equation we can obtain that

[TABLE]

Note $\lambda_{\ell}$ is a constant, the result follows.

(ii). Under $H_{0}$, for any $k_{1}\neq k_{2}$, we have that

[TABLE]

where the second inequality holds by the fact that $\{X_{t}\}$ is i.i.d., and the third by law of iterated expectation and the last by (13). For $k_{1}=k_{2}$ and $\ell_{1}\neq\ell_{2}$, we know that $\mathbb{E}[e_{\ell_{1}}(Z)e_{\ell_{2}}(Z)]=0$ by the orthogonality of eigenfunctions.

(iii). Note that

[TABLE]

By similar arguments as (ii), we know that if one element in $\{i_{1}-k_{1},i_{1},j_{1}-k_{1},j_{1},i_{2}-k_{2},i_{2},j_{2}-k_{2},j_{2}\}$ is not paired, using the M.D.S. property we can claim the expectation is zero. In addition, when $i_{1}<j_{1}$ and $i_{2}<j_{2}$, we know that ${i_{1}-k_{1}}<{i_{1}},j_{1}-k_{1}<{j_{1}}$, and ${i_{2}-k_{2}}<{i_{2}},j_{2}-k_{2}<{j_{2}}$. This suggests we must have $i_{1}-k_{1}=i_{2}-k_{2}$, and $j_{1}=j_{2}$. When $k_{1}=k_{2}$, for the nonnull expectation, we must have $i_{1}=i_{2}$; and when $k_{1}\neq k_{2}$, we must have $i_{2}=j_{1}-k_{1}$ and $i_{1}=j_{2}-k_{2}$. In either case, the index subset $\{i_{1}-k_{1},i_{1},j_{1}-k_{1},j_{1}\}$ is identical to $\{i_{2}-k_{2},i_{2},j_{2}-k_{2},j_{2}\}$.

∎

Lemma B.4**.**

Under $H_{0}$, $\mathbb{E}V_{n}(k)=O((n-k)^{-2})$.

Proof.

In view of B.2 and B.3 (ii), we first claim that except for the case that $(i_{1},i_{2})=(j_{1},j_{2})$ or $(i_{1},i_{2})=(j_{2},j_{1})$,

[TABLE]

By Hoeffding decomposition in (11), under $H_{0}$, we have $V(k)=0$, $h_{1}=0$ and $\mathbb{E}h_{2}=0$, $h^{(1)}=0$ and $\mathbb{E}h^{(2)}=0$. Therefore it suffices to consider $\mathbb{E}[h_{c}(Z_{i_{1}}^{(k)},\cdots,Z_{i_{c}}^{(k)})]$ for $c=3,4,$ and $k+1\leq i_{1}<\cdots,i_{c}\leq n$.

For $c=3$, for $s<u<v$, we have

[TABLE]

Clearly, by (14), $\mathbb{E}h_{31}=\mathbb{E}h_{32}=\mathbb{E}h_{34}=0$. Recall for $z_{\ell}=Z_{\ell}^{(k)},$ $\ell\in(s,u,v)$, we have

[TABLE]

where the nonnull expecation only appears in terms satisfying $(i_{3}-k,i_{1}-k)=(i_{1},i_{2})$, suggesting $s=u-k=v-2k$. Note there are at most $O(n-k)$ such $(s,u,v)$-tuples that can satisfy the constraint. This implies that in the summands of U-statistics $U_{n,k}(h^{(3)})$, there are only $O(n-k)$ terms that have nonnull expectation, hence $\mathbb{E}U_{n,k}(h^{(3)})=O((n-k)^{-2})$. Similar arguments also applies to $h_{4}$.

∎

Lemma B.5**.**

Under $H_{0}$, suppose 3.1 holds. Then, for $U_{n,k}(h^{(c)}),c=3,4,$ defined in (11), with any fixed $k$,

[TABLE]

for some uniform constant $C>0$ independent of $n$ and $k$.

Proof.

Note that $U_{n,k}(h^{(c)})$ is a U-statistic of order $c$ based on samples $Z_{i}^{(k)}=(X_{i},X_{i-k})$, $k+1\leq i\leq n$. This implies that $\{Z_{i}^{(k)}\}$ is $k$-dependent. Following the treatment in Janson, (2021), we define

[TABLE]

such that in above equation, $Z_{i_{1}}^{(k)},\cdots,Z_{i_{c}}^{(k)}$ are independent of each other. This implies that $U_{n,k}(h^{(c)};>k)$ are formed based on independent copies of $Z_{i}^{(k)}$. By standard arguments in Hoeffding decomposion, see e.g., Theorem 3 in Chapter 1.6 of Lee, (1990), we know that $\mathrm{Var}\left(U_{n,k}(h^{(c)};>k)\right)\leq C(n-k)^{-c}$. Furthermore, by Lemma 4.4 in Janson, (2021), when $\mathbb{E}h^{2}(\cdot)<\infty$ (which is ensure by Assumption 3.1), we have that $\mathrm{Var}\left(U_{n,k}(h^{(c)})-U_{n,k}(h^{(c)};>k)\right)\leq C(n-k)^{-3}$. Hence, by Cauchy-Schwarz inequality, we know that

[TABLE]

∎

Appendix C Technical Proofs

Proof of Theorem 3.1

(i). Leading term.

Under $H_{0}$, we have $V_{X}(k)=0$, $h_{1}(z)=0$ and $h^{(2)}(z,z^{\prime})=h_{2}(z,z^{\prime})=d_{\nu}(x,x^{\prime})d_{\nu}(y,y^{\prime})/6=\mathcal{K}(z,z^{\prime})/6$, with the kernel $\mathcal{K}(z,z^{\prime})$ defined in (12). Hence, by Hoeffding decomposition in Lemma B.1,

[TABLE]

By B.4 and B.5, we know that $\mathbb{E}[U_{n,k}(h^{(c)})]^{2}=O((n-k)^{-3})$ for $c=3,4$, hence by Markov inequality,

[TABLE]

Since $K$ is fixed, in the following proof, with minor abuse of notation, we can approximately consider

[TABLE]

(ii). Approximation for $\mathcal{K}(z,z^{\prime})$.

By Lemma B.2, let

[TABLE]

Let

[TABLE]

Then, we have that,

[TABLE]

where the second equality holds by Lemma B.3 (iii), and the inequality holds by Cauchy-Schwarz inequality $\mathbb{E}e^{2}_{\ell_{1}}(Z_{i}^{(k)})e^{2}_{\ell_{2}}(Z_{j}^{(k)})\leq[\mathbb{E}e^{4}_{\ell_{1}}(Z)]^{1/2}[\mathbb{E}e^{4}_{\ell_{2}}(Z)]^{1/2}\leq\sup_{\ell}\mathbb{E}[e_{\ell}^{4}(Z)]<\infty$, which is ensured by $\mathbb{E}(\mathcal{K}^{2}(Z,Z))<\infty$. The final convergence is ensured by $\sum_{\ell=1}^{\infty}\lambda_{\ell}=\mathbb{E}\mathcal{K}(Z,Z)<\infty$.

(iii). Joint convergence of $n\left\{\tilde{V}_{n}^{(M)}(k)\right\}_{k=1}^{K}$.

Note by (17), we have

[TABLE]

In view of B.3 (i), we know that $\left(\{e_{\ell}(Z_{t}^{(k)}\}_{\ell=1,\cdots,M;k=1,\cdots,K}),\mathcal{F}_{t-1}\right)$ forms a sequence of martingale difference for any fixed $M$ and $K$. Hence $\frac{1}{{n-K}}\sum_{i=K+1}^{n}e_{\ell}^{2}(Z_{i}^{(k)})\to_{p}\mathbb{E}e_{\ell}^{2}(Z)=1$ by the weak law of large numbers, and

[TABLE]

by CLT for martingale difference sequences and Cramér Wold device, where $\{G^{(k)}_{\ell}\}_{\ell=1,\cdots,M;k=1\cdots,K}$ is a sequence of i.i.d. standard normal distribution in view of B.3 (ii).

Then continuous mapping theorem implies that

[TABLE]

The final result then follows by letting $M\to\infty$ in view of (18).

Proof of Theorem 3.2

Denote

[TABLE]

and for $\{\xi_{k}\}$ defined in Theorem 3.1,

[TABLE]

By Proposition 6.3.9 of Brockwell et al., (1991), to show the weak convergence of $S_{n}(\zeta)$, it suffices to show that (i) for each $K$, $S_{n}^{K}(\zeta)\Rightarrow S^{K}(\zeta)$ as $n\to\infty$; (ii) $S^{K}(\zeta)\Rightarrow S(\zeta)$ as $K\to\infty$; (iii) for any $\epsilon>0$, $\lim_{K\to\infty}\lim_{n\to\infty}\mathbb{P}(\|R_{n}^{K}\|>\epsilon)=0$.

(i). Note that (i) follows by showing (a) the finite dimensional convergence of $\{S_{n}^{K}(\zeta_{i})\}_{i=1}^{M}$ for arbitrary finite $M$, $(\zeta_{1},\zeta_{2},\cdots,\zeta_{M})$ with $\zeta_{i}\in[0,\pi]$; (b) the tightness of $S_{n}^{K}(\zeta)$.

First, as $K$ and $M$ are finite, the proof of (a) follows directly from Theorem 3.1 and continuous mapping theorem.

Second, for (b), we note that $S_{n}^{K}(\zeta)=\sum_{k=1}^{K}(n-k){V}_{n}(k)\Psi_{k}(\zeta)$, and when $K$ is fixed, we only need to show the tightness of $(n-k){V}_{n}(k)\Psi_{k}(\zeta)$ for each $k$. Fix a complete orthonormal basis in $\mathbb{H}$, denoted by $e_{1},e_{2},\cdots$. By Lemma 7.1 in Panaretos and Tavakoli, (2013), it suffices to show that

[TABLE]

In fact, for any $h\in\mathbb{H}$, denote $W_{h}(k)=\int_{[0,\pi]}h(\zeta)\Psi_{k}(\zeta)\mathrm{d}\zeta$. Then

[TABLE]

where we note that under Assumption 3.1, for any $n\geq k+4$, we have $\mathbb{E}[(n-k){V}_{n}(k)]^{2}\leq C<\infty.$ Hence,

[TABLE]

Hence, (b) holds, and this completes the proof for (i).

(ii) Holds trivially.

(iii) Recall

[TABLE]

By Chebyshev inequality, we have

[TABLE]

where the equality holds by noting that $\langle\Psi_{k_{1}},\Psi_{k_{2}}\rangle=0$ if $k_{1}\neq k_{2}$ and the last inequality holds using $\|\Psi_{k}\|^{2}\leq k^{-2}$.

Hence, (iii) follows by letting $K\to\infty$. This completes the proof.

∎

Proof of 3.3

It suffices to show that $V_{n}(k)\to_{p}V(k)$. Recall the alternative representation in (6), in view of Theorem 1 (iii) in Arcones, (1998), it suffices to show that for any $i,j,q,r$, for some $0<\delta\leq 1$,

[TABLE]

For any $x\in\mathbb{R}$, by the elementary inequality that $\log(1+x)\leq x$ for $x\geq 0$, we have

[TABLE]

Therefore, it suffices to show that for some $\delta>0$,

[TABLE]

In fact,

[TABLE]

where the first inequality holds by Minkowski inequality, the second by Hölder inequality, and the third by triangle inequality.

∎

Proof of 4.1

(i). Leading term by $\tilde{V}_{n}^{*}(k)={n-k\choose 2}^{-1}\sum_{i=k+1}^{n}\sum_{j=i+1}^{n}\mathcal{K}(Z_{i}^{(k)},Z_{j}^{(k)})w_{i}(k)w_{j}(k)$.

Recall $\mathcal{K}(Z_{i}^{(k)},Z_{j}^{(k)})=d_{\nu}(X_{i},X_{j})d_{\nu}(X_{i-k},X_{j-k}),$ with $\mathcal{K}(z,z^{\prime})$ defined in (12). It is clear that

[TABLE]

so we first want to show that

[TABLE]

where the equality holds using the joint independence of $\{w_{i}(k)\}$.

By Cauchy-Schwarz and Hölder inequality, we have that

[TABLE]

Under Assumption 3.1, it is clear that $\frac{1}{(n-k-3)^{2}}\sum_{k+1\leq i\neq j\leq n}d^{4}_{\nu}(X_{i},X_{j})$ and $\frac{1}{(n-k-3)^{2}}\sum_{k+1\leq i\neq j\leq n}d^{4}_{\nu}(X_{i-k},X_{j-k})$ are bounded in probability. This implies we only need to show

[TABLE]

as

[TABLE]

is similar. By Minkowski inequality, we have

[TABLE]

For the first term, under $H_{0}$, we note that

[TABLE]

for distinct 4-tuples $(i_{1},i_{2},i_{3},i_{4})$. This implies that the first term is at most of order $O(n^{-1})$ in view of Assumption 3.1. Similarly, the second term is at most of order $O(n^{-1})$. Hence $\mathbb{E}[\tilde{a}_{ij}-d_{\nu}(X_{i},X_{j})]^{4}\to 0$, and by Markov inequality, (22) holds. Therefore (21) holds.

(ii). Approximation.

Denote

[TABLE]

We want to show that as $M\to\infty$,

[TABLE]

Recall $\mathcal{K}^{(M)}(z,z^{\prime})=\sum_{\ell=1}^{M}\lambda_{\ell}e_{\ell}(z)e_{\ell}(z^{\prime})$, we have

[TABLE]

where we note that $Z_{i}^{(k)}$ and $Z_{j}^{(k)}$ are independent for $|i-j|>k$, and that there are at most $kn$ terms in the summation of $\sum_{\begin{subarray}{c}k+1\leq i<j\leq n\\ j\leq i+k\end{subarray}}$.

Hence, (24) follows by that $\sum_{\ell=M+1}^{\infty}\lambda_{\ell}^{2}\to 0$ as $M\to\infty$.

(iii). Joint Convergence of $\{(n-k)\tilde{V}_{n}^{X,(M)*}(k)\}_{k=1}^{K}\to_{d^{*}}\{\sum_{\ell=1}^{M}\lambda_{\ell}[G_{\ell}^{(k)}]^{2}-1\}_{k=1}^{K}$ in probability.

By continuous mapping theorem, it suffices to show that, in probability,

[TABLE]

Note that $\mathbb{E}^{*}[\frac{1}{n-k}\sum_{i=k+1}^{n}w_{i}^{2}(k)e_{\ell}^{2}(Z_{i}^{(k)})]=\frac{1}{n-k}\sum_{i=k+1}^{n}e_{\ell}^{2}(Z_{i}^{(k)})\to_{p}1$, and $\mathrm{Var}^{*}[\frac{1}{n-k}\sum_{i=k+1}^{n}w_{i}^{2}(k)e_{\ell}^{2}(Z_{i}^{(k)})]=\frac{1}{(n-k)^{2}}\sum_{i=k+1}^{n}e_{\ell}^{4}(Z_{i}^{(k)})\mathrm{Var}(w_{i}^{2}(k))\to_{p}0$ by weak law of large numbers for $k-$dependent sequences and Slutsky’s theorem. Hence, by Chebyshev inequality,

[TABLE]

Next, note that for any fixed $\ell$ and $k$, $\sum_{i=k+1}^{n}[(n-k)^{-1/2}e_{\ell}(Z_{i}^{(k)})]^{2}\to_{p}1$ and

[TABLE]

which implies that Lyapunov central limit theorem holds in probability.

Furthermore, for any fixed $\ell_{1},\ell_{2}$ and $k_{1}\leq k_{2}$, using the independence between $\mathbf{w}(k_{1})$ and $\mathbf{w}(k_{2})$ for $k_{1}\neq k_{2}$, we have

[TABLE]

in view of B.3(ii).

Therefore, by Cramér-Wold device, we obtain the joint convergence.

Finally, in view of (21) and (24), the result follows.

∎

Proof of 4.2

The proof is similar to 3.2. By continuous mapping theorem, it suffices to show that

[TABLE]

where $\Rightarrow^{*}$ represents the weak convergence in $L_{2}[0,\pi]$ under bootstrap asymptotics.

Denote

[TABLE]

To show the weak convergence of $S_{n}^{*}(\zeta)$ in probability, it suffices to show in probability (i). for each $K$, $S_{n}^{K*}(\zeta)\Rightarrow S^{K}(\zeta)$ as $n\to\infty$; (ii). $S^{K}(\zeta)\Rightarrow S(\zeta)$ as $K\to\infty$; (iii). for any $\epsilon>0$, $\lim_{K\to\infty}\lim_{n\to\infty}\mathbb{P}^{*}(\|R_{n}^{K*}\|>\epsilon)=0$.

(i) We need to show that (a). $\{S_{n}^{K*}(\zeta_{i})\}_{i=1}^{M}\to_{d}\{S^{K}(\zeta_{i})\}_{i=1}^{M}$ in probability; (b). $S_{n}^{K*}(\zeta)$ is asymptotically tight conditional on the sample.

The proof of (a) follows directly from Theorem 4.1 and continuous mapping theorem. As for $(b)$, note that $K$ is fixed, we thus only need to show the tightness of $(n-k)V_{n}^{*}(k)\Psi_{k}(\cdot)$ for each $1\leq k\leq K$, which is easily ensured if $\mathbb{E}^{*}[(n-k)V_{n}^{*}(k)]^{2}$ is finite. In fact,

[TABLE]

and using Minkowski inequality,

[TABLE]

Therefore, $\mathbb{E}\mathbb{E}^{*}[(n-k)V_{n}^{*}(k)]^{2}<\infty$, which implies that in probability, $\mathbb{E}^{*}[(n-k)V_{n}^{*}(k)]^{2}$ is finite. And thus, (i) is proved.

(ii) is trivially satisfied.

(iii) Note that under $\mathbb{P}^{*}$, $V_{n}^{*}(k_{1})$ is independent of $V_{n}^{*}(k_{2})$ for $k_{1}\neq k_{2}$. Therefore, by Chebyshev’s inequality,

[TABLE]

By similar arguments used in proving (i)(b) above, and that $\|\Psi_{k}\|^{2}\leq k^{-2}$, we have $\mathbb{P}^{*}(\|R_{n}^{K*}\|>\epsilon)=O_{p}(K^{-1})$. Therefore, letting $K\to\infty$ the results follows.

∎

Proof of 4.3 By 3.3, it suffices to show that $CvM^{*}_{n}=O^{*}_{p}(1)$ in probability, where $O_{p}^{*}(1)$ is analogous to $O_{p}(1)$ but for bootstrap sample asymptotics.

Note that $CvM^{*}_{n}>0$, hence by Markov inequality, we only need to show $\mathbb{E}\mathbb{E}^{*}(\textsc{CvM}_{n}^{*})=O(1)$. In fact, using the fact that $\mathbb{E}^{*}V_{n}^{*}(k)=0$, and $\mathrm{Cov}^{*}(V_{n}^{*}(k),V_{n}^{*}(k^{\prime}))=0$ for $k\neq k^{\prime}$, we have

[TABLE]

By (27), we have $\mathbb{E}\mathbb{E}^{*}(\textsc{CvM}_{n}^{*})\leq C\sum_{k=1}^{n-4}k^{-2}\leq C.$ This implies that $\mathbb{E}^{*}(\textsc{CvM}_{n}^{*})$ is bounded in probability, the result follows.

∎