Granger Causality Testing in High-Dimensional VARs: a Post-Double-Selection Procedure

TL;DR

This paper introduces a new Granger causality test for high-dimensional VAR models using a post-double-selection approach with penalized least squares, validated through simulations and applied to volatility spillover networks.

Contribution

It proposes a novel LM test for high-dimensional VARs that remains valid after variable selection, addressing a key challenge in high-dimensional time series analysis.

Findings

Test performs well across various data generating processes

Causal relationships are clearer in high-dimensional models

Method is effective even without sparsity assumptions

Abstract

We develop an LM test for Granger causality in high-dimensional VAR models based on penalized least squares estimations. To obtain a test retaining the appropriate size after the variable selection done by the lasso, we propose a post-double-selection procedure to partial out effects of nuisance variables and establish its uniform asymptotic validity. We conduct an extensive set of Monte-Carlo simulations that show our tests perform well under different data generating processes, even without sparsity. We apply our testing procedure to find networks of volatility spillovers and we find evidence that causal relationships become clearer in high-dimensional compared to standard low-dimensional VARs.

Tables

Table 1. Table 1: Simulation results for the PDS-LM Granger causality test

DGP	Size/Power	$\rho$	$T$			50					100					200					500
			K	AIC     BIC    EBIC     ${\lambda }^{th}$   ${\lambda }^{TSCV}$					AIC     BIC    EBIC     ${\lambda }^{th}$   ${\lambda }^{TSCV}$					AIC     BIC    EBIC     ${\lambda }^{th}$   ${\lambda }^{TSCV}$					AIC     BIC    EBIC     ${\lambda }^{th}$   ${\lambda }^{TSCV}$
			10	6.7	6.1	5.3	6.8	6.9	6.7	6.1	5.9	6.5	6.4	5.5	4.5	4.5	4.7	5.4	4.3	4.1	4.3	4.2	4.3
1	Size	0	20	7.3	6.0	6.4	5.9	6.2	6.0	4.7	4.7	5.5	5.8	4.0	4.7	5.4	4.3	4.5	4.3	4.1	4.3	4.1	4.3
			50	7.0	5.9	5.7	6.4	5.5	7.4	6.4	6.4	6.2	6.6	6.9	5.9	5.9	6.2	6.1	6.8	6.5	6.7	6.8	6.7
			100	7.4	7.2	6.6	NA	6.8	7.2	4.9	4.9	5.9	5.4	6.5	4.7	4.9	5.0	5.2	5.8	3.9	4.0	5.0	4.1
			10	30.0	31.2	33.3	30.4	30.9	58.3	58.9	60.7	58.5	57.4	89.1	89.3	89.6	89.2	89.5	99.9	99.9	99.9	99.9	99.9
1	Power	0	20	22.8	26.4	30.2	25.8	24.2	52.8	55.1	57.4	53.9	54.5	85.6	88.0	89.0	86.6	86.1	99.9	100	100	99.9	99.9
			50	13.6	24.3	33.3	17.9	18.7	38.7	53.8	59.0	46.4	45.5	78.2	85.8	87.0	81.3	80.2	99.9	99.9	99.9	99.9	99.9
			100	10.4	20.3	32.1	NA	18.7	20.4	51.9	57.0	31.2	35.8	61.6	85.4	86.7	74.1	73.7	99.5	100	100	99.8	99.7
			10	6.2	5.5	6.1	6.1	5.9	5.0	4.9	5.4	5.0	5.0	4.7	4.6	4.9	4.6	5.2	3.9	3.9	3.9	3.7	3.7
2	Size	0	20	7.8	5.5	6.0	5.8	6.5	5.4	4.6	4.7	4.8	5.3	5.9	5.9	5.4	5.7	6.0	4.7	3.8	4.3	4.6	4.7
			50	7.3	5.9	4.7	6.0	7.1	7.9	6.3	6.5	7.3	6.5	7.3	6.1	6.0	6.6	6.6	6.3	6.6	6.4	6.1	6.0
			100	6.1	6.9	5.5	NA	6.6	6.8	5.1	6.0	5.3	4.8	5.4	5.5	5.8	4.6	4.9	6.1	4.3	5.0	4.9	4.5
			10	18.0	19.7	21.3	18.2	17.0	37.9	39.3	40.4	38.0	38.4	64.7	64.7	66.9	64.8	65.1	97.4	97.4	97.5	97.4	97.4
2	Power	0	20	16.0	19.6	24.6	18.8	16.9	35.4	39.8	44.4	37.4	36.7	64.4	67.4	69.7	66.0	65.0	97.2	97.5	97.5	97.5	97.4
			50	8.6	15.2	21.7	12.8	13.1	25.0	36.1	43.2	32.6	31.2	57.0	66.4	71.9	61.8	61.1	95.0	96.2	96.8	96.1	95.4
			100	9.2	14.1	25.1	NA	10.1	15.1	34.9	45.9	25.9	25.8	44.8	65.1	74.7	58.0	57.8	94.7	97.3	97.7	96.3	96.5
			10	5.2	5.0	5.6	5.7	4.6	5.6	4.9	5.7	5.7	6.2	4.0	4.1	6.1	4.1	4.4	4.1	4.0	3.8	3.9	4.1
3	Size	0	20	4.2	5.1	5.7	5.6	4.9	4.3	4.3	7.2	4.2	3.9	5.2	5.6	9.4	4.7	4.6	4.7	4.4	4.5	4.5	4.8
			50	7.5	6.3	7.4	6.9	6.6	6.4	7.0	9.5	5.9	5.5	6.9	6.9	12.4	6.0	6.6	4.9	5.3	6.6	5.3	5.4
			100	7.1	6.7	8.3	NA	7.0	6.2	5.7	8.3	5.7	6.0	4.7	6.2	10.7	4.2	4.7	4.4	5.1	6.5	5.1	4.7
			10	15.4	20.0	23.4	16.2	16.1	31.5	36.4	44.3	32.6	31.4	58.4	61.3	63.7	58.9	59.3	95.2	95.6	95.7	95.5	95.3
3	Power	0	20	13.4	18.7	26.0	13.8	13.8	29.5	37.0	47.8	30.0	29.8	56.5	62.0	69.8	56.8	55.5	94.0	94.6	94.8	94.4	94.3
			50	11.4	20.7	28.0	12.9	10.6	24.1	39.8	52.3	26.9	27.6	50.3	59.5	73.6	52.1	51.8	91.1	92.7	93.4	92.3	90.9
			100	7.9	18.7	26.9	NA	14.1	13.7	42.6	55.0	20.2	22.2	41.4	62.0	75.2	44.9	44.8	90.0	94.0	94.8	91.4	90.2

Table 2. Table B.1: Simulation results for the PDS-LM Granger causality test ( ρ = 0.7 𝜌 0.7 \rho=0.7 )

DGP	Size/Power	$\rho$	$T$			50					100					200					500
			K	AIC     BIC    EBIC     ${\lambda }^{th}$   ${\lambda }^{TSCV}$					AIC     BIC    EBIC     ${\lambda }^{th}$   ${\lambda }^{TSCV}$					AIC     BIC    EBIC     ${\lambda }^{th}$   ${\lambda }^{TSCV}$					AIC     BIC    EBIC     ${\lambda }^{th}$   ${\lambda }^{TSCV}$
			10	5.8	5.1	5.4	6.1	5.7	5.8	5.7	5.2	5.2	5.8	5.1	4.6	4.6	5.0	4.5	4.3	4.0	4.2	4.0	4.2
1	Size	$0.7$	20	7.3	6.4	6.1	6.0	7.2	6.0	6.1	6.2	6.3	6.5	5.1	5.5	5.5	5.8	5.9	5.0	4.9	4.9	4.4	5.0
			50	7.6	5.5	6.0	7.5	6.8	6.4	7.3	5.8	7.3	6.7	6.5	4.9	4.2	6.2	5.5	5.7	4.8	4.8	5.2	5.3
			100	6.4	8.1	6.6	NA	8.0	7.7	6.4	6.2	7.1	6.8	6.9	5.8	4.4	5.3	4.6	5.6	3.6	3.6	4.2	3.8
			10	19.0	18.2	18.8	20.3	19.5	33.7	34.2	34.3	35.2	33.4	56.6	57.1	56.8	57.8	55.6	93.7	94.2	94.2	94.3	93.3
1	Power	$0.7$	20	15.3	18.3	19.1	18.5	17.9	29.2	30.7	31.6	30.7	29.8	54.4	56.8	56.6	56.3	54.8	92.8	94.1	94.1	93.8	93.4
			50	10.3	14.1	19.3	15.2	14.6	24.0	30.4	32.7	31.1	30.3	50.1	57.1	58.8	56.7	51.7	90.7	92.0	92.5	92.1	91.4
			100	9.0	14.0	19.2	NA	15.4	13.6	29.2	33.6	25.4	24.1	34.1	53.6	55.7	49.4	48.9	88.1	94.0	94.5	92.4	91.1
			10	6.4	5.8	5.5	5.8	6.8	5.0	5.1	5.3	5.0	5.5	5.2	4.9	5.0	4.5	5.1	4.7	5.0	4.9	4.5	5.1
2	Size	0.7	20	7.5	6.5	5.6	5.9	5.5	4.8	5.5	5.3	5.8	5.2	5.1	5.0	4.8	5.0	5.8	4.8	4.5	4.9	4.4	6.3
			50	6.8	5.9	6.8	6.0	6.2	7.4	5.3	6.1	5.2	6.3	5.8	5.9	5.7	5.2	6.2	5.6	5.5	5.6	5.9	6.5
			100	8.1	6.5	6.8	6.9	6.5	8.1	5.8	6.2	7.0	5.9	6.7	6.1	6.3	5.5	5.9	5.3	4.4	4.0	4.5	5.0
			10	15.4	15.1	16.0	14.5	14.3	26.0	28.1	28.8	26.6	26.5	47.4	49.7	51.6	49.2	51.0	83.2	84.6	85.3	84.6	86.2
2	Power	0.7	20	12.6	15.8	19.3	15.7	15.5	26.5	27.9	30.0	26.7	26.9	48.4	51.6	53.2	51.1	51.8	83.1	85.3	86.7	85.4	85.9
			50	9.6	15.4	18.4	13.0	12.2	19.1	28.5	31.5	26.1	26.3	40.3	50.5	52.2	47.0	48.9	83.9	87.4	88.8	87.2	87.2
			100	10.3	16.9	21.3	10.7	11.5	12.6	29.3	32.7	21.0	21.4	32.5	53.0	56.0	45.4	49.8	78.7	88.6	89.1	85.4	85.4
			10	4.8	4.8	4.7	4.3	5.5	5.3	4.6	4.8	4.6	5.1	5.5	5.5	5.6	5.6	5.4	4.8	4.8	5.0	4.8	5.1
3	Size	0.7	20	5.6	5.4	5.2	5.1	5.6	5.4	5.8	5.6	5.2	5.2	4.4	4.0	4.0	4.1	4.5	4.5	4.0	4.2	4.2	4.4
			50	8.3	6.0	5.3	5.5	7.0	6.1	5.9	5.2	5.9	6.4	5.6	4.8	4.6	5.2	5.1	6.2	6.8	6.7	6.9	6.7
			100	7.6	6.3	6.3	5.6	6.8	7.2	5.2	5.2	6.1	6.4	5.6	5.8	5.8	6.2	4.7	3.6	4.1	3.7	4.4	4.6
			10	9.9	9.8	10.3	10.3	10.7	16.2	16.6	16.4	16.8	16.9	31.5	31.3	31.4	31.6	31.5	68.4	68.9	68.9	68.9	68.8
3	Power	0.7	20	8.1	8.8	8.2	9.1	9.5	15.7	16.4	16.2	16.4	16.8	29.2	29.0	29.0	28.7	28.3	66.5	68.5	69.4	68.8	67.9
			50	9.7	9.4	10.1	9.7	10.2	15.7	16.1	17.4	16.2	16.6	31.3	31.2	31.7	31.8	29.9	65.1	66.6	67.4	66.9	65.0
			100	8.4	9.4	9.3	9.7	9.9	12.0	15.3	17.7	14.6	15.6	26.4	30.9	31.9	30.2	29.3	66.2	67.8	68.3	68.3	66.5

Table 3. Table B.2: Simulation results for the PDS-WALD Granger causality test

DGP	Size/Power	$\rho$	$T$			50					100					200					500
			K	AIC     BIC    EBIC     ${\lambda }^{th}$   ${\lambda }^{TSCV}$					AIC     BIC    EBIC     ${\lambda }^{th}$   ${\lambda }^{TSCV}$					AIC     BIC    EBIC     ${\lambda }^{th}$   ${\lambda }^{TSCV}$					AIC     BIC    EBIC     ${\lambda }^{th}$   ${\lambda }^{TSCV}$
			10	7.1	6.6	6.0	7.0	7.6	7.0	6.4	6.1	6.8	7.1	5.5	4.3	4.6	4.7	6.9	4.3	4.1	4.3	4.3	4.8
1	Size	0	20	8.0	6.4	6.8	6.5	6.4	6.2	4.9	4.9	5.7	6.9	4.2	4.8	5.5	4.4	5.0	4.3	4.1	4.3	4.1	4.7
			50	7.2	6.2	6.2	7.2	6.2	7.6	6.4	6.6	6.4	7.1	7.1	5.9	5.9	6.3	7.0	6.8	6.5	6.7	6.8	7.1
			100	7.1	7.4	7.3	NA	7.3	8.0	5.2	5.1	6.1	5.9	6.6	4.7	4.9	5.0	5.3	5.8	3.9	4.0	5.2	4.5
			10	31.3	31.8	34.0	31.3	33.8	58.6	59.5	61.2	58.9	58.3	89.1	89.3	89.6	89.3	89.6	99.9	99.9	99.9	99.9	99.9
1	Power	0	20	23.9	27.6	30.9	27.3	25.3	53.4	55.9	58.4	54.6	55.0	85.6	88.2	89.1	86.8	86.1	99.9	100	100	99.9	99.8
			50	15.3	26.0	34.9	19.0	20.4	39.9	54.5	60.0	46.7	46.3	78.8	86.1	87.3	81.4	80.3	99.9	99.9	99.9	99.9	99.9
			100	13.0	21.9	33.1	NA	20.0	21.6	52.6	57.3	31.3	36.7	61.7	85.9	87.3	74.2	74.2	99.5	100	100	99.8	99.7
			10	6.5	6.1	6.5	6.4	6.6	5.1	5.1	5.8	5.1	6.2	5.1	4.8	5.2	4.8	6.3	3.9	3.9	3.9	3.7	4.2
2	Size	0	20	8.4	6.0	6.7	6.5	7.4	5.3	4.8	5.1	5.2	6.0	5.9	6.0	5.8	5.8	6.9	4.7	3.8	4.3	4.6	5.1
			50	8.3	6.2	5.1	6.7	7.6	8.4	6.4	6.7	7.4	7.4	7.6	6.4	6.2	6.8	7.7	6.4	6.6	6.4	6.2	6.4
			100	6.5	7.6	5.9	NA	7.2	7.4	5.3	6.2	5.6	5.2	5.4	5.6	5.9	4.7	5.1	6.1	4.3	5.1	5.0	4.8
			10	18.7	20.2	21.6	18.9	19.2	38.3	39.6	40.9	38.7	40.3	65.2	65.2	67.4	65.0	65.5	97.4	97.4	97.6	97.3	97.6
2	Power	0	20	16.7	20.9	25.4	19.7	19.0	35.8	40.2	44.9	38.1	38.1	64.8	67.5	69.8	66.2	65.5	97.2	97.5	97.4	97.5	97.4
			50	10.1	15.8	22.7	13.8	14.4	25.4	36.8	43.8	33.6	32.1	57.2	66.8	72.1	62.1	61.6	95.0	96.2	96.8	96.1	95.4
			100	10.0	14.6	25.9	NA	11.4	16.6	35.1	46.5	27.3	26.2	45.1	65.3	74.9	57.3	57.9	94.7	97.3	97.7	96.3	96.4
			10	5.5	5.6	6.1	5.8	5.2	6.0	5.1	5.9	5.9	6.3	3.9	4.1	6.1	4.2	4.6	4.1	4.1	3.9	4.0	4.2
3	Size	0	20	4.8	5.5	6.3	5.3	5.5	4.7	4.4	7.4	4.3	4.3	5.4	5.6	9.6	4.7	4.8	4.7	4.4	4.3	4.6	4.8
			50	7.9	6.7	7.8	7.4	7.0	6.7	7.2	9.6	6.1	5.6	7.0	6.9	12.7	6.1	6.8	5.0	5.3	6.6	5.3	5.4
			100	7.4	7.0	8.4	NA	7.7	7.0	6.0	8.6	5.9	6.0	4.7	6.3	10.8	4.2	4.8	4.4	5.1	6.7	5.1	4.8
			10	16.0	20.7	24.5	16.6	17.4	32.1	36.8	44.6	32.9	31.9	58.6	61.4	63.8	59.1	59.6	95.2	95.6	95.7	95.5	95.3
3	Power	0	20	14.3	19.9	27.2	14.2	14.8	29.9	37.8	48.8	30.6	30.3	56.7	62.2	70.0	57.2	55.7	94.1	94.6	94.8	94.5	94.3
			50	12.6	21.7	28.8	13.7	11.8	24.5	40.5	52.8	27.3	27.9	50.6	59.7	73.8	52.5	51.8	91.2	92.8	93.4	92.4	90.9
			100	9.8	20.0	27.7	NA	14.8	15.4	43.0	55.5	20.3	23.0	41.7	62.2	75.2	45.3	45.0	90.0	94.2	95.0	91.5	90.2
			10	6.2	5.5	5.9	6.4	5.6	5.9	5.9	5.4	5.4	5.2	5.1	4.7	4.7	5.1	4.4	4.5	4.1	4.2	4.1	4.2
1	Size	$0.7$	20	7.8	6.4	6.2	6.4	7.1	6.5	6.2	6.2	6.4	6.3	5.1	5.5	5.5	5.8	5.7	5.1	4.9	4.9	4.6	5.1
			50	9.0	6.4	6.9	8.3	6.4	6.9	7.5	6.1	7.5	6.6	6.6	5.1	4.3	6.5	5.5	5.7	4.8	4.9	5.2	5.3
			100	8.2	8.8	6.9	7.5	8.1	8.3	6.5	6.4	7.4	6.5	7.0	5.8	4.6	5.3	4.6	5.6	3.6	3.7	4.2	4.7
			10	19.7	18.9	19.8	21.1	19.3	34.7	35.2	35.2	35.9	33.9	56.7	57.3	57.1	57.8	55.7	93.7	94.2	94.2	94.3	93.9
1	Power	$0.7$	20	16.2	19.3	20.0	19.0	17.9	29.7	31.0	32.2	31.5	29.8	54.7	57.1	56.8	56.5	55.3	92.8	94.2	94.2	93.5	93.4
			50	11.3	14.9	20.1	16.4	15.1	24.9	30.8	33.7	31.9	30.8	50.6	57.3	59.0	56.9	52.2	90.7	92.0	92.5	92.1	91.4
			100	9.6	14.9	20.3	12.4	15.4	15.4	29.6	34.0	25.8	24.6	34.5	53.8	56.0	49.7	49.3	88.2	94.1	94.6	92.4	91.2

Table 4. Table B.3: Simulation results for the PDS-LM Granger causality test (Overspecified lag-length)

DGP	Size/Power	$\rho$	$T$			50					100					200					500
			K	AIC     BIC    EBIC     ${\lambda }^{th}$   ${\lambda }^{TSCV}$					AIC     BIC    EBIC     ${\lambda }^{th}$   ${\lambda }^{TSCV}$					AIC     BIC    EBIC     ${\lambda }^{th}$   ${\lambda }^{TSCV}$					AIC     BIC    EBIC     ${\lambda }^{th}$   ${\lambda }^{TSCV}$
1	Size	0	10	7.0	6.4	6.3	6.7	7.2	5.9	6.3	5.8	6.7	5.7	5.6	5.0	5.1	5.2	5.4	5.3	5.4	5.4	5.2	5.6
			20	8.4	7.5	6.7	8.3	7.3	6.3	6.0	5.6	4.9	6.0	5.8	5.6	5.6	5.7	5.7	5.0	3.7	3.7	4.7	4.3
			50	NA	5.4	4.6	NA	NA	8.7	5.0	5.4	6.7	5.9	7.4	5.2	4.9	6.6	6.1	5.7	4.7	5.0	4.7	4.9
			100	NA	NA	5.2	NA	NA	NA	5.0	5.2	NA	NA	7.0	4.3	4.0	5.9	5.5	5.9	5.1	5.3	5.6	4.3
1	Power	0	10	20.2	23.3	24.7	22.5	21.0	45.4	48.7	50.0	47.2	44.8	82.3	83.2	83.8	82.6	81.5	99.7	99.8	99.8	99.8	99.7
			20	13.8	17.4	22.6	16.5	15.1	33.0	44.3	47.5	38.2	38.1	72.9	78.9	79.4	75.2	74.4	99.6	99.6	99.6	99.6	99.6
			50	NA	NA	23.8	NA	NA	16.6	43.0	47.8	29.9	29.6	55.9	77.3	79.8	68.2	66.9	98.0	99.7	99.7	99.0	98.9
			100	NA	NA	23.9	NA	NA	NA	37.2	45.5	NA	NA	20.0	76.3	79.1	53.7	53.1	96.0	99.8	99.8	99.2	94.2
2	Size	0	10	5.4	5.0	4.3	5.2	6.0	6.8	6.9	6.4	6.8	5.9	5.8	5.3	6.0	5.4	5.8	5.1	5.1	5.0	5.1	4.8
			20	7.3	5.6	5.7	7.2	6.6	5.1	4.5	5.0	5.0	5.5	4.9	5.0	5.2	5.6	5.2	4.4	4.4	4.7	4.2	4.1
			50	NA	5.4	5.3	NA	NA	6.4	5.3	6.3	4.7	4.9	8.3	5.2	5.8	5.2	6.1	6.3	5.2	4.8	4.2	4.5
			100	NA	NA	5.7	NA	NA	NA	4.6	5.4	NA	4.6	7.6	5.0	5.9	5.0	6.5	6.1	5.1	4.9	4.5	5.4
2	Power	0	10	13.1	14.7	14.8	14.5	13.6	29.0	29.9	31.3	29.1	28.7	52.7	52.5	55.0	53.5	53.9	94.1	94.3	94.4	94.5	93.8
			20	10.8	15.7	18.7	15.4	13.7	24.0	28.0	34.2	26.8	26.2	52.5	54.9	58.2	54.3	51.7	92.6	92.7	93.0	92.8	91.9
			50	NA	NA	14.9	NA	NA	9.2	26.5	33.7	23.5	20.5	38.4	53.0	61.0	50.0	46.5	88.7	92.1	93.5	91.6	90.3
			100	NA	NA	20.3	NA	NA	NA	26.9	37.9	NA	5.7	13.4	54.8	62.7	42.3	13.6	83.9	93.6	95.2	91.4	89.2
3	Size	0	10	5.6	5.9	7.9	5.8	4.8	6.0	6.4	9.1	6.0	6.1	4.7	5.1	8.2	4.9	4.9	5.1	5.3	4.9	5.1	5.0
			20	6.4	5.5	7.7	6.3	5.8	5.3	5.7	8.2	3.8	4.6	5.8	5.3	9.0	5.1	5.6	4.2	3.5	3.6	3.5	3.7
			50	NA	6.8	7.7	NA	NA	7.9	6.8	11.8	5.9	5.0	6.2	6.2	13.1	5.2	5.6	4.8	4.5	5.5	4.6	5.1
			100	NA	NA	7.7	NA	NA	NA	6.9	11.4	NA	6.3	7.2	6.3	12.2	5.4	5.4	3.7	5.0	5.5	4.3	3.2
3	Power	0	10	10.0	12.2	16.5	9.8	10.1	22.1	24.8	31.7	22.0	21.7	43.2	44.0	50.0	43.8	43.4	87.4	87.4	87.8	87.3	86.8
			20	8.5	13.4	19.9	9.7	9.4	15.5	23.2	37.0	18.4	18.4	39.8	43.4	53.6	41.2	39.9	85.3	86.8	87.0	86.3	84.9
			50	NA	NA	20.5	NA	NA	10.4	26.3	40.1	18.4	17.9	32.2	45.8	59.0	38.3	35.5	78.3	81.6	82.6	79.9	79.4
			100	NA	NA	20.8	NA	NA	NA	26.4	40.2	NA	NA	14.1	46.5	63.8	27.3	30.0	73.0	83.9	85.6	80.9	77.2
1	Size	0.7	10	7.6	6.9	5.9	6.6	6.5	5.5	5.5	5.7	4.7	6.2	4.7	4.6	5.1	4.9	5.4	4.8	4.8	4.8	4.5	4.5
			20	6.9	7.3	7.4	7.7	6.9	6.3	7.1	5.8	7.2	6.0	5.1	4.7	4.9	5.6	5.6	5.4	5.6	5.3	6.0	5.1
			50	NA	NA	7.2	NA	NA	8.6	6.2	5.5	7.2	6.3	7.8	4.6	4.4	5.9	5.0	6.5	5.3	4.9	5.6	5.8
			100	NA	NA	6.5	NA	NA	NA	6.1	5.8	NA	5.2	6.3	5.1	4.5	4.5	5.3	5.9	4.3	4.6	3.9	4.3
1	Power	0.7	10	13.0	13.6	14.4	14.4	14.2	22.7	24.0	25.3	24.7	24.0	42.1	43.4	43.9	43.4	42.7	86.0	87.5	87.5	87.6	86.6
			20	10.2	11.5	12.5	12.8	13.2	19.6	21.1	20.9	21.1	20.4	39.4	42.3	43.2	42.0	40.9	84.8	86.8	86.7	86.4	85.0
			50	NA	NA	15.0	NA	NA	12.6	22.3	23.6	21.6	20.2	30.6	43.6	46.4	41.8	38.8	80.7	85.8	86.0	84.4	82.8
			100	NA	NA	13.2	NA	NA	NA	19.9	23.1	NA	NA	14.3	41.3	44.2	35.0	34.7	72.3	85.2	86.1	82.8	81.5

Table 5. Table C.1: Stocks used in Section 6

N.	Symbol	Issue name	N.	Symbol	Issue name
1	AAPL	APPLE INC	16	KO	COCA-COLA CO
2	AXP	AMERICAN EXPRESS CO	17	MCD	MCDONALD’S CORP
3	BA	BOEING CO	18	MMM	3M
4	CAT	CATERPILLAR	19	MRK	MERCK & CO
5	CSCO	CISCO SYSTEMS	20	MSFT	MICROSOFT CORPORATION
6	CVX	CHEVRON CORP	21	NKE	NIKE INC
7	DD	DOW CHEMICAL COMPANY	22	PFE	PFIZER INC
8	DIS	WALT DISNEY CO	23	PG	PROCTER & GAMBLE CO
9	GE	GENERAL ELEC	24	TRV	TRAVELERS COMPANIES INC
10	GS	GOLDMAN SACHS GROUP INC	25	UNH	UNITEDHEALTH GROUP INC
11	HD	HOME DEPOT INC	26	UTX	UNITED TECHNOLOGIES CORPORATION
12	IBM	INTL BUS MACHINE	27	V	VISA INC
13	INTC	INTEL CORP	28	VZ	VERIZON COMMUNICATIONS INC
14	JNJ	JOHNSON &JOHNSON	29	WMT	WALMART INC
15	JPM	JPMORGAN CHASE & CO	30	XOM	EXXON MOBIL CORPORATION

Equations

$${\mathbf{y}}_t={\mathbf{A}}_1{\mathbf{y}}_{t-1}+\cdots +{\mathbf{A}}_p{y}_{t-p}+{\mathbf{u}}_t,t=p+1,\dots ,T,$$

$$\mathcal{P}({\mathbf{y}}_{I,t}|{\mathbf{x}}_{J^c,t}^{(p)})=\mathcal{P}({\mathbf{y}}_{I,t}|{\mathbf{x}}_{t-1}^{(p)})$$

$$\left[\begin{array}{c}{\mathbf{y}}_{I,t}\\ {\mathbf{y}}_{J,t}\\ {\mathbf{y}}_{-(I\cup J),t}\end{array}\right]=\left[\begin{array}{ccc}{\mathbf{A}}_{I,I}& {\mathbf{A}}_{I,J}& {\mathbf{A}}_{I,-(I\cup J)}\\ {\mathbf{A}}_{J,I}& {\mathbf{A}}_{J,J}& {\mathbf{A}}_{J,-(I\cup J)}\\ {\mathbf{A}}_{I,-(I\cup J)}& {\mathbf{A}}_{-(I\cup J)}& {\mathbf{A}}_{-(I\cup J),-(I\cup J)}\end{array}\right]\left[\begin{array}{c}{\mathbf{y}}_{I,t-1}\\ {\mathbf{y}}_{J,t-1}\\ {\mathbf{y}}_{-(J\cup I),t-1}\end{array}\right]+{\mathbf{u}}_t,$$

$$\begin{array}{rl}& \mathcal{P}({\mathbf{y}}_{I,t}|{\mathbf{y}}_{t-1})={\mathbf{A}}_{I,I}{\mathbf{y}}_{I,t-1}+{\mathbf{A}}_{I,J}{\mathbf{y}}_{J,t-1}+{\mathbf{A}}_{I,-(I\cup J)}{\mathbf{y}}_{-(I\cup J),t-1},\\ & \mathcal{P}({\mathbf{y}}_{I,t}|{\mathbf{y}}_{J^c,t-1})={\mathbf{A}}_I^{\ast }{\mathbf{y}}_{J^c,t-1},\text{ where }{\mathbf{A}}_I^{\ast }=\underset{{\mathbf{A}}_I}{\min }\mathbb{E}\left[{\Vert {\mathbf{y}}_{I,t}-{\mathbf{A}}_I{\mathbf{y}}_{J^c,t-1}\Vert }_2^2\right].\end{array}$$

$${\mathbf{y}}_I={\mathbf{X}}^{\otimes }\mathbf{\beta }+{\mathbf{u}}_I={\mathbf{X}}_{GC}^{\otimes }{\mathbf{\beta }}_{GC}+{\mathbf{X}}_{-GC}^{\otimes }{\mathbf{\beta }}_{-GC}+{\mathbf{u}}_I,$$

$$\hat{\mathbf{\beta }}(\lambda )=\underset{\mathbf{\beta }}{\mathrm{argmin}}\left(\frac{1}{T}{\Vert \mathbf{y}-\mathbf{X}\mathbf{\beta }\Vert }_2^2+\lambda \sum _{m=1}^M\left|w_m{\beta }_m\right|\right),$$

$${\mathbf{y}}_I$$

$${\mathbf{x}}_{GC,j}$$

$$\begin{array}{rl}{\mathbf{\gamma }}_0& =\underset{\mathbf{\gamma }}{\mathrm{argmin}}\mathbb{E}{\Vert {\mathbf{y}}_{I,t}-{\mathbf{X}}_{-GC,t-1}^{\otimes \prime }\mathbf{\gamma }\Vert }_2^2={\left(\mathbb{E}{\mathbf{X}}_{-GC,t-1}^{\otimes }{\mathbf{X}}_{-GC,t-1}^{\otimes \prime }\right)}^{-1}\mathbb{E}{\mathbf{X}}_{-GC,t-1}^{\otimes }{\mathbf{y}}_{i,t},\\ {\mathbf{\gamma }}_j& =\underset{\mathbf{\gamma }}{\mathrm{argmin}}\mathbb{E}{\Vert x_{GC,j,t}-{\mathbf{x}}_{-GC,t-1}^{\prime }\mathbf{\gamma }\Vert }_2^2\\ & ={\left(\mathbb{E}{\mathbf{x}}_{-GC,t-1}{\mathbf{x}}_{-GC,t-1}^{\prime }\right)}^{-1}\mathbb{E}{\mathbf{x}}_{-GC,t-1}{\mathbf{x}}_{GC,j,t},j=1,\dots ,N_X,\end{array}$$

$${\lambda }^{IC}=\underset{\lambda }{\mathrm{argmin}}\ln \left(\frac{1}{T}{\Vert \mathbf{y}-\mathbf{X}\hat{\mathbf{\beta }}(\lambda )\Vert }_2^2\right)+\frac{C_T}{T}\left|\hat{S}(\lambda )\right|,$$

$$\Sigma$$

$$\begin{array}{rl}{\mathbf{E}}^{\otimes \prime }{\mathbf{u}}_I/\sqrt{T}& =\mathrm{vec}({\mathbf{E}}^{\prime }{\mathbf{U}}_I)/\sqrt{T}=\frac{1}{\sqrt{T}}\sum _{t=p+1}^T\mathrm{vec}({\mathbf{e}}_t{\mathbf{u}}_{I,t}^{\prime })\stackrel{d}{\to }N(0,\Omega ),\\ {\mathbf{E}}^{\prime }\mathbf{E}/T& =\frac{1}{T}\sum _{t=p+1}^T{\mathbf{e}}_t{\mathbf{e}}_t^{\prime }\stackrel{p}{\to }{\Sigma }_{GC|-GC}={\Sigma }_{GC,GC}-{\Sigma }_{GC,-GC}{\Sigma }_{-GC,-GC}^{-1}{\Sigma }_{-GC,GC},\\ {\mathbf{U}}_I^{\prime }{\mathbf{U}}_I/T& \stackrel{p}{\to }{\Sigma }_{u,I},\end{array}$$

$$\begin{split}&{\left\lVert\bm{X}_{-GC}^{\otimes}(\hat{\bm{\gamma}}_{0}-\bm{\gamma}_{0})\right\rVert}_{2}/\sqrt{T}\leq\delta_{T}T^{-1/4},&\max_{j=1,\ldots,N_{X}}{\left\lVert\bm{X}_{-GC}(\hat{\bm{\gamma}}_{j}-\bm{\gamma}_{j})\right\rVert}_{2}/\sqrt{T}\leq\delta_{T}T^{-1/4}.\end{split}$$

$${\tilde{\mathbf{y}}}_I={\tilde{\mathbf{X}}}_{GC}^{\otimes }{\mathbf{\beta }}_{GC}^{\text{textsc}pds}+{\tilde{\mathbf{X}}}_{{\hat{S}}^{\otimes }}^{\otimes }{\mathbf{\beta }}_{{\hat{S}}^{\otimes }}^{\text{textsc}pds}+{\tilde{\mathbf{u}}}_I.$$

$$\sqrt{T}({\hat{\mathbf{\beta }}}_{GC}^{\text{textsc}pds}-{\mathbf{\beta }}_{GC})\stackrel{d}{\to }\mathcal{N}\left(0,({\mathbf{G}}^{\prime }\mathbf{G}\otimes {\Sigma }_{GC|-GC}{)}^{-1}{\Omega }_{\mathbf{G}}({\mathbf{G}}^{\prime }\mathbf{G}\otimes {\Sigma }_{GC|-GC}{)}^{-1}\right),$$

$$LM\stackrel{d}{\to }{\chi }_{N_{GC}}^2\text{as }T\to \infty .$$

$${\mathbf{y}}_t=\left[\begin{array}{ccc}0.5& \dots & 0\\ ⋮& \ddots & ⋮\\ 0& \dots & 0.5\end{array}\right]{\mathbf{y}}_{t-1}+{\mathbf{ϵ}}_t,$$

$${\mathbf{y}}_t=\left[\begin{array}{ccc}(-1{)}^{|i-j|}{a}^{|i-j|+1}& \dots & (-1{)}^{|i-j|}{a}^{|i-j|+1}\\ ⋮& \ddots & ⋮\\ (-1{)}^{|i-j|}{a}^{|i-j|+1}& \dots & (-1{)}^{|i-j|}{a}^{|i-j|+1}\end{array}\right]{\mathbf{y}}_{t-1}+{\mathbf{ϵ}}_t,\text{with}a=0.4,$$

$${\mathbf{y}}_t=\left[\begin{array}{ccc}\mathbf{A}& \dots & 0\\ ⋮& \ddots & ⋮\\ 0& \dots & \mathbf{A}\end{array}\right]{\mathbf{y}}_{t-1}+{\mathbf{ϵ}}_t\text{with}\underset{5\times 5}{\underbrace{\mathbf{A}}}=\left[\begin{array}{ccc}0.15& \cdots & 0.15\\ ⋮& \ddots & ⋮\\ 0.15& \cdots & 0.15\end{array}\right].$$

$$y_{2,t}={\beta }_{GC}{y}_{1,t-1}+\sum _{j=2}^K{\beta }_j{y}_{j,t-1}+{ϵ}_{2,t}.$$

$$RV10_t\equiv \sum _{j=1}^M{r}_{j,t}^2,r_{j,t}=\ln P_{j,t}-\ln P_{j-1,t},$$

$${\mathbf{y}}_t=\mathbf{c}+{\mathbf{B}}^{(1)}{\mathbf{y}}_{t-1}+{\mathbf{B}}^{(2)}{\mathbf{y}}_{t-1}^{(week)}+{\mathbf{B}}^{(3)}{\mathbf{y}}_{t-1}^{(month)}+{\mathbf{ϵ}}_t,$$

$$H_0:{\beta }_{i,k}^{(1)}={\beta }_{i,k}^{(2)}={\beta }_{i,k}^{(3)}=0\text{vs.}{H}_1:{\beta }_{i,k}^{(1)},{\beta }_{i,k}^{(2)},{\beta }_{i,k}^{(3)}\ne 0,$$

$${\Vert \mathbf{\delta }\Vert }_2^2\le {\bar{s}}_T{\Vert {\tilde{\mathbf{X}}}_{-GC}^{\otimes }\mathbf{\delta }/\sqrt{T}\Vert }_2^2/{\varphi }_{T,\min }^2,$$

$$\begin{array}{rl}& \underset{\mathbf{\delta }:|S_{\mathbf{\delta }}|\le {\bar{s}}_T}{\min }{\Vert {\tilde{\mathbf{X}}}_{-GC}^{\otimes }\mathbf{\delta }/\sqrt{T}\Vert }_2^2/{\Vert \mathbf{\delta }\Vert }_2^2=\underset{\mathbf{\delta }:|S_{\mathbf{\delta }}|\le {\bar{s}}_T}{\min }{\mathbf{\delta }}^{\prime }\left({\mathbf{G}}_T^{\prime }{\mathbf{G}}_T\otimes {\mathbf{X}}_{-GC}^{\prime }{\mathbf{X}}_{-GC}/T\right)\mathbf{\delta }/{\Vert \mathbf{\delta }\Vert }_2^2\\ & =\underset{|S|\le {\bar{s}}_T}{\min }\underset{\mathbf{x}}{\min }{\mathbf{x}}^{\prime }\left({\mathbf{G}}_T^{\prime }{\mathbf{G}}_T\otimes {\mathbf{X}}_S^{\prime }{\mathbf{X}}_S/T\right)\mathbf{x}/{\Vert \mathbf{x}\Vert }_2^2=\underset{|S|\le {\bar{s}}_T}{\min }{\lambda }_{\min }({\mathbf{G}}_T^{\prime }{\mathbf{G}}_T\otimes {\mathbf{X}}_S^{\prime }{\mathbf{X}}_S/T)\\ & ={\lambda }_{\min }({\mathbf{G}}_T^{\prime }{\mathbf{G}}_T)\underset{|S|\le {\bar{s}}_T}{\min }{\lambda }_{\min }({\mathbf{X}}_S^{\prime }{\mathbf{X}}_S/T)\ge C{\varphi }_{T,\min },\end{array}$$

$$\begin{array}{rl}\sqrt{T}({\hat{\mathbf{\beta }}}_{GC}^{\text{textsc}pds}-{\mathbf{\beta }}_{GC})& =\underset{{\mathbf{B}}_T^{-1}}{\underbrace{{\left({\tilde{\mathbf{X}}}_{GC}^{\otimes }\mathcal{M}({\tilde{\mathbf{X}}}_{{\hat{S}}^{\otimes }}^{\otimes }){\tilde{\mathbf{X}}}_{GC}^{\otimes }/T\right)}^{-1}}}\\ & \times \underset{{\mathbf{a}}_T}{\underbrace{{\tilde{\mathbf{X}}}_{GC}^{\otimes }\mathcal{M}({\tilde{\mathbf{X}}}_{{\hat{S}}^{\otimes }})\left[{\tilde{\mathbf{X}}}_{-GC}^{\otimes }{\mathbf{\beta }}_{-GC}+{\tilde{\mathbf{u}}}_I\right]/\sqrt{T}}}\end{array}$$

$$\begin{array}{rl}{\mathbf{a}}_T& ={\tilde{\mathbf{E}}}^{\otimes \prime }{\tilde{\mathbf{u}}}_I/\sqrt{T}+\underset{{\mathbf{a}}_{T,1}}{\underbrace{{\Gamma }^{\otimes \prime }{\tilde{\mathbf{X}}}_{-GC}^{\otimes \prime }\mathcal{M}({\tilde{\mathbf{X}}}_{{\hat{S}}^{\otimes }}^{\otimes }){\tilde{\mathbf{X}}}_{-GC}^{\otimes }{\mathbf{\beta }}_{-GC}/\sqrt{T}}}+\underset{{\mathbf{a}}_{T,2}}{\underbrace{{\Gamma }^{\otimes \prime }{\tilde{\mathbf{X}}}_{-GC}^{\otimes \prime }\mathcal{M}({\tilde{\mathbf{X}}}_{{\hat{S}}^{\otimes }}^{\otimes }){\tilde{\mathbf{u}}}_I/\sqrt{T}}}\\ & +\underset{{\mathbf{a}}_{T,3}}{\underbrace{{\tilde{\mathbf{E}}}^{\otimes \prime }\mathcal{M}({\tilde{\mathbf{X}}}_{{\hat{S}}^{\otimes }}^{\otimes }){\tilde{\mathbf{X}}}_{-GC}^{\otimes }{\mathbf{\beta }}_{-GC}/\sqrt{T}}}-\underset{{\mathbf{a}}_{T,4}}{\underbrace{{\tilde{\mathbf{E}}}^{\otimes \prime }\mathcal{P}({\tilde{\mathbf{X}}}_{{\hat{S}}^{\otimes }}^{\otimes }){\tilde{\mathbf{u}}}_I/\sqrt{T}}}\end{array}$$

$${\Vert {\mathbf{a}}_{T,1}\Vert }_2\le \sqrt{T}\underset{{\Vert {\mathbf{A}}_{T,1,1}\Vert }_2}{\underbrace{{\Vert \mathcal{M}({\tilde{\mathbf{X}}}_{{\hat{S}}^{\otimes }}^{\otimes }){\tilde{\mathbf{X}}}_{-GC}^{\otimes }{\Gamma }^{\otimes }/\sqrt{T}\Vert }_2}}\underset{{\Vert {\mathbf{a}}_{T,1,2}\Vert }_2}{\underbrace{{\Vert \mathcal{M}({\tilde{\mathbf{X}}}_{{\hat{S}}^{\otimes }}^{\otimes }){\tilde{\mathbf{X}}}_{-GC}^{\otimes }{\mathbf{\beta }}_{-GC}/\sqrt{T}\Vert }_2}},$$

$$\begin{array}{rl}& \mathcal{M}\left({\mathbf{G}}_T^{\otimes }{\mathbf{X}}_{{\hat{S}}_X^{\otimes }}^{\otimes }\right)=\mathcal{M}\left({\mathbf{G}}_T\otimes {\mathbf{X}}_{{\hat{S}}_X}\right){\tilde{\mathbf{X}}}_{-GC}^{\otimes }{\Gamma }^{\otimes }\\ & =\mathcal{M}\left({\mathbf{G}}_T\otimes {\mathbf{X}}_{{\hat{S}}_X}\right)\left({\mathbf{G}}_T\otimes {\mathbf{X}}_{-GC}\Gamma \right)={\mathbf{G}}_T\otimes \mathcal{M}\left({\mathbf{X}}_{{\hat{S}}_X}\right){\mathbf{X}}_{-GC}\Gamma ,\end{array}$$

$$\begin{array}{rl}{\Vert {\mathbf{A}}_{T,1,1}\Vert }_2& \le {\Vert \mathcal{M}\left({\tilde{\mathbf{X}}}_{{\hat{S}}_X^{\otimes }}^{\otimes }\right){\tilde{\mathbf{X}}}_{-GC}^{\otimes }{\Gamma }^{\otimes }/\sqrt{T}\Vert }_2\le {\Vert {\mathbf{G}}_T\Vert }_2{\Vert \mathcal{M}({\mathbf{X}}_{{\hat{S}}_X}){\mathbf{X}}_{-GC}\Gamma /\sqrt{T}\Vert }_2\\ & \le {\Vert {\mathbf{G}}_T\Vert }_2\sum _{j=1}^{N_X}{\Vert \mathcal{M}({\mathbf{X}}_{{\hat{S}}_X}){\mathbf{X}}_{-GC}{\mathbf{\gamma }}_j/\sqrt{T}\Vert }_2\le {\Vert {\mathbf{G}}_T\Vert }_2\sum _{j=1}^{N_X}{\Vert \mathcal{M}({\mathbf{X}}_{\widehat{S_j}}){\mathbf{X}}_{-GC}{\mathbf{\gamma }}_j/\sqrt{T}\Vert }_2.\end{array}$$

$$\begin{array}{rl}{\Vert \mathcal{M}({\mathbf{X}}_{{\hat{S}}_j}){\mathbf{X}}_{-GC}{\mathbf{\gamma }}_j/\sqrt{T}\Vert }_2& =\underset{\gamma :{\gamma }_m=0,m\notin {\hat{S}}_j}{\min }{\Vert {\mathbf{X}}_{-GC}{\mathbf{\gamma }}_j-{\mathbf{X}}_{\widehat{S_j}}\mathbf{\gamma }\Vert }_2/\sqrt{T}\\ & \le {\Vert {\mathbf{X}}_{-GC}({\mathbf{\gamma }}_j-{\hat{\mathbf{\gamma }}}_j)\Vert }_2/\sqrt{T},j=0,\dots ,N_X,\end{array}$$

$$\begin{array}{rl}{\mathbf{\gamma }}_0& ={\left(\mathbb{E}{\mathbf{X}}_{-GC}^{\otimes \prime }{\mathbf{X}}_{-GC}^{\otimes }\right)}^{-1}\mathbb{E}{\mathbf{X}}_{-GC}^{\prime }({\mathbf{X}}_{GC}^{\otimes }{\mathbf{\beta }}_{GC}+{\mathbf{X}}_{-GC}^{\otimes }{\mathbf{\beta }}_{-GC}+{\mathbf{u}}_I)={\Gamma }^{\otimes }{\mathbf{\beta }}_{GC}+{\mathbf{\beta }}_{-GC},\end{array}$$

Full text

Keywords: Granger causality, Post-double-selection, vector autoregressive models, high-dimensional inference.

JEL codes: C55, C12, C32.

1 Introduction

Economics, statistics and finance have seen a rapid increase of applications involving time series in high-dimensional systems. Central to many of these applications is the vector autoregressive (VAR) model that allows for a flexible modelling of dynamic interactions between multiple time series. In this paper we develop a simple method to test for Granger causality in high-dimensional VARs (HD-VARs) with potentially many variables.

Many financial applications consider Granger causality analysis, especially for constructing high-dimensional networks. Networks of financial firms’ intedependencies are investigated in Basu et al. (2015), Gao et al. (2017), Demirer et al. (2018) and Barigozzi and Brownlees (2019). Similarly, spillovers and contagion among stock returns are investigated in networks using Granger causality analysis in Lin and Michailidis (2017), Vỳrost et al. (2015) and Corsi et al. (2018).

Most of the econometric literature has traditionally been focused on allowing for high dimensionality in VARs through the use of factor models (see e.g. Bernanke et al., 2005; Chudik and Pesaran, 2016) or Bayesian methods (Bańbura et al., 2010). For instance Billio et al. (2012) develops measures of connectedness to assess systemic risk propagation among institutions in the financial system using principal component analysis and Granger causality networks. Recent years have seen an increase in regularized, or penalized, estimation of sparse VARs based on popular methods from statistics such as the lasso (Tibshirani, 1996) and elastic net (Zou and Hastie, 2005), which impose sparsity by setting a (data-driven) selection of the coefficients to zero.

Compared to factor models, such sparsity-seeking methods have often an advantage of interpretability, as in many economic applications, it appears natural to believe that the most important dynamic interactions among a large set of variables can be adequately captured by a relatively small – but unknown – number of ‘key’ variables. As such, the use of these methods for estimating HD-VAR models has also increased significantly in recent years, see e.g. Nicholson et al. (2017), Basu et al. (2019), Billio et al. (2019), Wilms and Croux (2018), Korobilis and Pettenuzzo (2019)).

Regularized estimation theory for high-dimensional time series and VAR models is now well established, see among others Song and Bickel (2011), Basu and Michailidis (2015), Kock and Callot (2015), Davis et al. (2016), Medeiros and Mendes (2016), Audrino and Camponovo (2018) and Masini et al. (2019) and Wong et al. (2020); Kock et al. (2020) provide a recent review. However, performing inference on HD-VARs, such as testing for Granger causality, still remains a non-trivial matter. As is well known, performing inference after model selection (post-selection inference) is complicated as the selection step invalidates ‘standard’ inference where the uncertainty regarding the selection is ignored (see Leeb and Pötscher, 2005). Complexities introduced by the temporal and cross-sectional dependencies in the VAR mean that most recently developed post-selection inference methods are not automatically applicable.

Most existing literature on Granger causality testing in HD-VARs therefore has so far not considered post-selection inferential procedures. Wilms et al. (2016) propose a bootstrap Granger causality test in HD-VARs, but do not account for post-selection issues. Similarly, Skripnikov and Michailidis (2019) investigate the problem of jointly estimating multiple network Granger causal models in VARs with sparse transition matrices using lasso-type methods, but focus mostly on estimation rather than testing. Song and Taamouti (2019) focus on statistical procedures for testing indirect/spurious causality in high-dimensional scenarios, but consider factor models rather than regularized regression techniques. Lin and Michailidis (2017) consider high-dimensional multi-block VARs derived from a two-blocks recursive linear dynamical system and use a maximum likelihood (ML) estimator for Gaussian data. In order to obtain the ML estimates for the system transition matrices and the precision matrix, respectively the lasso and graphical lasso on the residuals are iterated until convergence. Krampe et al. (2018) develops bootstrap techniques for sparse VAR models combining a model-based bootstrap procedure and the de-sparsified lasso (see Van de Geer et al. (2014)) to perform inference on the autoregressive parameters. Chaudhry et al. (2017) look at de-biased estimators as in Javanmard and Montanari (2014), for Gaussian and sub-Gaussian VAR processes with a focus on Granger-causality and control of the false discovery rate.

In this paper we build on the post-double-selection approach proposed by Belloni et al. (2014b), to develop a valid post-selection test of Granger causality in HD-VARs. The finite-sample performance depends heavily on the exact implementation of the method. In particular, the tuning parameter selection in the penalized estimation is crucial. We therefore perform an extensive simulation study to investigate the finite-sample performance of the different ways to set up the test in order to be able to give some practical recommendations. In addition, we investigate the construction of networks of realized volatilities using a sample of 30 financial stocks modeled as a vector heterogeneous VAR (Corsi, 2009). We are able to demonstrate how our approach allows for obtaining much sharper conclusions than standard low-dimensional VAR techniques.

The remainder of the paper is as follows: Section 2 introduces the high-dimensional VAR model and Granger causality tests. In Section 3 we propose our estimation and inferential framework. Section 4 establishes the asymptotic properties of our method and discusses the assumptions required for the theory to hold. Section 5 reports the results of the Monte Carlo simulations. We apply our method in Section 6 to construct volatility spillover networks. Section 7 concludes. Proofs and supplemental results can be found in the appendix.

A few words on notation. For any $n$-dimensional vector $\bm{x}$, we let ${\left\lVert\bm{x}\right\rVert}_{p}=\left(\sum_{i=1}^{n}{\left\lvert x_{i}\right\rvert}^{p}\right)^{1/p}$ denote the $\ell_{p}$-norm. For any index set $S\subseteq\{1,\ldots,n\}$, let $\bm{x}_{S}$ denote the sub-vector of $\bm{x}_{t}$ containing only those elements $x_{i}$ such that $i\in S$. ${\left\lvert S\right\rvert}$ denotes the cardinality of the set $S$. We use $\xrightarrow{p}$ and $\xrightarrow{d}$ to denote convergence in probability and distribution, respectively.

2 High-dimensional Granger causality tests

Loosely speaking, the notion of Granger causality captures predictability given a particular information set (Granger, 1969, 1980). If the addition of variable $X$ to the given information set $\Omega$ alters the conditional distribution of another variable $Y$, and both $X$ and $\Omega$ are observed prior to $Y$, then $X$ improves predictability of $Y$, and is said to Granger cause $Y$ with respect to $\Omega$. Granger (1969) originally envisioned the information set $\Omega$ “be all the information in the universe” (p. 428), which is of course not a workable concept. Yet clearly the choice of information set has a major effect on the interpretation of the finding of (non-)Granger causality, as discussed in Granger (1980). In particular, spurious Granger causality from $X$ to $Y$ may be found when both $X$ and $Y$ are Granger caused by $Z$, but $Z$ is omitted from $\Omega$. As such, one might want to include as many potentially relevant variables in the information set as possible in order to avoid finding spurious causality due to omitted variables, thereby moving as much as possible towards the universal information set envisioned by Granger. However, conditioning on so many variables leads to obvious problems of high-dimensionality rendering many standard statistical techniques invalid.

In this paper we focus on testing Granger causality in mean using linear models, in which setup the VAR model is the natural tool to investigate this problem. However, to enlarge the information set means estimating a VAR with an increasing number of variables. The number of parameters in a VAR increases quadratically with the number of time series included; an unrestricted VAR($p$) has $K^{2}p$ coefficients to be estimated, where $K$ is the number of series and $p$ is the lag-length. As the time series dimension $T$ is typically fairly small for many economic applications, the data do not contain sufficient information to estimate the parameters and consequently standard least squares and maximum likelihood methods suffer from the curse of dimensionality, resulting in estimators with high variance that overfit the data.

2.1 Granger causality testing in VAR models

Let $\bm{y}_{1},\ldots,\bm{y}_{T}$ be a $K$-dimensional multiple time series process, where $\bm{y}_{t}=(y_{1,t},\ldots,y_{K,t})^{{}^{\prime}}$ is generated by a VAR($p$) process

[TABLE]

where for notational simplicity we assume the variables have zero mean; if not they can be demeaned prior to the analysis, or equivalently a vector of intercepts is added. $\bm{A}_{1},\ldots,\bm{A}_{p}$ are $K\times K$ parameter matrices and $\bm{u}_{t}$ is a martingale difference sequence (mds) of error terms. We consider weakly stationary VAR models, as formalized in Assumption 1 below.

Assumption 1**.**

The VAR model in (1) satisfies:

1. (a)

$\{\bm{u}_{t}\}_{t=1}^{T}$ is a weakly stationary mds with respect to $\mathcal{F}_{t}=\sigma(\bm{y}_{t},\bm{y}_{t-1},\bm{y}_{t-2},\ldots)$ $\bm{u}_{t}$ such that $\mathbb{E}(\bm{u}_{t}|\mathcal{F}_{t-1})=\bm{0}$ for all $t$ and $\bm{\varSigma}_{u}=\mathbb{E}(\bm{u}_{t}\bm{u}_{t}^{\prime})$ is positive definite. 2. (b)

All roots of $\det(\bm{I}_{K}-\sum_{j=1}^{p}\bm{A}_{j}z^{j})$ lie outside the unit disc, such that the lag polynomial is invertible.

In the VAR model (1) we are interested in testing whether variables in the set $J$ Granger cause variables in the set $I$ in mean, conditional on all the other variables, where $J,I\subset\{1,\ldots,K\}$ and $J\cap I=\emptyset$. Let $N_{I}={\left\lvert I\right\rvert}$ and $N_{J}={\left\lvert N_{J}\right\rvert}$ denote the number of variables in $I$ and $J$ respectively. We describe our proecdure here in general form for testing blocks of variables. For any sets $S_{1},S_{2}\subseteq\{1,\ldots,K\}$ of variables define the best linear predictor in $L_{2}$-norm of $\bm{y}_{S_{1},t}$ given $\bm{x}_{S_{2},t-1}^{(p)}=(\bm{y}_{S_{2},t-1}^{\prime},\ldots,\bm{y}_{S_{2},t-p}^{\prime})^{\prime}$ as $\mathcal{P}(\bm{y}_{S_{1},t}|\bm{x}_{S_{2},t-1}^{(p)})=\bm{\varGamma}^{*}\bm{x}_{S_{2},t-1}^{(p)}$, where $\bm{\varGamma}^{*}=\min_{\bm{\varGamma}}\mathbb{E}\left[\left.{\left\lVert\bm{y}_{S_{1},t}-\bm{\varGamma}\bm{x}_{S_{2},t-1}\right\rVert}_{2}^{2}\right.\right]$. Then we say that $\bm{y}_{J,t}$ does not Granger cause $\bm{y}_{I,t}$ conditionally on $\bm{x}_{J^{c},t}$ if

[TABLE]

for any value of $\bm{x}_{J^{c},t}$. In other words, conditional on $\bm{x}_{J^{c},t}$, addition of the lags of $\bm{y}_{J,t}$ to the information set does not improve predictability of $\bm{y}_{I,t}$. Note that Granger (non-)causality as defined in (2) is a property of the population. In the VAR (1) this means that testing for Granger causality can be done via testing the joint significance of the blocks of coefficients in the matrices $\bm{A}_{1},\ldots,\bm{A}_{p}$ corresponding to the impact of variables $J$ on $I$.

To illustrate, consider (1) with $p=1$ lag, and assume without loss of generality that the variables in $\bm{y}_{t}$ are ordered such that $\bm{y}_{t}=\left(\bm{y}_{I,t}^{\prime},\bm{y}_{J,t}^{\prime},\bm{y}_{-(I\cup J),t}^{\prime}\right)^{\prime}$, where $-(I\cup J))$ refers to all variables not in $J$ or $I$. Then we can write

[TABLE]

where $\bm{A}$ is partitioned conformably with the blocks in $\bm{y}_{t}$. In this case, the best linear predictors in (2) are given by

[TABLE]

For any arbitrary value of $\bm{y}_{t-1}$, these can only coincide if $\bm{A}_{I,J}=\bm{0}$. Hence, the null hypothesis of no Granger causality from $J$ to $I$ in the VAR($1$) model can be formulated in terms of $\bm{A}_{I,J}=\bm{0}$. This is easily extended to $p>1$ by simply testing if the $(I,J)$-block of all $p$ lag matrices is equal to zero.

In the remainder of the paper, we will be working with a stacked representation of (1) for the variables in $I$. Specifically, let $\bm{Y}=\left(\bm{y}_{p+1},\ldots,\bm{y}_{T}\right)^{\prime}$ and let $\bm{y}_{I}=\operatorname{vec}\left(\bm{Y}_{I}\right)$ denote the $N_{I}\times 1$ stacked vector containing all observations corresponding to the variables in $I$. Similarly, let $\bm{u}_{I}=\operatorname{vec}(\bm{U}_{I})$, where $\bm{U}=\left(\bm{u}_{p+1},\ldots,\bm{u}_{T}\right)^{\prime}$. Let $\bm{X}=\left(\bm{x}_{p}^{(p)},\ldots,\bm{x}_{T-1}^{(p)}\right)^{\prime}$ and $\bm{X}^{\otimes}=\bm{I}_{N_{I}}\otimes\bm{X}$, while defining the stacked parameter vector $\bm{\beta}=\operatorname{vec}((\bm{A}_{1},\ldots,\bm{A}_{p})^{\prime})$. Then we can write

[TABLE]

where $\bm{X}_{GC}^{\otimes}=\bm{I}_{N_{I}}\otimes\bm{X}_{GC}$, and $\bm{X}_{GC}=\left(\bm{x}_{J,p}^{(p)},\ldots,\bm{x}_{J,T-1}^{(p)}\right)^{\prime}$ contains those columns of $\bm{X}$ corresponding to the potentially Granger causing variables in $J$; $\bm{X}_{-GC}$ and $\bm{X}_{-GC}^{\otimes}$ are then defined similarly but containing the remaining variables.111Note that if $I=\{i\}$ for one particular value of interest, then (4) simply corresponds to a single equation from the VAR in (1). Testing for no Granger causality is then equivalent to testing $H_{0}:\bm{\beta}_{GC}=\bm{0}$ against $H_{1}:\bm{\beta}_{GC}\neq\bm{0}$.

Define $N_{J}={\left\lvert J\right\rvert}$ and $N_{I}={\left\lvert I\right\rvert}$. Note that $\bm{\beta}_{-GC}$ has $\left(K-N_{J}\right)\times N_{I}\times p$ elements, which we assume large through having a large number of variables $K$. On the other hand, throughout the paper we assume that $N_{J}$, $N_{I}$ and $p$ are small, or more precisely, fixed when sample size increases to infinity. As $\bm{\beta}_{GC}$ has $N_{GC}=N_{J}\times N_{I}\times p$ elements, these are also implied to be fixed. While theoretically it is possible to consider an increasing number of elements in $\bm{\beta}_{GC}$ (see Remark 6 for details), it would not be required for typical applications. $J$ and $I$ are under the researcher’s control and in most applications it is natural to consider a small number of variables of interest; often both $J$ and $I$ will only consist of a single variable, as in our application.

For $p$ it may appear more restrictive to assume it small. However, large $p$ in univariate regressions or small systems often arise from neglected dynamics with omitted variables (Hecq et al., 2016). As our HD-VAR attempts to include many more variables than typical small systems, we hope to alleviate the omitted variable issue, and thereby also directly making smaller $p$ much more realistic. Of course, $p$ is generally unknown in practice. However, in many applications it is possible to give a reasonable (and small) upper bound on $p$, which is sufficient for our algorithm. If not, $p$ has to be estimated. We discuss two ways in the next section.

Remark 1**.**

Our operational version of Granger causality only considers causality in mean. Additionally, one might argue that considering only linear models is a further restriction on the generality of the concept of Granger causality. However, in our high-dimensional approach linear models are less restrictive as would appear. First, the VAR does not have to be formulated for levels of variables of interest. In fact, in our application we formulate a VAR for (realized) variances, such that we are implicitly testing Granger causality in second moments rather than first moments. Second. the linear VAR model in many cases provides a good approximation to a general nonlinear process via a Wold-type representation argument, see e.g. Meyer and Kreiss (2015). Finally, non-linear transformations (such as powers) of the original variables can be added to (4), by which general functional forms can be approximated (even if one then strictly loses the VAR equivalence). While in small systems this is infeasible as it increases the dimensionality dis-proportionally, our high-dimensional approach can handle this without any conceptual issues. In fact, Belloni et al. (2014a) explicitly motivate their high-dimensional linear approach as an approximation to a general function; their arguments apply here as well.

3 Inference after selection by the lasso

In this section we introduce our inferential procedure to the Granger causality tests in high-dimensional VARs. We first discuss the lasso, which we use in the initial stage to select relevant variables. Next we discuss how naive use of the lasso introduces post-selection problems for inference, and we propose our algorithm to remedy this.

3.1 The lasso estimator

As $\bm{\beta}$ is high-dimensional when $Kp$ is large relative to $T$, least squares estimation is not appropriate, and a structure must be imposed on $\bm{\beta}$ to be able to estimate it consistently. We assume sparsity of $\bm{\beta}$; that is, we assume that $\bm{\beta}$ can accurately be approximated by a coefficient vector with a (significant) portion of the coefficients equal to zero.

The sparsity assumption validates the use of variable selection methods, thereby reducing the dimensionality of the system without having to sacrifice predictability. For a general $n$-dimensional vector of responses $\bm{y}$ and $n\times M$-dimensional matrix of covariates $\bm{X}$, the (weighted) lasso simultaneously performs variable selection and estimation of the parameters by solving

[TABLE]

where $\lambda$ is a non-negative tuning parameter determining the strength of the penalty, and $\{w_{m}\}_{m=1}^{M}$ are non-negative weights corresponding to the parameters in $\bm{\beta}$. For the standard lasso the weights are either equal to one, or equal to zero (if this parameter should not be penalized). The notation $\hat{\bm{\beta}}(\lambda)$ highlights that the solution to the minimization problem depends on $\lambda$, which has te be selected as well (see Section 3.3). When no confusion can arise, we simply write $\hat{\bm{\beta}}$.

One may also consider the adaptive lasso (Zou, 2006) with parameter-specific weights $w_{j}$ in (5) based on an initial estimation of $\bm{\beta}$, which is able to delete more irrelevant variables. However, for our purpose such oracle properties are not very relevant; we wish to eliminate the effects of the other “nuisance” variables on the relation between the variables tested for Granger causality, but we do not need to identify which of these nuisance variables matter.

Theoretical properties of lasso estimation in stable VAR models have now been studied extensively. We here non-exhaustively mention some of the key results for our setting; see Kock et al. (2020) for a thorough review. Kock and Callot (2015) derive oracle properties of the adaptive lasso for VAR models. Basu and Michailidis (2015) establish restricted eigenvalue conditions for VAR models and show their sufficiency for estimation consistency. Medeiros and Mendes (2016) relax the Gaussianity assumptions of these papers by considering conditionally heteroskedastic errors, and demonstrate that the adaptive lasso retains oracle properties in time series settings. Finally, Masini et al. (2019) derive bounds on estimation errors in approximately sparse VAR models under very general conditions, allowing for heavy tails and dependence in the error terms. In particular, they show that several commonly used volatility processes in financial research satisfy these assmuptions, thereby formally establishing the suitability of the lasso for many financial applications of VAR models.

3.2 Post-Double-Selection Granger causality test

3.2.1 The need for post selection inference

One might be tempted to simply perform the (adaptive) lasso as in (5) on (4), setting $w_{GC}=0$, and then testing whether $\bm{\beta}_{GC}=0$, potentially after re-estimating the model by OLS on only the selected variables. However, this ignores the fact that the final, selected, model is random and a function of the data. The randomness contained in the selection step means the post-selection estimators do not converge uniformly to a normal distribution, as the potential omitted variable bias from omitting (weakly) relevant variables in the selection step is too large to maintain uniformly valid inference.

In a sequence of papers (see e.g. Leeb and Pötscher, 2005), Leeb and Pötscher address these issues, showing that distributions of post-selection estimators only converge point-wise but not uniformly in the parameter space to normal distributions. Therefore, “standard” asymptotics fail to deliver a proper approximation of finite-sample behavior due to the presence of small, hard to detect parameters, whose omitted variable bias is too large to ignore asymptotically. As such, post-selection based on oracle properties is only appropriate if one a priori rules out small parameters conditions (via beta-min conditions, see e.g. van de Geer et al., 2011) thus obtaining a sharp separation of non-zero from zero coefficients. This is typically far too strong to be reasonable in applications, and methods explicitly accounting for selection are required.

Several approaches to valid post-selection inference, also referred to as honest inference, have been developed in recent years based on various philosophies, such as simultaneous inference across models (Berk et al., 2013), inference conditional on selected models (Lee et al., 2016), or debiasing (desparsifying) the lasso estimates (Van de Geer et al., 2014; Zhang and Zhang, 2014). We focus on the double selection approach developed by Belloni, Chernozhukov and co-authors; see e.g. (Belloni et al., 2014a) for an overview. This approach is tailored for the lasso, easy to implement, and can be extended to dependent data.

Belloni et al. (2014) develop a post-double-selection approach to construct uniform inference for treatment effects in partially linear models with high-dimensional controls using the lasso. Two initial lasso estimations of both the outcome and the treatment variable on all the controls are performed, and a final post-selection least squares estimation is conducted of the outcome variable on the treatment variable and all the controls selected in at least one of the two steps. The double variable selection step substantially diminishes the omitted variable bias and ensures the errors of the final model are (close enough to) orthogonal with respect to the treatment. The authors proved uniform validity of the procedure under a wide range of DGPs, including heteroskedastic and non-Gaussian errors.

Chernozhukov et al. (2020) extend the analysis of estimation and inference for highly-dimensional systems in regressions, allowing for (weak) temporal and cross-sectional dependency. Regularization techniques for dimensionality reduction are applied iteratively in the system and the overall penalty is jointly chosen by a block multiplier bootstrap procedure. Oracle properties and bootstrap consistency of the test procedure are derived. Furthermore, simultaneous valid inference is obtained via algorithms employing least square or least absolute deviation after (double) lasso selection step(s). Although our approach is closely related to that of Chernozhukov et al. (2020), it differs in a number of ways. Our method is simpler and faster to implement as it does not rely on bootstrap methods. Also, Chernozhukov et al. (2020) focus on general systems of equations and general ways of performing inference, which is different from our specific focus on Granger causality and VAR models. Third, we consider a different set of assumptions to establish the validity of our method, where we specifically focus on the relevance of these assumptions for applications in financial econometrics.

3.2.2 High-dimensional Granger causality test

We here describe how to implement the post-double-selection procedure in a VAR context. Let $\bm{x}_{GC,j}$, $j=1,\ldots,N_{X}$, where $N_{X}=pN_{J}$, denote the $j$-th column of $\bm{X}_{GC}$ and consider the partial regressions:

[TABLE]

where $\bm{\gamma}_{j},\;j=0,\ldots,N_{X}$, are the best linear prediction coefficients222Note that Assumption 2(a) implies that $\left(\mathbb{E}\bm{x}_{-GC,t-1}\bm{x}_{-GC,t-1}\right)^{-1}$ and hence $\left(\mathbb{E}\bm{X}_{-GC,t-1}^{\otimes}\bm{X}_{-GC,t-1}^{\otimes\prime}\right)^{-1}$ exist.

[TABLE]

where $\bm{X}_{-GC,t}^{\otimes}=\bm{I}_{N_{I}}\otimes\bm{x}_{-GC,t-1}$. As the errors $\bm{e}_{0},\ldots,\bm{e}_{N_{X}}$ are orthogonal to $\bm{X}_{-GC}$, partialling out the effects of these variables would allow for a valid test of Granger causality. Of course, (6) and (7) are still high-dimensional and cannot be estimated by least squares. However, we can select the relevant variables from lasso estimation of (6) and (7) and collect all these for the final estimation of $\bm{y}_{I}$ on $\bm{X}_{GC}^{\otimes}$ plus only those relevant variables.

Intuitively, this works because to cause omitted variable bias on the coefficients of $\bm{X}_{GC}$, a particular variable in $\bm{X}_{-GC}$ must have a nonzero coefficient in both (6) and one of the regressions in (7). If its coefficient is zero in (6), it has no effect on $\bm{y}_{I}$ and is therefore not wrongfully omitted. If it has a zero coefficient in all regressions in (7), it is not correlated with any variables of interest, and omitting it will not result in a bias. By including all variables that are selected in at least a single of these regressions, we essentially allow for “one free mistake” by the lasso in failing to select a relevant variable. That is, omitted variable bias will only occur if the lasso fails to select a relevant variable in both regressions simultaneously. As the probability of this occurring decreases quadratically, this is sufficient to be negligible asymptotically and allow for uniformly valid inference. We provide a formal justification in Section 4.

We now state the details of our algorithm which executes the post-double-section along the lines described above, and conclude this section with some remarks.

Remark 2**.**

We perform the initial regressions in terms of $\bm{X}_{GC}$ amd $\bm{X}_{-GC}$ instead of $\bm{X}_{GC}^{\otimes}$ and $\bm{X}_{-GC}^{\otimes}$. The two are equivalent, as the Kronecker product essentially just copies the columns of $\bm{X}$ both in the dependent and explanatory variables. Running the initial regressions in terms of $\bm{X}^{\otimes}$ therefore essentially means running the same regression $N_{I}$ times, which is pointless as the selected variables remain the same in terms of the columns of $\bm{X}_{-GC}$. We therefore perform the regressions just once for each column in $\bm{X}_{GC}$. The construction of $\hat{S}_{X}^{\otimes}$ ensures that for any selected column $\bm{x}_{-GC,m}$, we select every column of $\bm{X}_{-GC}^{\otimes}$ in which $\bm{x}_{-GC,m}$ appears.

Remark 3**.**

The feasible generalized least squares (FGLS) estimation in Step [2] is needed when $N_{I}>1$ to account for the correlation between equations of the VAR, and the fact we do not have the same selected regressors in each equation, as those coming from (6) differ. Note that if $N_{I}=1$, FGLS estimation collapses to the familiar form of the LM statistic. In that case one regresses $\hat{\bm{\xi}}$ by OLS onto the variables retained by the previous regularization steps plus the Granger causality variables, and retain the residuals $\hat{\bm{\nu}}=\hat{\bm{\xi}}-\bm{X}_{\hat{S}\;\cup\;GC}^{\otimes}\hat{\bm{\beta}}^{*}$, obtaining $R^{2}=1-\hat{\bm{\nu}}^{\prime}\hat{\bm{\nu}}/\hat{\bm{\xi}}^{\prime}\hat{\bm{\xi}}$.

Remark 4**.**

Our lasso estimation of an HD-VAR can be interpreted as a general, data-driven, approach to Granger causality testing which encompasses the theory-driven ‘standard’ approach in low-dimensional VARs. In particular, the lasso can be interpreted as imposing (approximate) sparsity over a high-dimensional information set, with the extent and location of the sparsity, or irrelevance, determined in a data-driven way. Conversely, testing Granger causality in a low-dimensional setting can then be interpreted as a priori assuming an extreme degree of sparsity over the same information set; in other words, it amounts to assuming that none of the additional series are relevant.

Remark 5**.**

Given that we essentially have $N_{GC}=N_{J}\times N_{I}\times p$ steps of selection, it would be more appropriate to refer to our method as “post-$N_{GC}$-selection” approach. For expositional simplicity however we stick to the post-double-selection name, as this is the common name for such a procedure, and conveys the essence of our method equally well.

Remark 6**.**

Although the lasso regressions can handle increasing $N_{J}$, $N_{I}$ or $p$ with any issues, inference becomes more complicated when $N_{GC}$ increases with the sample size as the proposed LM statistic (or similarly a Wald test) will not have a limit distribution anymore. In such a case one could use recently developed Gaussian approximations of maxima of high-dimensional vectors (Chernozhukov et al., 2013; Zhang and Wu, 2017) to base a test statistic on $\max_{m=1,\ldots,N_{GC}}{\left\lvert\hat{\beta}_{GC,m}^{\textsc{pds}}\right\rvert}$, where $\hat{\bm{\beta}}_{GC}^{\textsc{pds}}$ are the coefficients of $\bm{X}_{GC}^{\otimes}$ in a regression of $\bm{y}_{I}$ on $\bm{X}_{GC}$ and $\bm{X}_{\hat{S}^{\otimes}}^{\otimes}$ as in (8). However, the critical values of this test statistic have to be simulated, which complicates the testing. As we argued in Section 2 that a fixed $N_{GC}$ is a reasonable assumption for typical Granger causality applications, we do not pursue this route.

Remark 7**.**

In Step [1] we propose not to consider the GC variables in the first regularization and insert them back at Step [2]. Alternatively, the GC variable(s) can be left in the regression, such that, we regress on the full $\bm{X}^{\otimes}$ matrix. In this case there are then two further possibilities by either penalizing these variables or not. Simulations for these two alternatives have been carried out and in practice we do not find significant differences among the three in terms of size and power. The approach proposed in Step 1 delivers the best results in terms of size.

Remark 8**.**

When the time series length is of same magnitude as the number of covariates, information criteria and time series cross-validation tend to break down and select too many covariates in order to perform a post-selection by OLS. To overcome this issue we propose to place a lower bound on the penalty to ensure that in each selection regression at most $c\,TN_{I}$ variables are selected, for some $0<c<1$. In our simulation and empirical studies we set $c=0.5$. Note that, as we have $N_{GC}$ selection steps, the possibility remains that different variables are selected in each steps, making the number of variables in the union $\hat{s}$ still too large to perform the post-selection OLS, although this problem is likely to occur far less often. This can be addressed by ensuring that fewer than $N_{I}T/N_{GC}=T/N_{X}$ variables are selected in each selection step. We do not impose this stricter bound in general, as it will often be much too strict. Instead, we recommend to only address this issue if it arises in practice by an ad-hoc increase of the lower bound on the penalty.333Although it happens less often, the theoretical plug-in method for the tuning parameter occasionally also selects too many variables to make the post-OLS estimation infeasible. However, for this method no easy solution is available for bounding the penalty. One could increase the constant in the plug-in expression, thus strengthening the penalty, but this would be a rather ad-hoc adjustment. In particular, imposing the lower bound for the other methods only limits the allowed range of the tuning parameter, forcing the minimization to choose another (local) minimum that can still be far away from the boundary and justified graphically. For the plug-in method it is however difficult to justify the right amount of the increase, as the tuning parameter will be fixed to that value, and thus the chosen increase is rather arbitrary.

Remark 9**.**

Although our Granger causality test has a $\chi^{2}$ distribution under the null hypothesis asymptotically, in smaller samples the test might still suffer from the usual small-sample approximation error. As such we propose a finite-sample correction to the test in Step [3b], which in our simulation studies improved the size of our test.

Remark 10**.**

Instead of the (adaptive) lasso, other estimators can be used in Step [1] as long as they deliver a sparse coefficient vector. For instance, the elastic net of Zou and Hastie (2005) that adds an $\ell_{2}$-penalty in addition to the $\ell_{1}$-penalty of the lasso can be used. The additional penalty ensures that the elastic net is strictly convex, and as a consequence tends to select highly correlated variables as a group together, whereas the lasso would tend to select only one of these variables (Zou and Hastie, 2005). Given the typically strong correlations between many economic variables, this appears particularly useful for our context. However, we used the elastic net for both the simulations and the empirical application, and in both cases we found that the results are widely comparable to those of lasso. Therefore we chose to omit them from the paper; they are available upon request.

Remark 11**.**

One can also perform a standard Wald test of Granger causality instead of the LM test, by regressing the variables of interest on $\bm{X}_{GC}^{\otimes}$ and $\bm{X}_{\hat{S}^{\otimes}}^{\otimes}$, and testing for the significance of the coefficients of $\bm{X}_{GC}^{\otimes}$. While asymptotically the LM and Wald tests behave equally, differences might arise in small samples. We investigated the Wald version of the test in simulations as well, with results reported in Appendix B, Table B.2. In general, differences between the two methods are negligible. However, for the Wald test, occasionally we run into the problem described in Remark 8, where even with the imposed lower bound on the penalty, too many variables are selected for performing a post-selection OLS. For this reason we prefer the LM version.

3.3 Tuning parameter selection

Appropriate selection of the lasso tuning parameter $\lambda$ in (5) is crucial to achieve good performance. Many different data-driven methods exist giving wildly varying results. We provide a systematic comparison of several popular methods discussed in the literature in a simulation study. To the best of our knowledge, this is the first such comparison in the context of post-selection inference. We now introduce the methods considered in our study.

One option is to minimize an information criterion (IC) to determine an appropriate data-driven $\lambda$. Let $\hat{S}(\lambda)=\left\{m\in\{1,\ldots,Kp\}:{\left\lvert\hat{\beta}_{m}(\lambda)\right\rvert}>0\right\}$ denote the set of active variables in the lasso solution for a given $\lambda$. For a generic response vector $\bm{y}$ and predictor matrix $\bm{X}$, the value $\lambda^{IC}$ is found as

[TABLE]

where $C_{T}$ is the penalty specific to each criterion. We consider the Akaike information criterion (AIC) by Akaike (1974) with $C_{T}=2$, the Bayesian information criterion (BIC) by Schwarz (1978) with $C_{T}=\ln(T)$, and the Extended Bayesian information criterion (EBIC) by Chen and Chen (2008) with $C_{T}=\ln(T)+2\gamma\ln(Kp)$ with $\gamma=0.5$ proposed by Chen and Chen (2012) who argue that BIC fails to select the correct variables when the number of parameters is larger than the sample size.

An alternative approach is to plug in estimates of theoretically optimal values (see e.g. Bickel et al., 2009; Belloni and Chernozhukov, 2013; Belloni et al., 2011). The lasso requires that $\lambda\geq c{\left\lVert\bm{X}^{\prime}\bm{u}\right\rVert}_{\infty}/T$ for some constant $c>0$ with “high probability”. The central limit theorem motivates a Gaussian approximation where one chooses \lambda^{th}=\frac{2c\hat{\sigma}}{\sqrt{T}}\Phi^{-1}\bigg{(}1-\frac{\alpha}{2N}\bigg{)} for a small $\alpha=o(1)$, where $\Phi^{-1}(\cdot)$ is the inverse of standard Gaussian cumulative distribution function and $\hat{\sigma}$ is an estimate the variance of $\bm{u}$. In this paper we set $\alpha=0.05/\ln(T)$ and $c=0.5$, while we follow Belloni et al. (2012) in the estimation of $\sigma$. Specifically, we obtain an initial (conservative) estimate by least squares estimation of $\bm{y}$ on the five most correlated regressors. This estimate is then updated iteratively, for details see Belloni et al. (2012).

Perhaps the most popular way to choose the tuning parameter is cross-validation (CV), although CV is not always appropriate in the time series setup without modifications (Bergmeir et al., 2018). To estimate the tuning parameter with CV in a time series setup (TSCV) we use an expanding window out-of-sample forecasting scheme and minimize its squared forecasting error. The rolling window is set up with $80\%$ of the sample for training and $20\%$ for testing. Cross-validation is appealing since it does not require any plug-in estimates, however, as observed in Chetverikov et al. (2020) it typically yields small values of $\lambda$ thus still gaining fast convergence rate but at the price of less variable selection.

Remark 12**.**

Although we assume $p$ fixed, in practice it may still need to be estimated if no reasonable value (or upper bound) can be given. As $p$ determines the number of selection regressions to be conducted, it has to be determined a priori and cannot be integrated in the lasso estimation. It can still be determined though by a (separate) lasso-type algorithm. For example, one may estimate (4) with a large initial lag length $p^{*}$, and let $p$ be determined as the largest lag for which variables are selected, possibly also varying the lag length over variables. For this approach the hierarchical penalties of Nicholson et al. (2020) provide a better option than the regular lasso, as the regular lasso tends to select occasional “spurious” high lags, which would have a significant impact on the testing procedure. Alternatively one may marginalize the VAR to a collection of univariate AR($p$) processes, and select the lag length by minimizing an information criterion on the residual covariance matrix. As marginalization increases the lag length, such an approach would yield a simple to compute upper bound on $p$.

4 Asymptotic Properties

In this section we derive the asymptotic properties of our method. We first present and discuss our general high-level assumption under which the properties are derived, and then state our main results.

Assumption 2**.**

Let $\delta_{T}$ and $\Delta_{T}$ denote sequences such $\delta_{T},\Delta_{T}\rightarrow 0$ as $T\rightarrow\infty$. Then assume that the following conditions are satisfied:

1. (a)

Population Eigenvalues: Let $\bm{e}_{t}=\left(e_{1,t},\ldots,e_{N_{X},t}\right)^{\prime}$, $\bm{E}=\left(\bm{e}_{1},\ldots,\bm{e}_{T}\right)^{\prime}$ and $\bm{E}^{\otimes}=\bm{I}_{N_{I}}\otimes\bm{E}$. Define

[TABLE]

Then there exists a constant $c_{L}>0$ not depending on $T$ and $k$ such that $\lambda_{\min}(\bm{\varSigma})>c_{L}$, where $\lambda_{\min}(\bm{\varSigma})$ denotes the minimum eigenvalue of $\bm{\varSigma}$. 2. (b)

Limit Behavior: Let

[TABLE]

where $\bm{\varOmega}=\operatorname*{plim}_{T\rightarrow\infty}\left(\bm{E}^{\otimes\prime}\bm{u}_{I}\bm{u}_{I}^{\prime}\bm{E}^{\otimes}\right)/T$. 3. (c)

Empirical Process: We have with probability at least $1-\Delta_{T}$ that ${\left\lVert\bm{X}_{-GC}^{\prime}\bm{u}_{i}/\sqrt{T}\right\rVert}_{\infty}\leq\gamma_{T}$ for all $i\in I$ and ${\left\lVert\bm{X}_{-GC}^{\prime}\bm{e}_{j}/\sqrt{T}\right\rVert}_{\infty}\leq\gamma_{T}$ for all $j=1,\ldots,N_{X}$, with $\bm{e}_{j}$ the $j$-th column of $\bm{E}$, for some deterministic sequence $\gamma_{T}$ subject to the restrictions in (h). 4. (d)

Boundedness: The (Granger causality) parameters of interest are bounded, that is, there exists a fixed constant $C>0$ such that ${\left\lVert\bm{\beta}_{GC}\right\rVert}_{1}\leq C$. 5. (e)

Consistency: The initial estimators $\hat{\bm{\gamma}}_{j}$ are consistent in the prediction sense; specifically, with probability at least $1-\Delta_{T}$ we have that

[TABLE] 6. (f)

Sparsity: Let $S_{j}=\{m:\gamma_{m,j}\neq 0\}$ denote the sets of active variables in (6) and (7) and let $s={\left\lvert S_{0}\right\rvert}+\sum_{j=1}^{N_{X}}{\left\lvert S_{j}\right\rvert}$. Let $\hat{s}$ be as defined in Algorithm 1. Then both the DGP and the initial estimators are sufficiently sparse; in particular, we have that with probability at least $1-\Delta_{T}$, $\max(s,\hat{s})\leq\bar{s}_{T}$ for a deterministic sequence $\bar{s}_{T}$ subject to the restrictions in (h). 7. (g)

Sparse Eigenvalues: for any $\bm{\gamma}\in\mathbb{R}^{(K-N_{J})p}$ with ${\left\lVert\bm{\gamma}\right\rVert}_{0}\leq\bar{s}_{T}$, we have with probability at least $1-\Delta_{T}$ that ${\left\lVert\bm{\gamma}\right\rVert}_{2}^{2}\leq{\left\lVert\bm{X}_{-GC}\bm{\gamma}/\sqrt{T}\right\rVert}_{2}^{2}/\phi_{T,\min}^{2}$, where $\phi_{T,\min}>0$ is subject to the restrictions in (h). 8. (h)

Rate Conditions: The deterministic sequences $\bar{s}_{T},\gamma_{T}$ and $\phi_{T,\min}$ introduced above satisfy the restriction $\bar{s}_{T}\gamma_{T}/\phi_{T,\min}\leq\delta_{T}T^{1/4}$.

Assumption 2 is a high-level assumption that allows for much flexibility on the underlying DGP and the used estimators in the first step. We now discuss each part in turn. Part (a) assumes that the minimum eigenvalue of $\bm{\varSigma}$ is bounded. This is required for application of lasso methods, as well as for the inverse covariance matrix $\bm{\varSigma}^{-1}$ and the projection coefficients in (6) and (7) to exist. Part (b) assumes that a central limit theorem and weak law of large numbers hold. Essentially this require that the process is sufficiently well-behaved in terms of moments and dependence allowed. Although for convenience we assume martingale difference errors in Assumption 1, (b) holds under much weaker conditions such as mixing errors; see e.g. Davidson (1994, Chapter 14).

Part (c) is closely related to (b), but additionally controls the tail behavior of the empirical process. Results of this kind are standard in the lasso literature and can be derived using a variety of tail bounds depending on the properties of the random variables of interest, see e.g. Kock and Callot (2015) and Medeiros and Mendes (2016) for results relevant to VAR and time series models. Of particular interest for financial applications, Masini et al. (2019, Proposition 2) show that this condition is satisfied for VAR models with general weakly dependent erorrs that include many popular multivariate volatility models. The boundedness assumption in (d) is not very restrictive, and with $N_{GC}$ fixed follows directly if the parameter space of $\bm{\beta}$ is a compact set.

Part (e) imposes an appropriate consistency rate on the predictions coming from the first-stage estimator. Such prediction consistency is a standard result for lasso estimators; in particular, Wong et al. (2020) obtain it for a very general class of VAR models allowing for conditional heteroskedasticity and dependence in the error terms. Adamek et al. (2020) derive consistency of the lasso under misspecified time series models, and show that their setting covers (among others) the first-step regressions of the relevant predictors in $\bm{X}_{GC}$ on the other regressions, which are inherently misspecified in a VAR setup due to the missing lags; see their Remark 3 for further details.

Next to consistency, we also require sparsity of the DGP and the estimator, as controlled by part (f). The assumption of exact sparsity in the DGP for the initial regressions can be relaxed to approximate sparsity as in Belloni et al. (2014b). For the sake of expositional clarity we do not work under that assumption here but stick to the simpler exact sparsity. Sparsity of the first-stage estimator is needed in our framework as we perform OLS on the selected variables from the first-stage regressions. If the selected variables are not sparse enough, too many variables will be selected for OLS to be feasible. Sparsity of lasso estimators is analysed in Belloni and Chernozhukov (2013), while Kock and Callot (2015) and Medeiros and Mendes (2016) provide results for adaptive lasso for time series. Importantly, we do not require consistent model selection; the selection method used is allowed to make “persistent” mistakes, allowing for both variables to be incorrectly included and relevant variables to be missed, as long as the estimator remains sufficiently sparse and consistency is guaranteed. Unlike Belloni et al. (2014b), we allow for the order of sparsity of the estimator to differ from the true sparsity thereby opening the way for conservative selection procedures.

Given the assumptions above, the eigenvalue assumption in (g) becomes almost superfluous, as it is generally needed to establish (e) and (f) for lasso-type estimators; see e.g. Belloni and Chernozhukov (2013) and Medeiros and Mendes (2016) for details. It requires that for sufficiently sparse vectors, the eigenvalues of the subset of the Gram matrix corresponding to their non-zero support do not decrease to zero too fast. Such assumptions are standard in the lasso literature in various guises as restricted eigenvalue conditions, and can typically be derived by making similar conditions on the population covariance matrix $\bm{\varSigma}_{-GC,-GC}$ coupled with a convergence result of the Gram matrix $\bm{X}_{-GC}^{\prime}\bm{X}_{GC}$ to $\bm{\varSigma}_{-GC,-GC}$. Basu et al. (2015), Masini et al. (2019) and Wong et al. (2020) establish the plausibility of such restricted eigenvalue conditions for various VAR models. We state the condition here explicitly as it is needed directly in the proofs.

Finally, note that the restrictions on tail behavior (via $\gamma_{T}$), sparsity (via $\bar{s}_{T}$) and minimum eigenvalues (via $\phi_{T,\min}$) are meaningless if no rates on these sequences are imposed. Part (h) therefore is the key part which connects all assumptions with explicit rates needed for the validity of the PDS method. The restrictions here represents a trade-off between sparsity, thickness of tails and minimum eigenvalues. For example, if, as often assumed $\phi_{T,\min}$ is fixed and $\bm{u}_{t}$ is Gaussian, tails are sufficiently thin that $\gamma_{T}$ can be chosen as roughly the order of $\sqrt{\ln(K^{2}p)}$ (cf. Kock and Callot, 2015, Lemma 4), leaving room for either almost exponentially large $K$ relative to $T$, or a fairly non-sparse model. On the other hand, if only $m$ moments of $\bm{u}_{t}$ exist, $\gamma_{t}$ should be taken roughly of the order $(K^{2}p)^{2/m}$ (Masini et al., 2019, Lemma 2), requiring polynomial growth of $K$ compared to $T$ and sparser models.

The most restrictive and crucial assumption needed on the underlying DGP for satisfying Assumption 2 is the sparsity of the underlying DGP formulated in part (f). The plausibility of this assumption highly depends on the specific application. In many financial applications sparsity (or its approximate version) is natural, for example in portfolio selection when the number of assets is large and the estimation of high-dimensional volatility matrices in financial risk assessment (see Fan et al. (2011) for an overview), as well as in our investigation of Granger causality in networks of realized volatilities in Section 6. The volatility of one particular stock is likely to have specific channels of contagion rather than affecting the whole stock market at the same time. Shocks to one asset therefore likely propagate through the system via specific channels, which corresponds to sparse lag polynomials. One might worry about systemic shocks affecting many assets; however, the dense covariance matrix $\bm{\varSigma}_{u}$ can accommodate simultaneous common shocks. Moreover, the dynamic of such shocks can generally well be captured through a sparse combination of the most important and most affected assets. Similarly, in macroeconomic applications it has been found that a few important variables can capture the effects of unobserved common factors, leading sparse models to perform as well as common factors (De Mol et al., 2008; Smeekes and Wijler, 2018).

We are now ready to state our main asymptotic result of this section in Theorem 1 which establishes the asymptotic normality of the post-lasso (generalized) least squares estimator. Here we slightly deviate from the LM test in Algorithm 1; after the double selection procedure carried out in Step [1], we regress the transformed outcome variables $\widetilde{\bm{y}}_{t}=(\bm{G}_{T}\otimes\bm{I}_{T})\bm{y}_{I}$ on both the Granger causing $\widetilde{\bm{X}}_{GC}^{\otimes}=(\bm{G}_{T}\otimes\bm{I}_{T})\bm{X}_{GC}^{\otimes}$ and selected variables $\widetilde{\bm{X}}_{\hat{S}^{\otimes}}^{\otimes}=(\bm{G}_{T}\otimes\bm{I}_{T})\bm{X}_{\hat{S}^{\otimes}}^{\otimes}$

[TABLE]

The transformation by the matrix $\bm{G}_{T}$ allows for the GLS esitmation needed in the LM procedure by taking $\bm{G}_{T}=\hat{\bm{\varSigma}}_{u,I}^{-1/2}$, while OLS is performed with $\bm{G}_{T}=\bm{I}_{N_{I}}$. In the latter case the theorem provides the foundation for the Wald test discussed in Remark 11 (minus the required variance estimation for that test). We state this result separately as it is interesting in its own right, and can be used to establish validity of other tests such as the Wald test.

Theorem 1**.**

Let $\hat{\bm{\beta}}_{GC}^{\textsc{pds}}$ denote the OLS estimator of $\bm{\beta}_{GC}^{\textsc{pds}}$ in 8. Let $\bm{G}_{T}$ be any matrix satisfying with probability at least $1-\Delta_{T}$ that $0<c_{1}\leq\lambda_{\min}(\bm{G}_{T}^{\prime}\bm{G}_{T})\leq{\left\lVert\bm{G}_{T}^{\prime}\bm{G}_{T}\right\rVert}_{\max}\leq c_{2}<\infty$, where $c_{1},c_{2}$ are constants not depending on $T$. Then, uniformly in all DGPs that satisfy Assumption 2, we have as $T\rightarrow\infty$,

[TABLE]

where $\bm{\varOmega}_{\bm{G}}=\operatorname*{plim}_{T\rightarrow\infty}\left[(\bm{G}_{T}^{\prime}\bm{G}_{T}\otimes\bm{E}^{\prime})\bm{u}_{I}\bm{u}_{I}^{\prime}(\bm{G}_{T}^{\prime}\bm{G}_{T}\otimes\bm{E})\right]/T$.

Theorem 1 establishes the asymptotic normality of the post-double-selection OLS estimators. The statement ‘uniformly in all DGPs that satisfy Assumption 2’ should be interpreted as the theorem holding uniformly over a parameter space that is defined such that Assumption 2 holds for all parameters in that parameter space. Importantly, no beta-min conditions on the smallest magnitude of parameters are required, thus alleviating the post-selection inference problem. We refer to Comments 3.4 and 3.5 in Belloni et al. (2014b) for further details regarding the uniformity. The limit distribution of the LM test now follows straightforwardly from Theorem 1, and is stated in the corollary below.

Theorem 2**.**

Let $\bm{\beta}_{GC}=\bm{0}$. Then, uniformly in all DGPs that satisfy Assumption 2 and for which $\bm{\varOmega}=\bm{\varSigma}_{u,I}\otimes\bm{\varSigma}_{GC|-GC}$, we have that

[TABLE]

Theorem 2 establishes the limiting distribution of the PDS-LM test under an additional condition on the (co)variances of the partial regression errors, which is satisfied if the errors are iid. To allow for heteroskedaticity the LM test has to be modified, which would only lead to more cumbersome proofs without adding any novelty specific to the high-dimensional case. Therefore we focus on the homoskedastic case here, although we do consider a heteroskedasticity-robust version of the test in the volatility application in Section 6.444Note that this is no different for the Wald test, for which the variance estimation has to be adjusted as well.

5 Monte-Carlo Simulations

We now evaluate the finite-sample performance of our proposed Granger causality test. We consider three Data Generating Processes (DGPs) inspired by Kock and Callot (2015):

[TABLE]

The diagonal VAR in DGP1 respects the sparsity assumption while in DGP2 the entries are set to decrease exponentially fast in the distance from the main diagonal and hence the sparsity assumption is not met. DGP2 could be empirically motivated by looking e.g. at financial interconnectedness. Financial institution, such as banks, lend to and borrow from one another becoming interconnected through interbank credit exposures. The financial distress experienced by one bank is likely to be most heavily transmitted the closer the connections are as well as less transmitted, the weaker the connections. DGP3 is a block-diagonal system. Such a structure is motivated by e.g. typical quarterly macroeconomic models capturing business cycle dynamic and monetary and fiscal policy effects. One such example is DSGE models, where the dynamic of the economy through time is monitored on quarterly frequency. Note that as written above, DGP1 satisfies the null of no Granger causality from unit 2 to 1, while DGP2 and DGP3 do not. Therefore, we adapt DGP 1 for the power analysis by setting the coefficient in position $(2,1)$ equal to 0.2. Conversely, we set the same coefficient equal to zero for DGP2 and DGP3 for the size analysis.

We choose our series of interest as $I=\{2\}$ and $J=\{1\}$, therby focusing on the case where we have single variables of interest for both elements of the test. Here we consider for simplicity $p=1$ lag, namely the same lag-length as in the DGPs, so $j=1$. The equation of interest can then be written as

[TABLE]

Hence, for each DGP we test $H_{0}:\beta_{GC}=0$ against $H_{1}:\beta_{GC}\neq 0$ using our proposed PDS-LM test.

Table 1 reports the size and power of the test for 1000 replications by using different combinations of time series length $T=(50,100,200,500)$ and number of variables in the system $K=(10,20,50,100)$ and a fixed lag-length $p=1$. All the rejection frequencies are reported using a burn-in period of fifty observations. For each scenario, AIC, BIC and EBIC are compared with the theoretical choice of the tuning parameter $\lambda^{th}$ and time series cross validation $\lambda^{TSCV}$ as described in Subsection 3.3.

Simulations are also reported for different types of covariance matrices of the error terms. We employ a Toepliz-version for calculating the covariance matrix as $\Sigma_{i,j}=\rho^{|i-j|}$ by using two scenarios of correlation: $\rho=(0,0.7)$. The first case corresponds to no correlation, and is equivalent to set $\bm{\varSigma}=\bm{I}_{K}$.

In the Appendix we provide some additional simulation results. First, Table B.1 reports the simulation results for all three DGPs using $\Sigma_{i,j}=0.7^{|i-j|}$. Second, we investigate the Wald version of our test in Table B.2. Third, in Table B.3 we investigate the effects of miss-specification of the lag length by estimating the over-specified VAR$(p+1)$ instead of the true-order VAR$(p)$.555For both the Wald test and the over-specified VAR$(p+1)$ we report the simulations for $\Sigma_{i,j}=0.7^{|i-j|}$ and DGP1 only. Results for the other DGPs are available upon request. Fourth, in Table B we report the results for the size of a bivariate Granger causality test for a non-sparse DGP when using a standard Wald ($F$) test. This test is obviously sensitive to omitted variable bias, and our goal is to demonstrate its effect. Finally, although all results reported here use the finite sample correction in Step 3b of the algorithm, we also investigated the differences with Step 3a. We comment on these results in the next subsection. All results not reported in this paper are available from the authors upon request.

[FIGURE:]

Our proposed approach shows a good performance in terms of size and (unadjusted) power for all DGPs considered. Both for the setting of no correlation and high correlation of errors, sizes are in the vicinity of 5% and power is increasing with the sample size $T$.

Only moderate size distortion is visible in large systems for small samples (e.g. $K\geq 50,\;T=50$). As expected, the test procedure works remarkably well for the sparse DGP1 in high dimensions. However, size properties under the non-sparse DGP2 do not deviate much from its sparse counterpart, although for both DGP2 and DGP3 we do observe a slight deterioration of size when the dimension of the system increases.

Interestingly, the three different information criteria show substantially different behavior. EBIC, due to its very stringent nature, tends to perform well only in very large systems, while it is essentially equal to a bivariate Granger causality test in small systems. We have to add though that the good performance of AIC in particular is somewhat inflated by the imposed lower bound on the penalty; unreported simulations show that without the lower bound AIC performs significantly worse, often selecting too many variables rendering the post-OLS estimation infeasible. The one advantage of using EBIC as information criterion to tune $\lambda$ in the $K>>T$ settings when $T$ is small (e.g. $T=50,100$) is the possibility to avoid the lower bound on the penalty. However, since this comes at a price of more size distortion, we recommend the use of BIC instead, along with the lower bound on the penalty. When comparing the different choices of the tuning parameter we can narrow down the best performing ones (in terms of size and power) to BIC and $\lambda^{th}$. However, in terms of computational time, estimating the tuning parameter using information criteria is considerably faster.

Comparing our test to the bivariate VAR in Table B, it is clear that our proposed PDS-LM is very robust to omitted variable bias, unlike the bivariate test, whose size distortions increase with both the sample size and the number of variables, with sizes of 45% observed for the sample sizes we consider in our application in Section 6. There we will also further elaborate on this difference between our method and the bivariate test. Table 1 shows that for sample sizes smaller than $T=500$, rarely the power exceeds 90%. However, one must keep in mind that the powers are not size-adjusted, and thus the high reported power of the low-dimensional test is an artefact of the huge size distortions rather than genuine power. It also seems unreasonable to expect that PDS-LM test has vey high power if $T$ is small; we are still considering large systems with many parameters to estimate, and there seems to be no way around this if one desires to test Granger causality in large systems with many (control) variables. In that sense we may fully expect the bivariate test to also have higher size-corrected power; yet with all its disadvantages and sensitivity to omitted variables this is not a good comparison. All in all, we believe our test still has sufficiently adequate power properties to be useful in practice.

The results of robustness to misspecification of the lag length order with $p=2$ instead of $p=1$ are reported in Table B.2 in Appendix B. As the size distortions across the range of considered DGPs are only marginally higher for large $K$ and $T$ comparatively small, the test appears to be quite robust to this misspecification. Again, BIC seems to be the best choice for tuning the penalty for all DGPs. Unreported simulations (available upon request) further show that the finite sample adjustment for the test performed in Step 3b of the algorithm is able to substantially reduce size distortions in small samples compared to the asymptotic version of Step 3a.

6 Networks in Realized Volatilities

6.1 Realized Variances

We first investigate the volatility transmission in stock return prices using the daily realized variances of 30 US assets.666We would like to thank Marcelo C. Meideiros for providing us with the high frequency data on stock prices that we have used to construct the realized variances. See Table C.1 for the stocks considered. The R package HDGCvar is available on the GitHub page of the corresponding author (https://github.com/Marga8). Both the computational simplicity and the theoretical foundations make realized volatility measures (realized variance, bi-power variation, median realized variance, etc.) very attractive among practitioners and academics for modelling time varying volatilities and monitoring financial risk. We have considered 10-minute realized variances

[TABLE]

using $j=1,\ldots,M$ intraday 10 minutes stock prices $P_{j,t}$. We consider 10 minute returns as this is the frequency that minimizes for our sample the microstructure noise (McAleer and Medeiros (2008)).777To determine the optimal frequency, we computed realized variances using different frequencies of 1, 5, 10, 15, 30, 65 and 130 minutes, in addition to the estimation using daily returns. The latter estimation has the advantage of being unbiased but the drawback of being very noisy (Pooter et al. (2008)). To find an optimal trade-off between bias and variance (Martens, 2004, see e.g)), mean, variances and mean squared errors (MSE) were computed for each estimation frequency in a similar way as Pooter et al. (2008), and it was found that the frequency of 10 minutes minimizes the MSE. We investigate the period from March 2008 until February 2017 (2236 trading days).

Given the time series of realized volatilities as defined in (9), we employ a multivariate version of the heterogeneous autoregressive model (VHAR) of Corsi (2009) to model their joint behavior (see also Cubadda et al. (2019)). To formally define the VHAR model, we log-transform the series and we stack the logarithmic RV into a vector $y_{t}$. The VHAR specification is given by the following model:

[TABLE]

where $\bm{y}_{t}^{(week)}=\frac{1}{5}\sum_{j=0}^{4}\bm{y}_{t-j}$ and $\bm{y}_{t}^{(month)}=\frac{1}{22}\sum_{j=0}^{21}\bm{y}_{t-j}$ are the vectors containing the average volatility over the last 5 (week) and 22 (month) days. Granger causality in this context represents contagion, or spillover, of volatility from one asset to another. To test for the null hypothesis of no Granger causality / no volatility spillovers from $y_{k,t}$ to $y_{i,t}$ against the alternative of spillovers, we test

[TABLE]

where $\beta_{i,k}^{(1)}$ is the $(i,k)$-th element of $B^{(1)}$. We perform this test for every $(i,k)$-pair to obtain the full $29\times 29$ network of spillover effects. As heteroskedasticity is likely present in these data, we robustify the PDS-LM procedure by implementing the heteroskedasticity-robust LM test such as for example described in Wooldridge (2015, Ch. 8). The full algorithm for the heteroskedasticity-robust PDS-LM test is given in Appendix C.888In the presence of heteroskedasticity, one might prefer the Wald version of the test, as this can be corrected in the standard way by using heteroskedasticty-robust standard errors. Empirically we found hardly any differences between the LM and Wald versions.

We now report the results of our spillover tests for the volatility network. We use BIC to select the tuning parameter of the lasso, and perform the Granger causality tests with a 1% significance level.999We do not perform a correction for multiple testing, as this would only qualitatively affect our results. Moreover, our goal is not to identify exactly the set of spillovers, but to get a feeling of the relations between two variables at a time. As such, we believe a multiple testing correction is not needed, though it can be easily implemented. Figure 1 reports the transmission networks of volatilities estimated with the high-dimensional HVAR (PDS-LM HVAR), bivariate Granger causality tests (BiHVAR) for each pair of stocks, Granger causality tests from a full-system VAR (FullHVAR). The latter is feasible because of our large time series dimension with $T=2236$. For all methods we consider heteroskedasticity-robust variants.

While our PDS-LM method identifies a volatility transmission network which consists of 54 connections and the FullHVAR test picks up 44 connections, the BiHVAR tests detect a network consisting of 803 connections. These differences between our PDS-LM HVAR and the BiHVAR results are in line with our simulation results, confirming that bivariate Granger causality testing in VAR models is seriously affected by omitted variable bias in high-dimensional systems. Given the huge sample size $(T=2236)$ relative to the number of stocks, the FullHVAR is a feasible option, and it is reassuring how similar our PDS-LM HVAR performs compared to the FullHVAR. The similarity is visualized in Figures 1LABEL:sub@fig:im1 and 1LABEL:sub@fig:im3, where the connections picked by both methods are highlighted in red. Of the 54 spillovers identified by the PDS-LM HVAR, 43 are also identified by the FullHVAR, while only 1 of the identified spillovers by the FullHVAR is not picked up by the PDS-LM HVAR.

We also consider a shorter time span, namely the period 2016-2017. Considering a shorter time period makes it more likely that relations remain stable over time. In particular, the chosen period avoids two major events that occurred previously and caused substantial instability on financial markets, namely the global financial crisis of 2008 and the U.S. debt-ceiling crisis of 2011 (Baker et al., 2019). It also allows us to study the performance differences among the three methods in a smaller sample of $T=284$ trading days, where the FullHVAR suffers from the curse of dimensionality. We present the results for the PDS-LM HVAR, BiHVAR and FullHVAR in Figure 2. The number of significant transmissions is 91 for the PDS-LM HVAR, 85 for the BiHVAR and only 5 for the FullHVAR. Hence, the FullHVAR breaks down in this setting due to the small sample size and curse of dimensionality. On the other hand, while superficially the PDS-LM HVAR and the BiHVAR appear to perform similarly, they identify mostly different spillovers. The red lines in Figures 2LABEL:sub@fig:im4 and 2LABEL:sub@fig:im5 show the common connections, which are only 14 out of 91 for the PDS-LM HVAR (85 for the BiHVAR).

As a next step, we use our identified networks to find clusters of closely connected stocks, or communities as they are called in graph theory. Communities are groups of densely connected nodes with fewer connections across groups. In order to represent volatility spillover communities in the graph we use the Newman and Girvan (2004) algorithm based on edge-betweenness. The edge betweenness for edge $e$ is defined as $\sum_{s,t\neq e}\frac{\sigma_{st}(e)}{\sigma_{st}},$ where $\sigma_{st}$ is total number of shortest paths from node $s$ to node $t$ and $\sigma_{st}(e)$ is the number of shortest paths passing through $e$. The edge with the highest betweenness is sequentially removed and the betweenness is recalculated at each step until the best partitioning of the network is found.

Figure 3 reports the graphs of the clustered network for the full sample analysis for the PDS-LM VAR, BiHVAR and FullHVAR respectively. The results for the PDS-LM VAR and FullHVAR show similar spillover clustering behavior, as expected. One large big-industry cluster, containing – among others – assets such as Johnson & Johnson (J&J), IBM, Nike (NKE) and Intel (INTC) dominates the picture being surrounded by small clusters containing 1 to 4 stocks. The PDS-LM VAR and FullHVAR resepctively identify 4 and 6 isolated stocks, which do not have any connections to others. Instead, the BiHVAR finds one single cluster containing all stocks. This reinforces our finding that bivariate Granger causality testing is not informative in high-dimensional systems.

The clusters for the analysis done on the smaller 2016-2017 sample are reported in Figures 4a, 4b and 4c. The patterns highlighted in the spillover network graphs re-occur in the clusters. PDS-LM HVAR in Figure 4a picks up two main clusters of volatility spillovers containing 12 and 6 assets. In addition, four medium size clusters and three single-stock clusters are found. The difference between PDS-LM HVAR and BiHVAR is also reflected in the identified clusters. BiHVAR in Figure 4b shows only one big conglomerate cluster of stocks linked to three two-stock clusters and 6 single-stock clusters. Finally, the breakdown of FullHVAR shows clearly in the non-informative, mostly unconnected single-stock clusters in Figure 4c.

6.2 Realized Variances & Covariances

In this subsection we extend our investigation to allow for spillovers from realized correlations to variances. While our application in Section 6.1 was only high-dimensional when we considered the shorter subsample, including correlations, which are of the order $K^{2}$, put a significantly larger strain on estimation, making the standard full VAR no option. As elaborated later, it appears quite reasonable to expect changing correlations to also have an affect on the volatilities. By ignoring these in Section 6.1, we are exposing ourselves again to a potential omitted variable bias. However, our method can directly incorporate these, as we demonstrate here.

While we remain mostly interested in contagion between the 30 realized volatilities, we add the $\frac{30\times 29}{2}=435$ realized correlations between all these assets as control variables. By maintaining our focus on the relations between the variances, our results are directly comparable to Section 6.1 and can be interpreted by assessing how the network changes when the correlations are added as controls in the VAR. Moreover, it also avoids the difficulties of trying to visualize the results from all $(30\times 435)^{2}$ possible connections in the large VAR considered here. In the same way that the realized variances employ high frequency data to estimate the integrated variance, the realized covariance (RC hereafter) estimates the integrated covariance of a multivariate diffusion process. Working with the full RC time-varying matrix is important for portfolio allocation and risk management. For a set of $n$ intra-day asset returns at day $t$ observed at $j=1,\ldots,M$ stacked in a column vector $\bm{r}_{j,t}$, the realized covariance is obtained such as $\bm{RC}_{t}=\sum_{j=1}^{M}\bm{r}_{j,t}\bm{r}_{j,t}^{\prime}$. Note that the realized variances are obviously on the diagonal of RC and that the RC matrix is positive definite when $M>n,$ namely when the number of high frequency observations per day is larger than the number of series. We investigate the same period as before and construct 10-minute realized covariances. Several studies have also proposed a Lasso framework on RC, see for instance Callot et al. (2017) and Brito et al. (2018), although their focus is more on portfolio allocation and forecasting.

There are two main ways to work with the RC matrix. The first approach stacks realized variances and covariances in a single vector. For instance, Bauer and Vorkink (2011) consider the matrix log transformation of $RC_{t}$ series, a matrix that they call the log-space volatility. The drawback of that log transform is that the interpretation of the original series, in our case the volatilities, is lost as the combinations involve nonlinear transforms of both realized variances and covariances. This is not compatible with the aim of this paper.

The second approach uses the log realized volatilities and the correlations separately, as done by for instance Oh and Patton (2016). The underlying idea, following the DCC model of Engle (2002), is to decompose $RC_{t}^{(d)}=D_{t}^{(d)}R_{t}^{(d)}D_{t}^{(d)}$ with $D_{t}^{(d)}$ a diagonal matrix with the square root of the individual realized variance and $R_{t}^{(d)}$ the realized correlation matrix. Oh and Patton (2016) use the HVAR model structure for each realized volatilities, they consequently assume no Granger causality across volatilities.

We propose something which is, to some extent, in between these two approaches. We look at two separate objects as in the DCC model, but stack the log of the realized variances $\bm{y}_{1t}^{(d)^{\prime}}$ and $z$-transforms $\bm{y}_{2t}^{(d)}=\mathrm{arc}\tanh\left(\mathrm{vech}(\bm{R}_{t}^{(d)})\right)$ of the realized correlations in a larger vector $\bm{y}_{t}^{(d)}=(\underset{1\times 30}{\underbrace{\bm{y}_{1t}^{(d)^{\prime}}}},\underset{1\times 435}{\underbrace{\bm{y}_{2t}^{(d)^{\prime}}}})^{\prime}$ on which we estimate a VHAR of dimension 465.

In this HVAR each of the 465 equations depends on 1395 dynamic parameters plus the constant. We focus on the 30 equations corresponding to $\bm{y}_{1t}^{(d)}$ volatilities and consequently the bivariate causalities between these realized volatilities as in the previous section. Figure 5a reports a total of 113 connections, which is about twice the connections in Figure 1a. In red we highlighted the 31 common connections with Figure 1a. Interestingly, adding more variables therefore allows us to uncover more relations. It seems that this allows us to uncover partial effects that were previously obscured by counteracting effects of the correlations. Importantly, the number of connections is still far less than compared to the BiHVAR in Figure 1b, and the PDS-LM HVAR is still able to deliver a clear picture of the causal connections when the system considered is high-dimensional. While the different connections found here obviously also lead to a different clustering, Figure 5b shows that the clustering is quite similar, certainly regarding qualitative conclusions.

7 Conclusion

We propose an LM test in order to test for Granger causality in high-dimensional VAR models. We employ a post-double selection procedure using the lasso to select the set of relevant covariates in the system. The double selection step allows to substantially reduce the omitted variable bias and thereby allowing for valid post-selection inference on the parameters.

We provide an extensive simulation study to evaluate the performance of our method in finite samples, paying particular attention to the tuning of the penalty parameter. We compare different information criteria, time series cross-validation and a plug-in method based on theoretical arguments, and find that generally BIC and the theoretically tuned penalty perform best. However, to use information criteria in systems with a significantly larger number of variables than observations, a lower bound on the penalty parameter is needed to prevent too many variables being selected. The simulations also show that, when properly tuned, our proposed PDS-LM test attains good results both for size and power under different DGPs. Especially, it is shown to be robust both to non-sparse settings as well as to lag-length overspecification.

We also empirically investigate the usefulness of our method in a study where we apply our PDS-LM method to a high-dimensional VHAR process in order to construct a contagion network of volatility spillovers for 30 large capital stocks, also accounting for effects from changing correlations. We find that by increasing the information set through considering a high-dimensional VAR model instead of bivariate models, we are able to obtain more realistic effects than in low-dimensional models. Furthermore, even when the sample size is not large enough to use standard full-system VAR techniques, our method remains reliable and delivers accurate results.

Note that unlike Belloni et al. (2014a), we do not give a “truly” causal interpretation to the established Granger causalities. In how far Granger causality is a useful concept to study true causality is (and has long been) open to debate, see for example (Eichler, 2013) and the references therein. Moreover, though it appears desirable and in line with Granger’s (1969) original intentions to make the information set as large as possible, it is well known in the literature on graphical models (see Eichler, 2013) for causality that considering only the full model is not sufficient for establishing true causal relations from Granger causal ones. For instance, one-period Granger causality in systems with more than two variables cannot capture indirect causal chains spanning over multiple periods. However, the analysis of the full model is a necessary ingredient for any study of causality in a graphical framework. It would therefore be an interesting avenue for further research to study how the method proposed here could fit into such a graphical framework.

Appendix A Proofs

Lemma A.1**.**

Let $\bm{X}_{-GC}$ satisfy Assumption 2(g). Then with probability at least $1-\Delta_{T}$, we have that

[TABLE]

for any $\bm{\delta}=(\bm{\delta}_{1}^{\prime},\ldots,\bm{\delta}_{N_{I}}^{\prime})^{\prime}$ such that ${\left\lvert S_{\delta}\right\rvert}\leq\bar{s}_{T}$, where $S_{\delta}=\bigcup_{i=1}^{N_{I}}\{m:\delta_{i,m}\neq 0\}$.

Proof of Lemma A.1.

As before, let $\bm{X}_{S}$ denote the submatrix containing those columns of $\bm{X}_{-GC}$ corresponding to the elements in $S$. It follows from Assumption 2(a) that for any $\bm{\gamma}$ satisfying ${\left\lvert S_{\bm{\gamma}}\right\rvert}\leq\bar{s}_{T}$, we have that $\lambda_{\min}(\bm{X}_{S_{\bm{\gamma}}}^{\prime}\bm{X}_{S_{\bm{\gamma}}}/T)\geq\phi_{T,\min}^{2}$. Then, with probability $1-\Delta_{T}$ we have that

[TABLE]

as $\lambda_{\min}(\bm{G}_{T}^{\prime}\bm{G}_{T})\geq C>0$ by assumption. Without loss of generality we may then absorb the constant $C$ into $\phi_{T,\min}$. ∎

Proof of Theorem 1.

Our proof follows along the lines of the proof of Theorem 1 and 2 of Belloni et al. (2014b), with the main distinction that we consider multiple variables of interest, and multiple “treatments” instead of a single one for each. For this purpose we first define some notation. Let $\bm{\varGamma}=(\bm{\gamma}_{1},\ldots,\bm{\gamma}_{N_{X}})$ and $\hat{\bm{\varGamma}}=(\hat{\bm{\gamma}}_{1},\ldots,\hat{\bm{\gamma}}_{N_{X}})$, and let $\bm{\gamma}^{\otimes}=\bm{I}_{N_{I}}\otimes\bm{\varGamma}$. Furthermore, let $\mathcal{P}(\bm{A})=\bm{A}(\bm{A}^{\prime}\bm{A})^{-1}\bm{A}^{\prime}$ denote the projection on the space spanned by $\bm{A}$ and let $\mathcal{M}(\bm{A})=\bm{I}-\mathcal{P}(\bm{A})$ denote the corresponding residual-maker. By standard partitioned regression algebra we get

[TABLE]

where $\widetilde{\bm{X}}^{\otimes}=\bm{G}_{T}^{\otimes}\bm{X}^{\otimes}=\bm{G}_{T}\otimes\bm{X}$. We will now show that $\bm{a}_{T}=\widetilde{\bm{E}}^{\otimes\prime}\widetilde{\bm{u}}_{I}/\sqrt{T}$ and $\bm{B}_{T}=\widetilde{\bm{E}}^{\otimes\prime}\widetilde{\bm{E}}^{\otimes}/T+o_{p}(1)$. Given these results, the limit distribution then follows directly from Assumption 2(b).

We first consider $\bm{a}_{T}$. Note that from (7) we have that $\widetilde{\bm{X}}_{GC}^{\otimes}=\bm{G}_{T}\otimes\left[\bm{X}_{-GC}\bm{\varGamma}+\bm{E}\right]=\widetilde{\bm{X}}_{-GC}^{\otimes}\bm{\varGamma}^{\otimes}+\widetilde{\bm{E}}^{\otimes}$, and therefore we can write

[TABLE]

We will now show that the terms $\bm{a}_{T,1},\ldots,\bm{a}_{T,4}$ vanish. For $\bm{a}_{T,1}$, note that

[TABLE]

where for any matrix $\bm{M}$, the norm ${\left\lVert\cdot\right\rVert}_{p}$ represents the induced $l_{p}$-matrix norm ${\left\lVert\bm{M}\right\rVert}_{p}=\sup_{\bm{x}\neq 0}{\left\lVert\bm{M}\bm{x}\right\rVert}_{p}/{\left\lVert\bm{x}\right\rVert}_{p}$. As $\hat{S}_{j}\subseteq\hat{S}_{X}$ for all $j=1,\ldots,N_{X}$ and $\hat{S}_{X}^{\otimes}\subseteq\hat{S}^{\otimes}$, the space spanned by $\bm{G}_{T}^{\otimes}\bm{X}_{\hat{S}_{X}^{\otimes}}=\bm{G}_{T}\otimes\bm{x}_{\hat{S}_{X}}$ is a subspace of the space spanned by $\bm{G}_{T}^{\otimes}\bm{X}_{\hat{S}^{\otimes}}^{\otimes}$, and therefore ${\left\lVert\mathcal{M}\left(\bm{G}_{T}^{\otimes}\bm{X}_{\hat{S}^{\otimes}}^{\otimes}\right)\bm{y}\right\rVert}_{2}\leq{\left\lVert\mathcal{M}\left(\bm{G}_{T}^{\otimes}\bm{X}_{\hat{S}_{X}^{\otimes}}^{\otimes}\right)\bm{y}\right\rVert}_{2}$ for any compatible matrix $\bm{G}_{T}^{\otimes}$ and vector $\bm{y}$. Then, using that

[TABLE]

we find that

[TABLE]

Then,

[TABLE]

as $\hat{S}_{j}=\{m:\hat{\gamma}_{m,j}\neq 0\}$ and therefore the constraint in the minimization is satisfied. It then follows from Assumption 2(e) that ${\left\lVert\bm{A}_{T,1,1}\right\rVert}_{2}\leq N_{GC}\delta_{T}T^{-1/4}$ with probability $1-\Delta_{T}$.

For $\bm{a}_{T,1,2}$, from the definition of the best linear predictor it directly follows that

[TABLE]

such that we can substitute $\bm{\beta}_{-GC}=\bm{\gamma}_{0}-\bm{\varGamma}^{\otimes}\bm{\beta}_{GC}$ in $\bm{a}_{T,1,2}$ to find

[TABLE]

where the negligibility of the second term follows directly from the result above plus Assumption 2(d). As $\hat{S}_{0}\subseteq\hat{S}^{\otimes}$, the first term can be bounded by

[TABLE]

where we use that, for $\hat{\bm{\beta}}_{\bm{G}_{T}}=\mathcal{M}(\bm{G}_{T}\bm{X})\bm{G}_{T}\bm{y}$,

[TABLE]

It then follows directly that ${\left\lVert\bm{a}_{T,1}\right\rVert}_{2}=O_{p}(\delta_{T}^{2})=o_{p}(1)$.

For $\bm{a}_{T,2}$, let $\bm{\gamma}_{j}$ denote the $j$-th column of $\bm{\varGamma}^{\otimes}$ and define the noiseless generalized least squares estimator

[TABLE]

for any compatible index set $S$, and let $\check{\bm{\varGamma}}_{S}^{\otimes}=\left(\check{\bm{\gamma}}_{1,S}^{\otimes},\ldots,\check{\bm{\gamma}}_{N_{X},S}^{\otimes}\right)$, such that $\mathcal{M}\left(\widetilde{\bm{X}}_{S}^{\otimes}\right)\widetilde{\bm{X}}_{-GC}^{\otimes}\bm{\varGamma}^{\otimes}=\widetilde{\bm{X}}_{-GC}^{\otimes}\left(\bm{\varGamma}^{\otimes}-\check{\bm{\varGamma}}_{S}^{\otimes}\right)$. Then, with probability $1-\Delta_{T}$,

[TABLE]

Here inequality $(i)$ uses that

[TABLE]

from the dual norm inequality, where $\bm{A}$ is a generic $m\times n$ matrix $\bm{A}$ with $j$-th row denoted as $\bm{a}_{j\cdot}$, and a $n\times 1$ vector $\bm{x}$. Letting ${\left\lVert\bm{A}\right\rVert}_{\max}=\max_{i,j}{\left\lvert a_{ij}\right\rvert}$, Step $(ii)$ follows from the fact that ${\left\lVert\widetilde{\bm{X}}_{-GC}^{\otimes\prime}\widetilde{\bm{u}}_{I}/\sqrt{T}\right\rVert}_{\infty}={\left\lVert(\bm{G}_{T}^{\prime}\bm{G}_{T}\otimes\bm{X}_{-GC}^{\prime})\bm{u}_{I}/\sqrt{T}\right\rVert}_{\infty}\leq{\left\lVert\bm{G}_{T}^{\prime}\bm{G}_{T}\right\rVert}_{\max}{\left\lVert\bm{X}_{-GC}^{\otimes\prime}\bm{u}_{I}/\sqrt{T}\right\rVert}_{\infty}\leq\gamma_{T}$ by Assumption 2(c), while $(iii)$ follows from bounding the $l_{1}$-norm by the $l_{2}$-norm and applying Lemma A.1. $(iv)$ follows from the definition of $\check{\bm{\gamma}}_{\hat{S}}$ as minimizer of the sum of squares and $(v)$ from the properties of the Kronecker product. Finally $(vi)$ follows from Assumption 2(e).

For $\bm{a}_{T,3}$, define $\check{\bm{\gamma}}_{0,S}=\operatorname*{arg\,min}_{\bm{\gamma}:\gamma_{m}=0,m\notin S}{\left\lVert\widetilde{\bm{X}}_{-GC}^{\otimes}\bm{\gamma}_{0}-\widetilde{\bm{X}}_{-GC}^{\otimes}\bm{\gamma}\right\rVert}_{2}^{2}$ analogously to (A.4). Then we have with probability $1-\Delta_{T}$

[TABLE]

Inequality $(i)$ follows from the fact that $\bm{\beta}_{-GC}=\bm{\gamma}_{0}-\bm{\varGamma}^{\otimes}\bm{\beta}_{GC}$, while $(ii)$ follows from the definition of $\tilde{\bm{\gamma}}_{0,S}$ and (A.4). For the first term in (iii) we use (A.5) whereas for the second term we apply it twice to get

[TABLE]

for any $p\times n$ matrix $\bm{B}$. Step $(iv)$ follows from Assumption 2(d) and the results for $\bm{a}_{T,2}$, while Step $(v)$ applies the same arguments as used therein to ${\left\lVert\check{\bm{\gamma}}_{0,\hat{S}}-\bm{\gamma}_{0}\right\rVert}_{1}$. Finally, Step $(vi)$ follows analogoulsy to Step $(ii)$ for $\bm{a}_{T,2}$ by noting that ${\left\lVert\widetilde{\bm{X}}_{-GC}^{\otimes\prime}\widetilde{\bm{e}}_{j}^{\otimes}/\sqrt{T}\right\rVert}_{\infty}\leq{\left\lVert\bm{G}_{T}^{\prime}\bm{G}_{T}\right\rVert}_{\max}{\left\lVert\bm{X}_{-GC}^{\prime}\bm{e}_{j}/\sqrt{T}\right\rVert}_{\infty}\leq\gamma_{T}$, plus using the bound from Assumption 2c).

Finally, we consider $\bm{a}_{T,4}$. We get

[TABLE]

where step $(i)$ follows from (A.5). For $(ii)$ we bound the $l_{1}$-norm with the $l_{2}$-norm, using that $\widetilde{\bm{X}}_{\hat{S}^{\otimes}}^{\otimes}$ contains a subset of the variables in $\widetilde{\bm{X}}_{-GC}^{\otimes}$ and therefore ${\left\lVert\widetilde{\bm{X}}_{\hat{S}^{\otimes}}^{\otimes\prime}\widetilde{\bm{u}}_{I}^{\otimes}\right\rVert}_{\infty}\leq{\left\lVert\widetilde{\bm{X}}_{-GC}^{\otimes\prime}\widetilde{\bm{u}}_{I}^{\otimes}\right\rVert}_{\infty}$ and apply Assumption 2(c). Step $(iii)$ follows from the Cauchy-Schwarz inequality, bounding the $l_{2}$-norm by the $l_{\infty}$-norm and reasoning as for Step $(ii)$ that ${\left\lVert\widetilde{\bm{X}}_{\hat{S}^{\otimes}}^{\otimes\prime}\widetilde{\bm{e}}_{j}^{\otimes}\right\rVert}_{\infty}\leq{\left\lVert\widetilde{\bm{X}}_{-GC}^{\otimes\prime}\widetilde{\bm{e}}_{j}^{\otimes}\right\rVert}_{\infty}$. Finally $(iv)$ follows from Assumption 2(c) and Lemma A.1.

We next consider $\bm{B}_{T}$. Using that $\widetilde{\bm{X}}_{GC}^{\otimes}=\widetilde{\bm{X}}_{-GC}^{\otimes}\bm{\varGamma}^{\otimes}+\widetilde{\bm{E}}^{\otimes}$, we write

[TABLE]

These terms can be handled as the terms for $\bm{a}_{T}$. In particular, with probability $1-\Delta_{T}$, ${\left\lVert\bm{B}_{T,1}\right\rVert}_{2}\leq{\left\lVert\bm{A}_{T,1,1}\right\rVert}_{2}^{2}\leq\delta_{T}^{2}T^{-1/2}$, ${\left\lVert\bm{B}_{T,2}\right\rVert}_{2}\leq\delta_{T}^{2}T^{-1/2}$ using the same steps as for $\bm{a}_{T,2}$, and ${\left\lVert\bm{B}_{T,3}\right\rVert}_{2}\leq\delta_{T}T^{-1/2}$ analogously to $\bm{a}_{T,4}$.

This shows that $\bm{a}_{T}=\widetilde{\bm{E}}^{\otimes\prime}\widetilde{\bm{u}}_{I}/\sqrt{T}$ and $\bm{B}_{T}=\widetilde{\bm{E}}^{\otimes\prime}\widetilde{\bm{E}}^{\otimes}/T+o_{p}(1)$. It then follows directly from Assumption 2(b) that

[TABLE]

Proof of Theorem 2.

By partitioned regression algebra, we find that

[TABLE]

Note that $\bm{a}_{T}^{*}$ and $\bm{B}_{T}^{*}$ are special cases of their counterparts in the proof of Theorem 1 with $\bm{G}_{T}=\hat{\bm{\varSigma}}_{T,I}^{-1/2}$. We now show that this choice of $\bm{G}_{T}$ satisfies the conditions of Theorem 1. We do this by proving that $\bm{G}_{T}$ converges to $\bm{G}=\bm{\varSigma}_{u,I}^{-1/2}$, and this satisfies the conditions in the theorem.

Consider one particular element $(i,j)$ of $\hat{\bm{\varSigma}}_{u,I}$, say $\hat{\sigma}_{u,ij}$. Let $\hat{S}_{0,I_{i}}$ denote the variables selected in $\hat{S}_{0}$ corresponding to the equation for variable $\bm{y}_{I_{i}}$, where $I=\{I_{1},\ldots,I_{N_{I}}\}$, and let $\hat{S}_{i}=\left(\bigcup_{j=1}^{N_{X}}\hat{S}_{j}\right)\cup\hat{S}_{0,i}$ denote all variables selected that are relevant for $\bm{y}_{I_{i}}$. We can then write

[TABLE]

where (under $H_{0}$) we write $\bm{y}_{I_{i}}=\bm{X}_{-GC}\bm{\beta}_{-GC,i}+\bm{u}_{I_{i}}$.

We can use the same reasoning as used in the proof of Theorem 1 to prove that the terms $d_{T,ij,k}$, $k=1,2,3$ are negligible. Let $\gamma_{0,i}$ denote the sub-vector of $\bm{\gamma}_{0}$ corresponding to unit $i$, and note that under the null hypothesis $\bm{\gamma}_{0,i}=\bm{\beta}_{-GC,i}$. Define $\bar{\bm{\gamma}}_{i,S}=\operatorname*{arg\,min}_{\bm{\gamma}:\gamma_{m}=0,m\notin S}{\left\lVert\bm{X}_{-GC}\gamma_{0,i}-\bm{X}_{-GC}\bm{\gamma}\right\rVert}_{2}^{2}$. Then, as $\hat{S}_{0,i}\subseteq\hat{S}_{i}$, we get that

[TABLE]

Similarly,

[TABLE]

Finally,

[TABLE]

where all bounds hold with probability at least $1-\Delta_{T}$. Similarly, by Assumption 2(b) we know that there exist a sequence $\delta_{T}\rightarrow 0$, such that with probability $1-\Delta_{T}$, we have that ${\left\lvert\bm{u}_{I_{i}}^{\prime}\bm{u}_{I_{j}}-\sigma_{u,ij}\right\rvert}\leq\delta_{T}$. As $\bm{\varSigma}{u,I}$ only contains a finite number ($N_{I}^{2}$) elements, we may then conclude that with probability at least $1-\Delta_{T}$, it holds that ${\left\lVert\hat{\bm{\varSigma}}_{u,I}-\bm{\varSigma}\right\rVert}_{2}\leq\delta_{T}$.

We have that $0<c_{L}\leq\lambda_{\min}(\bm{\varSigma}_{u,I})\leq\lambda_{\max}(\bm{\varSigma}_{u,I})\leq c_{u}<\infty$, where the lower bound follows from Assumption 2(a) and the upper bound from the fact that $\bm{\varSigma}_{u,I}$ has a finite number of elements. As

[TABLE]

the established the consistency of $\hat{\bm{\varSigma}}_{u,I}$ then directly yields that $C_{1}-\delta_{T}\leq\lambda_{\min}(\hat{\bm{\varSigma}}_{u,I})\leq\lambda_{\min}(\hat{\bm{\varSigma}}_{u,I})\leq C_{2}+\delta_{T}$ with probability $1-\Delta_{T}$. It then also follows that with probability $1-\Delta_{T}$ we can find a $0<C_{1}\leq C_{2}<\infty$ such that $c_{1}\leq 1/\lambda_{\max}(\hat{\bm{\varSigma}}_{u,I})=\lambda_{\min}(\hat{\bm{\varSigma}}_{u,I}^{-1})\leq\lambda_{\max}(\hat{\bm{\varSigma}}_{u,I}^{-1})=1/\lambda_{\min}(\hat{\bm{\varSigma}}_{u,I})\leq c_{2}$, such that the conditions of Theorem 1 are satisfied for $\bm{G}_{T}=\hat{\bm{\varSigma}}_{u,I}^{-1/2}$.

With this choice of $\bm{G}_{T}$ and $\bm{\varOmega}=\bm{\varSigma}_{u,I}\otimes\bm{\varSigma}_{GC|-GC}$, we have that $\bm{\varOmega}_{\bm{G}}=\bm{I}_{N_{X}}\otimes\bm{\varSigma}_{GC|-GC}$. Letting $\bm{Z}_{N_{GC}}\sim\mathcal{N}(\bm{0},\bm{I}_{N_{GC}})$, it then follows that

[TABLE]

Appendix B Additional Simulation Results

[FIGURE:]

[FIGURE:]

[FIGURE:]

Simulation results for the bivariate Granger causality test

DGP Size/Power $\rho$ ${K\char 92\relax T}$ 50 100 200 500

10 5.9 6.6 7.8 11.8

2 Size 0 20 5.6 5.9 7.8 11.8

50 4.3 7.0 9.7 14.5

100 5.5 6.7 8.9 13.9

• •

Notes: Size is reported for DGP 2, as described in Section 5, for 1000 replications. $T=(50,100,200,500)$ is the time series length, $K=(10,20,50,100)$ the number of variables in the system, the lag length is fixed to $p=1$. $\rho$ indicates the correlation employed to simulate the time series with the Toeplitz covariance matrix.

Appendix C Additional Material for the Empirical Application