High-dimensional GARCH process segmentation with an application to Value-at-Risk

TL;DR

This paper introduces a method for detecting multiple change points in high-dimensional GARCH models, improving Value-at-Risk estimation by accounting for structural breaks in both individual volatilities and correlations.

Contribution

It proposes a novel consistent change point detection method for high-dimensional GARCH panel data, addressing both individual and cross-sectional dependence changes.

Findings

Method demonstrates high accuracy in simulations

Effective in real-world Value-at-Risk application

Available as an R package segMGarch

Abstract

Models for financial risk often assume that underlying asset returns are stationary. However, there is strong evidence that multivariate financial time series entail changes not only in their within-series dependence structure, but also in the cross-sectional dependence among them. In particular, the stressed Value-at-Risk of a portfolio, a popularly adopted measure of market risk, cannot be gauged adequately unless such structural breaks are taken into account in its estimation. We propose a method for consistent detection of multiple change points in high-dimensional GARCH panel data set where both individual GARCH processes and their correlations are allowed to change over time. We prove its consistency in multiple change point estimation, and demonstrate its good performance through simulation studies and an application to the Value-at-Risk problem on a real dataset. Our methodology is implemented in the R package segMGarch, available from CRAN.

Figures

Tables

Table 1. Table 1: (M0) Size of the change point test at α = 0.05 𝛼 0.05 \alpha=0.05 when T = 1000 𝑇 1000 T=1000 .

	Gaussian ${\epsilon }_{i,t}$		$t_{10}$-distributed ${\epsilon }_{i,t}$
$N$	(M0)(M0.1)	(M0)(M0.2)	(M0)(M0.1)	(M0)(M0.2)
$50$	0.01	0.05	0.03	0.13
$100$	0.02	0.09	0.02	0.3

Table 2. Table 2: (M1)  (i) The number of estimated change points (%) and the accuracy in change point location (% of | η ^ 1 − η 1 | < log 2 ⁡ T subscript ^ 𝜂 1 subscript 𝜂 1 superscript 2 𝑇 |\widehat{\eta}_{1}-\eta_{1}|<\log^{2}T ) when α = 0.05 𝛼 0.05 \alpha=0.05 , T = 500 𝑇 500 T=500 , N = 50 𝑁 50 N=50 (left) and N = 100 𝑁 100 N=100 (right) with Gaussian innovations.

		$\hat{B}$					accuracy (%)		$\hat{B}$					accuracy (%)
	$\varrho$	0	1	2	3	$\ge 4$	${\eta }_1$	${\eta }_2$	0	1	2	3	$\ge 4$	${\eta }_1$	${\eta }_2$
(M1)(M1.1)	1	0	0	100	0	0	100	100	0	0	100	0	0	90	99
	0.75	0	0	100	0	0	100	100	0	0	100	0	0	100	100
	0.5	0	0	100	0	0	100	100	0	0	100	0	0	97	99
	0.25	0	7	92	1	0	100	93	0	91	9	0	0	100	6
(M1)(M1.2)	1	0	0	100	0	0	100	100	0	0	100	0	0	100	100
	0.75	0	2	98	0	0	96	100	0	0	100	0	0	100	100
	0.5	0	38	62	0	0	60	100	0	2	98	0	0	98	100
	0.25	0	95	5	0	0	2	100	0	80	20	0	0	32	84
(M1)(M1.3)	1	0	89	11	0	0	3	100	0	68	31	1	0	13	100
	0.75	0	89	11	0	0	3	100	0	85	15	0	0	3	100
	0.5	0	94	6	0	0	0	100	0	88	12	0	0	3	100
	0.25	0	97	3	0	0	0	100	1	96	3	0	0	0	98

Table 3. Table 3: (M1)  (ii) The number of estimated change points (%) and the accuracy in change point location when α = 0.05 𝛼 0.05 \alpha=0.05 , T = 500 𝑇 500 T=500 , N = 50 𝑁 50 N=50 (left) and N = 100 𝑁 100 N=100 (right) with t 10 subscript 𝑡 10 t_{10} innovations.

		$\hat{B}$					accuracy (%)		$\hat{B}$					accuracy (%)
	$\varrho$	0	1	2	3	$\ge 4$	${\eta }_1$	${\eta }_2$	0	1	2	3	$\ge 4$	${\eta }_1$	${\eta }_2$
(M1)(M1.1)	1	0	5	74	12	8	92	100	0	7	72	11	10	90	100
	0.75	0	2	86	8	4	98	100	1	1	67	20	11	96	99
	0.5	0	2	78	14	6	98	100	3	1	72	12	12	97	100
	0.25	0	28	58	6	8	100	70	1	47	39	8	5	99	40
(M1)(M1.2)	1	0	31	67	2	0	67	100	0	0	100	0	0	90	99
	0.75	0	11	86	3	0	86	100	0	0	100	0	0	100	100
	0.5	0	5	93	2	0	100	100	0	5	93	2	0	94	100
	0.25	0	3	94	2	0	94	99	0	91	9	0	0	100	6
(M1)(M1.3)	0	0	99	1	0	0	0	100	0	100	0	0	0	0	100
	0.75	0	99	1	0	0	0	100	0	100	0	0	0	0	100
	0.5	0	100	0	0	0	0	100	0	100	0	0	0	0	100
	0.25	0	99	1	0	0	0	100	1	100	0	0	0	0	100

Table 4. Table 4: (M2)(M2.1) The number of estimated change points (%) and the accuracy in change point location when α = 0.05 𝛼 0.05 \alpha=0.05 , T = 500 𝑇 500 T=500 , N = 50 𝑁 50 N=50 (left) and N = 100 𝑁 100 N=100 (right).

		$\hat{B}$					accuracy (%)		$\hat{B}$					accuracy (%)
	$\varrho$	0	1	2	3	$\ge 4$	${\eta }_1$	${\eta }_2$	0	1	2	3	$\ge 4$	${\eta }_1$	${\eta }_2$
(M2.1.1)	1	0	0	98	2	0	100	100	0	0	100	0	0	100	100
	0.75	0	0	100	0	0	100	100	0	0	99	1	0	100	100
	0.5	0	0	95	5	0	100	100	0	0	100	0	0	100	100
	0.25	0	18	73	8	1	100	80	0	97	3	0	0	100	1
(M2.1.2)	1	0	1	97	2	0	99	100	0	0	99	1	0	100	100
	0.75	0	8	90	2	0	92	100	0	0	99	1	0	100	100
	0.5	0	50	47	3	0	46	100	0	7	85	8	0	93	100
	0.25	1	88	10	1	0	5	99	3	79	18	0	0	27	78

Table 5. Table 5: (M2)(M2.2) The number of estimated change points (%) and the accuracy in change point location when α = 0.05 𝛼 0.05 \alpha=0.05 , T = 500 𝑇 500 T=500 , N = 50 𝑁 50 N=50 (left) and N = 100 𝑁 100 N=100 (right).

		$\hat{B}$					accuracy (%)		$\hat{B}$					accuracy (%)
	$\varrho$	0	1	2	3	$\ge 4$	${\eta }_1$	${\eta }_2$	0	1	2	3	$\ge 4$	${\eta }_1$	${\eta }_2$
(M2)((M2.2))(M2.2.1)	1	0	0	100	0	0	100	100	0	0	100	0	0	100	100
	0.75	0	0	98	2	0	100	100	0	0	98	2	0	100	100
	0.5	0	0	99	1	0	100	100	0	0	99	1	0	100	100
	0.25	0	5	92	2	1	95	100	0	15	74	11	0	100	83
(M2)((M2.2))(M2.2.2)	1	0	0	99	1	0	100	100	0	0	100	0	0	100	100
	0.75	0	0	97	3	0	100	100	0	0	99	1	0	100	100
	0.5	0	0	98	2	0	100	100	0	0	98	2	0	100	100
	0.25	0	45	51	4	0	48	100	0	12	79	9	0	92	93

Table 6. Table 6: (M3) The number of estimated change points (%) and the accuracy in change point location when α = 0.05 𝛼 0.05 \alpha=0.05 , T = 500 𝑇 500 T=500 , N = 50 𝑁 50 N=50 (left) and N = 100 𝑁 100 N=100 (right).

		$\hat{B}$					accuracy (%)		$\hat{B}$					accuracy (%)
	$\varrho$	0	1	2	3	$\ge 4$	${\eta }_1$	${\eta }_2$	0	1	2	3	$\ge 4$	${\eta }_1$	${\eta }_2$
(M3)(M3.1)	1	0	5	34	52	9	91	97	0	0	45	44	11	95	97
	0.75	0	3	36	49	12	95	90	0	0	33	60	7	95	96
	0.5	0	4	37	52	7	89	89	0	3	17	64	16	96	88
	0.25	0	14	42	39	5	77	53	0	6	32	48	14	92	70
(M3)(M3.2)	1	0	13	87	0	0	87	99	0	4	96	0	0	96	98
	0.75	0	22	78	0	0	79	96	0	7	93	0	0	93	97
	0.5	0	74	24	2	0	29	93	0	38	61	1	0	64	97
	0.25	6	89	5	0	0	8	70	0	93	7	0	0	8	92

Table 7. Table 7: Three change points detected by our methodology (in bold) from the training data spanning from 1 January 2007 to 28 April 2016, and the corresponding periods of stationarity (Periods 1–4). Period 5 covers the 12-month period following the bankruptcy of Lehman Brothers with its length chosen in accordance with Basel 2.5. It also presents the single change point detected by our methodology (in italic) from the data panning from 6 August 2009 to 31 December 2020 and the corresponding periods of stationarity (Periods 6 and 7). For each period, we report the 95% and 99% VaR of an equally weighted portfolio of the N = 79 𝑁 79 N=79 stocks.

	Period range (mm/dd/yyyy)		Value-at-Risk
	From	To	95%	99%
Period 1	01/01/2007	09/07/2008	0.01836	0.02600
Period 2	09/08/2008	05/04/2009	0.05326	0.07508
Period 3	05/05/2009	08/05/2009	0.01904	0.02765
Period 4	08/06/2009	04/28/2016	0.03217	0.04554
Period 5	09/08/2008	09/08/2009	0.04416	0.06254
Period 6	08/06/2009	02/24/2020	0.02523	0.03880
Period 7	02/25/2020	12/31/2020	0.03711	0.05560

Table 8. Table 8: Results from sVaR backtesting on the test data (from 1 May 2016 to 31 December 2020) using the historical simulation approach calibrated on five different stress periods b = 1 , … , 5 𝑏 1 … 5 b=1,\ldots,5 from Table  7 . From left to right, the table reports: The number of sVaR failures (PoF) and the number of days until the first failure in sVaR (TFF), with the corresponding p 𝑝 p -values in brackets from the respective two-sided tests in ( 13 ); traffic light test (traffic); p 𝑝 p -values from the DQ test (DQ) in ( 14 ).

Scenario			5-day rebalance ($h_{\text{reb}}=5$)				20-day rebalance ($h_{\text{reb}}=20$)
${\mu }_{\text{target}}$	$T_{\text{sv}}$	Period	PoF	TFF	traffic	DQ	PoF	TFF	traffic	DQ
0.25%	250	1	17 (0.003)	14 (0.132)	yellow(0.014)	0	13 (0.244)	195 (0.451)	green(0.143)	0
		2	0 (1)	NA (1)	green(1)	0.157	0 (1)	NA (1)	green(1)	0.157
		3	20 (0.002)	195 (0.451)	yellow(0.001)	0	17 (0.022)	195 (0.451)	yellow(0.014)	0
		4	7 (0.435)	198 (0.439)	green(0.817)	0	3 (0.016)	723 (0.003)	green(0.995)	0
		5	1 (0)	723 (0.003)	green(1)	0.284	2 (0.004)	723 (0.003)	green(0.999)	0
0.25%	500	1	20 (0)	117 (0.871)	red(0)	0	20 (0)	117 (0.871)	red(0)	0
		2	2 (0.03)	480 (0.034)	green(0.991)	0.513	3 (0.102)	475 (0.036)	green(0.965)	0.764
		3	23 (0)	117 (0.871)	red(0)	0	24 (0)	117 (0.871)	red(0)	0
		4	7 (0.927)	473 (0.036)	green(0.515)	0	8 (0.641)	470 (0.037)	green(0.365)	0
		5	5 (0.476)	473 (0.036)	green(0.805)	0	4 (0.248)	475 (0.036)	green(0.906)	0
0.5%	250	1	37 (0)	14 (0.132)	red(0)	0	34 (0)	14 (0.132)	red(0)	0
		2	2 (0.004)	732 (0.003)	green(0.999)	0.352	1 (0)	891 (0.001)	green(1)	0.272
		3	39 (0)	14 (0.132)	red(0)	0	35 (0)	14 (0.132)	red(0)	0
		4	14 (0.146)	195 (0.451)	yellow(0.087)	0	12 (0.387)	198 (0.439)	green(0.222)	0
		5	7 (0.435)	596 (0.011)	green(0.817)	0	4 (0.050)	596 (0.011)	green(0.983)	0.594
0.5%	500	1	27 (0)	117 (0.871)	red(0)	0	37 (0)	117 (0.871)	red(0)	0
		2	4 (0.248)	480 (0.034)	green(0.906)	0	4 (0.248)	475 (0.034)	green(0.906)	0
		3	29 (0.022)	117 (0.871)	red(0)	0	33 (0)	66 (0.696)	red(0)	0
		4	17 (0.001)	173 (0.545)	yellow(0.001)	0	15 (0.006)	173 (0.545)	yellow(0.004)	0
		5	5 (0.476)	480 (0.034)	green(0.805)	0	7 (0.927)	346 (0.117)	green(0.515)	0.003
0.75%	250	1	43 (0)	14 (0.132)	red(0)	0	44 (0)	14 (0.132)	red(0)	0
		2	7 (0.435)	596 (0.011)	green(0.817)	0	4 (0.05)	596 (0.011)	green(0.983)	0.498
		3	44 (0)	14 (0.132)	red(0)	0	38 (0)	14 (0.132)	red(0)	0
		4	23 (0.007)	14 (0.132)	red(0)	0	21 (0.001)	14 (0.132)	yellow(0.001)	0
		5	10 (0.809)	596 (0.011)	green(0.447)	0	6 (0.250)	596 (0.011)	green(0.900)	0.001
0.75%	500	1	37 (0)	66 (0.696)	red(0)	0	43 (0)	66 (0.696)	red(0)	0
		2	6 (0.765)	173 (0.545)	green(0.669)	0.002	6 (0.765)	345 (0.118)	green(0.669)	0.001
		3	31 (0.000)	66 (0.696)	red(0)	0	34 (0)	66 (0.696)	red(0)	0
		4	18 (0.00)	173 (0.545)	red(0)	0	21 (0)	173 (0.545)	red(0)	0
		5	10 (0.242)	173 (0.545)	green(0.145)	0	10 (0.242)	345 (0.118)	green(0.145)	0

Table 9. Table 9: Results from sVaR backtesting using the DCC models on the test data (from 1 May 2016 to 31 December 2020) calibrated using five different stress periods b = 1 , … , 5 𝑏 1 … 5 b=1,\ldots,5 from Table  7 , with T sv = 650 subscript 𝑇 sv 650 T_{\text{sv}}=650 used throughout. From left to right, the table reports: The number of sVaR failures (PoF) and the number of days until the first failure in sVaR (TFF), with the corresponding p 𝑝 p -values in brackets from the respective two-sided tests in ( 13 ); traffic light test (traffic); p 𝑝 p -values from the DQ test (DQ) in ( 14 ).

Scenario		5-day rebalance ($h_{\text{reb}}=5$)				20-day rebalance ($h_{\text{reb}}=20$)
${\mu }_{\text{target}}$	Period	PoF	TFF	traffic	DQ	PoF	TFF	traffic	DQ
0.25%	1	10 (0.065)	12 (0.110)	yellow(0.041)	0.029	9 (0.137)	12 (0.110)	yellow(0.085)	0.036
	2	3(0.281)	383 (0.084)	green(0.897)	0.813	3(0.281)	383 (0.084)	green(0.897)	0.855
	3	16 (0)	12 (0.110)	red(0)	0	17 (0)	12 (0.110)	red(0)	0
	4	12 (0.011)	12 (0.110)	yellow(0.008)	0.008	12 (0.011)	12 (0.110)	yellow(0.008)	0.008
	5	10 (0.065)	12 (0.110)	yellow(0.019)	0.016	12 (0.110)	12 (0.110)	yellow(0.008)	0.006
0.5%	1	9 (0.137)	12 (0.110)	yellow(0.085)	0.040	9 (0.137)	12 (0.110)	yellow(0.085)	0.042
	2	5 (0.909)	295 (0.186)	green(0.605)	0.439	3 (0.281)	383 (0.084)	green(0.897)	0.906
	3	22 (0)	12 (0.110)	red(0)	0	22 (0)	12 (0.110)	red(0)	0
	4	10 (0.065)	12 (0.110)	yellow(0.041)	0.030	11 (0.028)	12 (0.110)	yellow(0.019)	0.017
	5	8 (0.265)	12 (0.110)	yellow(0.161)	0.039	8 (0.265)	12 (0.110)	yellow(0.161)	0.046
0.75%	1	10 (0.065)	12 (0.110)	yellow(0.041)	0.040	9 (0.137)	12 (0.110)	yellow(0.085)	0.044
	2	6 (0.751)	295 (0.186)	green(0.430)	0.248	4 (0.564)	383 (0.084)	green(0.771)	0.890
	3	23 (0)	12 (0.110)	red(0)	0	23 (0)	12 (0.110)	red(0)	0
	4	12 (0.011)	12 (0.110)	yellow(0.008)	0.009	12 (0.011)	12 (0.11)	yellow(0.008)	0.01
	5	9 (0.137)	12 (0.110)	yellow(0.085)	0.044	8 (0.265)	12 (0.110)	green(0.161)	0.046

Table 10. Table 10: (M4) Power of the change point test at α = 0.05 𝛼 0.05 \alpha=0.05 and the accuracy of η ^ 1 subscript ^ 𝜂 1 \widehat{\eta}_{1} (% of | η ^ 1 − η 1 | < log 2 ⁡ T subscript ^ 𝜂 1 subscript 𝜂 1 superscript 2 𝑇 |\widehat{\eta}_{1}-\eta_{1}|<\log^{2}T ) when N = 50 𝑁 50 N=50 , T = 1000 𝑇 1000 T=1000 and η 1 = [ T / 2 ] subscript 𝜂 1 delimited-[] 𝑇 2 \eta_{1}=[T/2] .

	(M4)(M4.1)		(M4)(M4.2)		(M4)(M4.3)		(M4)(M4.4)
$\varrho$	power	accuracy (%)	power	accuracy (%)	power	accuracy (%)	power	accuracy (%)
1	1.00	100	1.00	100	1.00	100	0.96	84
0.75	1.00	100	1.00	100	1.00	100	0.91	79
0.5	1.00	100	1.00	99	1.00	100	0.57	47
0.25	0.96	91	0.83	73	0.99	97	0.12	7
	(M4)(M4.5)		(M4)(M4.6)		(M4)(M4.7)		(M4)(M4.8)
$\varrho$	power	accuracy (%)	power	accuracy (%)	power	accuracy (%)	power	accuracy (%)
1	1.00	100	1.00	100	1.00	100	1.00	98
0.75	1.00	100	1.00	100	1.00	100	1.00	98
0.5	1.00	100	1.00	100	1.00	100	0.98	93
0.25	0.99	93	1.00	98	1.00	99	0.54	43

Table 11. Table 11: (M4) Power of the change point test at α = 0.05 𝛼 0.05 \alpha=0.05 and the accuracy of η ^ 1 subscript ^ 𝜂 1 \widehat{\eta}_{1} (% of | η ^ 1 − η 1 | < log 2 ⁡ T subscript ^ 𝜂 1 subscript 𝜂 1 superscript 2 𝑇 |\widehat{\eta}_{1}-\eta_{1}|<\log^{2}T ) when N = 50 𝑁 50 N=50 , T = 1000 𝑇 1000 T=1000 and η 1 = [ 0.9 ​ T ] subscript 𝜂 1 delimited-[] 0.9 𝑇 \eta_{1}=[0.9T] .

	(M4)(M4.1)		(M4)(M4.2)		(M4)(M4.3)		(M4)(M4.4)
$\varrho$	power	accuracy (%)	power	accuracy (%)	power	accuracy (%)	power	accuracy (%)
1	1.00	100	1.00	99	1.00	100	0.43	37
0.75	1.00	100	1.00	99	1.00	100	0.27	20
0.5	1.00	100	0.89	84	1.00	100	0.09	6
0.25	0.92	90	0.32	28	0.96	95	0.08	2
	(M4)(M4.5)		(M4)(M4.6)		(M4)(M4.7)		(M4)(M4.8)
$\varrho$	power	accuracy (%)	power	accuracy (%)	power	accuracy (%)	power	accuracy (%)
1	1.00	100	1.00	100	1.00	100	0.96	84
0.75	1.00	100	0.99	97	1.00	100	0.89	79
0.5	1.00	100	0.87	86	1.00	100	0.61	50
0.25	0.63	62	0.10	8	1.00	100	0.16	10

Table 12. Table 12: (M4) Power of the change point test at α = 0.05 𝛼 0.05 \alpha=0.05 and the accuracy of η ^ 1 subscript ^ 𝜂 1 \widehat{\eta}_{1} when N = 100 𝑁 100 N=100 , T = 1000 𝑇 1000 T=1000 and η 1 = [ T / 2 ] subscript 𝜂 1 delimited-[] 𝑇 2 \eta_{1}=[T/2] .

	(M4)(M4.1)		(M4)(M4.2)		(M4)(M4.3)		(M4)(M4.4)
$\varrho$	power	accuracy (%)	power	accuracy (%)	power	accuracy (%)	power	accuracy (%)
1	1.00	100	1.00	100	1.00	100	1.00	96
0.75	1.00	100	1.00	100	1.00	100	1.00	93
0.5	1.00	100	1.00	100	1.00	100	0.94	83
0.25	1.00	98	1.00	90	1.00	100	0.46	36
	(M4)(M4.5)		(M4)(M4.6)		(M4)(M4.6)		(M4)(M4.8)
$\varrho$	power	accuracy (%)	power	accuracy (%)	power	accuracy (%)	power	accuracy (%)
1	1.00	100	1.00	100	1.00	100	1.00	100
0.75	1.00	100	1.00	100	1.00	100	1.00	100
0.5	1.00	100	1.00	100	1.00	100	1.00	99
0.25	1.00	100	1.00	100	1.00	100	0.98	93

Table 13. Table 13: (M4) Power of the change point test at α = 0.05 𝛼 0.05 \alpha=0.05 and the accuracy of η ^ 1 subscript ^ 𝜂 1 \widehat{\eta}_{1} when N = 100 𝑁 100 N=100 , T = 1000 𝑇 1000 T=1000 and η 1 = [ 9 ​ T / 10 ] subscript 𝜂 1 delimited-[] 9 𝑇 10 \eta_{1}=[9T/10] .

	(M4)(M4.1)		(M4)(M4.2)		(M4)(M4.3)		(M4)(M4.4)
$\varrho$	power	accuracy (%)	power	accuracy (%)	power	accuracy (%)	power	accuracy (%)
1	1.00	100	1.00	100	1.00	100	0.87	78
0.75	1.00	100	1.00	100	1.00	100	0.65	57
0.5	1.00	100	1.00	100	1.00	100	0.18	15
0.25	1.00	100	0.81	76	1.00	100	0.07	1
	(M4)(M4.5)		(M4)(M4.6)		(M4)(M4.7)		(M4)(M4.8)
$\varrho$	power	accuracy (%)	power	accuracy (%)	power	accuracy (%)	power	accuracy (%)
1	1.00	100	1.00	100	1.00	100	1.00	98
0.75	1.00	100	1.00	100	1.00	100	1.00	97
0.5	1.00	100	1.00	100	1.00	100	0.93	84
0.25	1.00	100	0.52	52	1.00	100	0.31	26

Equations

$$r_{i,t}$$

$${\mathbf{\Omega }}_i(t)=({\omega }_i(t),{\alpha }_{i,1}(t),\dots ,{\alpha }_{i,p}(t),{\beta }_{i,1}(t),\dots ,{\beta }_{i,q}(t){)}^{\top }\in {\mathbb{R}}^{1+p+q},$$

$${\mathbf{\Theta }}_i(t)=({\sigma }_{i,1}(t),\dots ,{\sigma }_{i,i-1}(t),{\sigma }_{i,i+1}(t),\dots ,{\sigma }_{i,N}(t){)}^{\top }\in {\mathbb{R}}^{N-1}.$$

$$\underset{1\le i\le N}{\min }\underset{t\in \mathbb{Z}}{\inf }{\omega }_i(t)>{ϵ}_1\text{ and }\underset{1\le i\le N}{\max }\underset{t\in \mathbb{Z}}{\sup }{\omega }_i(t)\le {\Xi }_1<\infty .$$

$$\underset{1\le i\le N}{\max }\underset{t\in \mathbb{Z}}{\sup }\left(\sum _{j=1}^p{\alpha }_{i,j}(t)+\sum _{k=1}^q{\beta }_{i,k}(t)\right)\le 1-{ϵ}_2.$$

$$\left(\int |f_{i,i^{\prime }}(u,v)-f_{i,i^{\prime }}(u(1+a),v)|dudv\right)\vee \left(\int |f_{i,i^{\prime }}(u,v)-f_{i,i^{\prime }}(u,v(1+a))|dudv\right)\le Ka$$

$$\underset{1\le i<i^{\prime }\le N}{\sup }\underset{\begin{array}{c}G\in \sigma ({\mathbf{r}}_{i,i^{\prime },u}:u\ge t+k),\\ H\in \sigma ({\mathbf{r}}_{i,i^{\prime },u}:u\le t)\end{array}}{\sup }|\mathbb{P}(G\cap H)-\mathbb{P}(G)\mathbb{P}(H)|\le M{\alpha }^k,$$

$$U_{i,t}$$

$${\check{h}}_{i,t}$$

$$U_{i,t}^2=\frac{h_{i,t}}{{\check{h}}_{i,t}}\cdot \frac{r_{i,t}^2}{h_{i,t}}=\frac{h_{i,t}}{{\check{h}}_{i,t}}{\epsilon }_{i,t}^2=\frac{h_{i,t}}{{\check{h}}_{i,t}}+\frac{h_{i,t}}{{\check{h}}_{i,t}}({\epsilon }_{i,t}^2-1),$$

$$g_1({\mathbf{r}}_{i,t}^{t-p},{\mathbf{h}}_{i,t-1}^{t-q})\equiv {\left\{g_0({\mathbf{r}}_{i,t}^{t-p},{\mathbf{h}}_{i,t-1}^{t-q})\right\}}^2=U_{i,t}^2=\frac{r_{i,t}^2}{{\check{h}}_{i,t}}.$$

$$U_{i{i}^{\prime },t}=g_2({\mathbf{r}}_{i,t}^{t-p},{\mathbf{h}}_{i,t-1}^{t-q},{\mathbf{r}}_{i^{\prime },t}^{t-p},{\mathbf{h}}_{i^{\prime },t-1}^{t-q})=(U_{i,t}+s_{i,i^{\prime }}{U}_{i^{\prime },t}{)}^2$$

$$U_{i{i}^{\prime },t}=$$

$$+2s_{i,i^{\prime }}\sqrt{\frac{h_{i,t}{h}_{i^{\prime },t}}{{\check{h}}_{i,t}{\check{h}}_{i^{\prime },t}}}\left\{{\sigma }_{i,i^{\prime }}(t)+({\epsilon }_{i,t}{\epsilon }_{i^{\prime },t}-{\sigma }_{i,i^{\prime }}(t))\right\},$$

$$\{U_{i,t}^2,1\le i\le N,U_{i{i}^{\prime },t},1\le i<i^{\prime }\le N;1\le t\le T\}.$$

$$x_{j,t}=U_{i{i}^{\prime },t}\text{with}j\equiv j(i,i^{\prime })=(N-i/2)(i-1)+i^{\prime },1\le i\le i^{\prime }\le N,$$

$${\mathcal{X}}_{s,c,e}^j=\sqrt{\frac{(c-s+1)(e-c)}{e-s+1}}\left(\frac{1}{c-s+1}\sum _{t=s}^c{x}_{j,t}-\frac{1}{e-c}\sum _{t=c+1}^e{x}_{j,t}\right)\text{for}s\le c<e.$$

$${\mathcal{D}}_{s,e}(c,m)=\sqrt{\frac{m(2d-m)}{2d}}\left(\frac{1}{m}\sum _{j=1}^m|{\mathcal{X}}_{s,c,e}^{(j)}|-\frac{1}{2d-m}\sum _{j=n+1}^d|{\mathcal{X}}_{s,c,e}^{(j)}|\right)$$

$${\mathcal{T}}_{s,e}=\underset{s\le c<e}{\max }\underset{1\le m\le d}{\max }{\mathcal{D}}_{s,e}(c,m),$$

$$\hat{\eta }=\arg \underset{s\le c<e}{\max }\underset{1\le m\le d}{\max }{\mathcal{D}}_{s,e}(c,m).$$

$$x_{j,t}=f_{j,t}+z_{j,t},1\le j\le d,1\le t\le T.$$

$$\displaystyle\widetilde{\Pi}_{b}=\Big{\{}j:\,\Delta_{j,b}:=|f_{j,\eta_{b}+1}-f_{j,\eta_{b}}|\neq 0\Big{\}}.$$

$$\underset{1\le j\le d}{\max }\underset{1\le s<e\le T}{\max }\frac{1}{\sqrt{e-s+1}}\left|\sum _{t=s}^e{z}_{j,t}\right|=O_p(\sqrt{\log (T)}).$$

$$\frac{T^{7\gamma /4-3/2}{\underline{\Delta }}_{d,T}}{\sqrt{d\log (T)}}\to \infty \text{ as }T\to \infty ,\text{where}{\underline{\Delta }}_{d,T}:=\underset{1\le b\le B}{\min }|{\tilde{\Pi }}_b{|}^{-1/2}\sum _{j\in {\tilde{\Pi }}_b}{\Delta }_{j,b}.$$

$$\mathbb{P}\left\{\hat{B}=B;|{\hat{\eta }}_b-{\eta }_b|<c_0{\rho }_{d,T}\text{ for }b=1,\dots ,B\right\}\to 1$$

$$\displaystyle F_{i}=\max\left[1,\frac{\min(0.99,\sum_{j=1}^{p}\widehat{\alpha}_{i,j}+\sum_{k=1}^{q}\widehat{\beta}_{i,k})}{\max\big{\{}0.01,1-(\sum_{j=1}^{p}\widehat{\alpha}_{i,j}+\sum_{k=1}^{q}\widehat{\beta}_{i,k})\big{\}}}\right].$$

$${\hat{h}}_{i,t}={\hat{\omega }}_i+\sum _{j=1}^p{\hat{\alpha }}_{i,j}{r}_{i,t-j}^2+\sum _{k=1}^q{\hat{\beta }}_{i,k}{\hat{h}}_{i,t-k}$$

$$r_{i,t}^{\ell }=(h_{i,t}^{\ell }{)}^{1/2}{\epsilon }_{i,t}^{\ell },\text{where}{h}_{i,t}^{\ell }={\hat{\omega }}_i+\sum _{j=1}^p{\hat{\alpha }}_{i,j}(r_{i,t-j}^{\ell }{)}^2+\sum _{k=1}^q{\hat{\beta }}_{i,k}{h}_{i,t-k}^{\ell }$$

$$\displaystyle\Big{\{}g_{1}(\mathbf{r}^{\ell,t-p}_{i,t},\mathbf{h}^{\ell,t-q}_{i,t-1}),\,1\leq i\leq N,\,g_{2}(\mathbf{r}^{\ell,t-p}_{i,t},\mathbf{h}^{\ell,t-q}_{i,t-1},\mathbf{r}^{\ell,t-p}_{i^{\prime},t},\mathbf{h}^{\ell,t-q}_{i^{\prime},t-1}),\,1\leq i<i^{\prime}\leq N;\,1\leq t\leq T\Big{\}}$$

$${\mathbf{r}}_t$$

$$\text{where }{h}_{i,t}$$

$${\text{VaR}}_t(\alpha )=-F^{-1}(\alpha |{\mathcal{G}}_t),$$

Full text

Keywords: Value-at-Risk; stress period selection; data segmentation; multivariate GARCH; high dimensionality

11footnotetext: School of Mathematics, University of Bristol, UK. Email: haeran.cho@bristol.ac.uk

Supported by the Engineering and Physical Sciences Research Council grant no. EP/N024435/1.22footnotetext: School of Computing and Mathematical Sciences, University of Greenwich, UK. Email: kkorkas@yahoo.co.uk

1 Introduction

The increased financial uncertainty during the recent global economic crisis has confirmed the close volatility linkage between asset markets. For example, there is strong evidence that the economy and oil prices (Hamilton, , 2003), foreign exchange rates (Baillie, , 1991), equity markets (Baele, , 2003) and crude oil and agricultural commodities (Du et al., , 2011) are related. These co-movements are naturally expected since the rate of information influences the volatility in asset returns and, therefore, the information flow from one asset market can be incorporated into another related market (Ross, , 1989). Hence, good understanding of the correlations among multiple markets is crucial for policy makers, financial institutions and investors.

The joint modelling of financial returns as a multivariate generalised autoregressive conditional heteroscedastic (GARCH) process has attracted considerable attention in the literature. A partial list includes the vectorised multivariate GARCH model (Bollerslev et al., , 1988), Baba-Engle-Kraft-Kroner (BEKK) model (Engle and Kroner, , 1995), constant conditional correlation (CCC) model (Bollerslev, , 1990), dynamic conditional correlation (DCC) model (Engle, , 2002), generalised orthogonal GARCH model (Van der Weide, , 2002), full-factor multivariate GARCH model (Vrontos et al., , 2003) and conditionally uncorrelated components-based multivariate volatility processes (Fan et al., , 2008); for a survey of multivariate GARCH modelling and inference, see Bauwens et al., (2006).

The assumption that the underlying dynamics remain unchanged is restrictive considering that the fundamentals driving an economy, the asset markets in particular, exhibit sudden changes or regimes switches. Empirical evidence of change points (a.k.a. structural breaks or breakpoints) in various macroeconomic and financial time series is well documented, such as those in exchange rates (Pesaran and Timmermann, , 2007), equity and bond returns (Cappiello et al., , 2006), commodity (Ewing and Malik, , 2013) and share (Barigozzi et al., , 2018) prices. Further, Diebold and Inoue, (2001) and Mikosch and Stărică, (2004) noted that stochastic regime switching may be confused with long-range dependence. In this paper, we observe the importance of accounting for structural breaks in the volatilities and correlations of a multi-asset portfolio. In particular, we show that the stressed Value-at-Risk, a popular measure of market risk widely adopted by financial institutions, under-estimates the exposure of a portfolio when a pre-selected period of a fixed length is used as the stress period, compared to when the most volatile period is identified by change point analysis.

The problem of testing for at most a single structural break in univariate conditional heteroscedastic models has been studied in Kokoszka and Leipus, (2000), Kokoszka and Teyssière, (2002), Lee et al., (2003), Berkes et al., (2004) and De Pooter and Van Dijk, (2004), to name a few. For multiple change point detection, Fryzlewicz and Subba Rao, (2014) proposed the BASTA (binary segmentation for transformed ARCH), a two-stage procedure for detecting change points in the conditional variance of univariate series, while Andreou and Ghysels, (2003) studied change point detection in the co-movement of bivariate returns. More recent change point methods for (conditional) covariance structure of multivariate data include Aue et al., (2009), Dette et al., (2018) and Barassi et al., (2020) on testing for a single change point and Wang et al., (2021) on multiple change point detection. We also note some recent works on homogeneity pursuit for panel data with factor structures (Ke et al., , 2016), where the focus lies in identifying a group structure in the panel data by casting the problem as that of change point detection.

In this paper, we propose a methodology for multiple change point detection in multivariate, possibly high-dimensional GARCH processes. It simultaneously segments high-dimensional GARCH processes by identifying ‘common’ change points, each of which can be shared by a subset or all of the component time series as a change point in their within-series and/or cross-sectional dependence structure. The methodology consists of two stages: The first stage transforms the $N$-dimensional time series into the $N(N+1)/2$-dimensional panel data consisting of empirical residual series and their cross-products, whereby change points in the complex structure (conditional variance and covariance) are made detectable as change points in the level of the transformed data, at the price of the increased dimensionality. A number of methodologies have been investigated for change point analysis in the mean of high-dimensional panel data, such as Horváth and Hušková, (2012), Jirak, (2015), Cho and Fryzlewicz, (2015) and Wang and Samworth, (2018). Among many, we adopt the double CUSUM Binary Segmentation procedure proposed in Cho, (2016) for the second stage, which achieves consistency in multiple change point estimation while permitting within-series and cross-sectional correlations. Extending the mixing property originally derived for univariate, time-varying ARCH processes in Fryzlewicz and Subba Rao, (2011) to that of time-varying bivariate GARCH processes (Proposition 1), we establish the consistency of the combined methodology (Theorem 1).

The rest of the paper is organised as follows. Section 2 introduces a time-varying multivariate GARCH model which provides a framework for the theoretical treatment of our methodology. Section 3 describes the proposed two-stage methodology and investigates its theoretical properties, and simulation studies in Section 4 confirm its good finite sample performance. In Section 5, we apply the methodology to a real financial data set and show the importance of accounting for structural breaks in a portfolio of assets in risk management. Section 6 concludes the paper, and the Appendix contains the proofs of theoretical results and additional simulation results. Our methodology is implemented in the R package segMGarch, available from CRAN.

Notation

By $\mathbb{R}_{+}$, we denote the set of positive real numbers. For any set $\Pi\subset\{1,\ldots,N\}$, we denote its cardinality by $|\Pi|$. For given observations $\mathbf{r}_{t}\in\mathbb{R}^{N}$, we denote by $\mathcal{F}_{t}$ the $\sigma$-algebra $\sigma\{\mathbf{r}_{s},\,s\leq t\}$. Also, we use the notations $a\vee b=\max(a,b)$ and $a\wedge b=\min(a,b)$. The notation $a_{n}\asymp b_{n}$ indicates that $a_{n}$ is of the order of $b_{n}$, i.e. $a_{n}=O(b_{n})$ and $b_{n}=O(a_{n})$, and $a_{n}\gg b_{n}$ indicates that $a_{n}^{-1}b_{n}\to 0$. We denote a vector of zeros by $\mathbf{0}$ whose dimension should be clear from the context.

2 Time-varying multivariate GARCH model

We consider the following time-varying multivariate GARCH (tv-MGARCH) model denoted by $\mathbf{r}_{t}=(r_{1,t},\ldots,r_{N,t})^{\top}$, $t=1,\ldots,T$:

[TABLE]

The independent innovations $\bm{\varepsilon}_{t}=(\varepsilon_{1,t},\ldots,\varepsilon_{N,t})^{\top}$ satisfy $\mathbb{E}(\bm{\varepsilon}_{t})=\mathbf{0}_{N}$ and $\text{var}(\bm{\varepsilon}_{t})=\bm{\Sigma}_{\varepsilon}(t)=[\sigma_{i,i^{\prime}}(t)]_{i,i^{\prime}=1}^{N}$ with $\sigma_{i,i}(t)=1$ and $\sigma_{i,i^{\prime}}(t)=\sigma_{i^{\prime},i}(t)$ for all $1\leq i,i^{\prime}\leq N$ and $t$. We denote the vector of parameters involved in modelling the within-panel conditional variance of $r_{i,t}$ by

[TABLE]

and that involved in modelling the cross-correlations of $\varepsilon_{i,t}$ with $\varepsilon_{i^{\prime},t},\,i^{\prime}\neq i$, by

[TABLE]

Then, the following assumption specifies how change points enter into the structure of $\mathbf{r}_{t}$.

1. (A1)

The functions $\bm{\Omega}_{i}(t)$ and $\bm{\Theta}_{i}(t)$ are piecewise constant in $t$ and share $B\equiv B_{T}$ change points $\eta_{b}\equiv\eta_{b,T},\,b=1,\ldots,B$ (satisfying $0\equiv\eta_{0}<\eta_{1}<\ldots<\eta_{B}<\eta_{B+1}\equiv T$) across $i=1,\ldots,N$, such that at any $\eta_{b}$, there exists $\Pi_{b}\subset\{(i,i^{\prime}):\,1\leq i\leq i^{\prime}\leq N\}$ with $|\Pi_{b}|\geq 1$, where $(i,i^{\prime})\in\Pi_{b}$ iff either $\bm{\Omega}_{i}(\eta_{b})\neq\bm{\Omega}_{i}(\eta_{b}+1)$ (then $(i,i)\in\Pi_{b}$) or $\sigma_{i,i^{\prime}}(\eta_{b})\neq\sigma_{i,i^{\prime}}(\eta_{b}+1)$.

Under Assumption (A1), the tv-MGARCH process $\mathbf{r}_{t}$ is approximately piecewise stationary over each segment $[\eta_{b}+1,\eta_{b+1}],\,b=0,\ldots,B$, with boundary effects around the change points. Assumption (A1) does not rule out that $|\Pi_{b}|<N(N+1)/2$, i.e. it is permitted that only a subset of the components of $\bm{\Omega}_{i}(t)$ and $\bm{\Theta}_{i}(t)$ undergoes a change at each $\eta_{b}$.

Under model (1), it is assumed that the (unconditional) correlations across the components of $\mathbf{r}_{t}$ are attributed to those of $\bm{\varepsilon}_{t}$, and that the conditional variance of each component series is separately modelled as a univariate GARCH process, and the model is reduced to the CCC model of Bollerslev, (1990) over each stationary segment $[\eta_{b}+1,\eta_{b+1}]$. We adopt this simplistic approach in order to specify the parameters ($\bm{\Omega}_{i}(t)$ and $\bm{\Theta}_{i}(t)$) that introduce structural changes and to motivate the data transformation introduced later for change point analysis in high-dimensional panel data, rather than to provide a complete model for time-varying cross-sectional conditional correlations of the panel data; the latter task can be accomplished once stationary segments are identified by the estimated change points. We note that the literature on multivariate GARCH processes, as those cited in Introduction, considers a relatively lower-dimensional applications ($N\leq 8$), whereas we consider GARCH modelling of both simulated and real datasets of higher dimensions ($N$ up to $100$) in this paper, which renders the simplistic approach essential for computational feasibility.

We further assume the following assumption on the model (1).

1. (A2)

The dimensionality $N$ satisfies $N\asymp T^{\theta}$ for some $\theta\in[0,\infty)$. 2. (A3)

For some $\epsilon_{1}>0$, $\Xi_{1}<\infty$ and all $T$, we have

[TABLE] 3. (A4)

For some $\epsilon_{2}\in(0,1)$ and all $T$, we have

[TABLE]

Assumption (A2) indicates that the dimensionality can either be fixed or increase with $T$ at a polynomial rate. Assumptions (A3)–(A4) guarantee that between any two consecutive change points, each $r_{i,t}$ admits a well-defined solution a.s. and is weakly stationary (see e.g. Theorem 4.35 of Douc et al., (2014)).

Boussama et al., (2011) shows that there exists a unique, strictly stationary and strongly mixing solution for the BEKK multivariate GARCH models introduced in Engle and Kroner, (1995), see their Theorem 2.4. Relatively few results exist on mixing for non-stationary processes, but such results are essential in investigating the properties of the change point detection methodology developed under (1). Adopting the tools developed in Fryzlewicz and Subba Rao, (2011), who investigate the mixing rate of univariate, time-varying ARCH processes, we derive that any pair of the components of a tv-MGARCH process, $\mathbf{r}_{ii^{\prime},t}=(r_{i,t},r_{i^{\prime},t})^{\top}$, $1\leq i<i^{\prime}\leq N$, is strong mixing at a geometric rate under the following Lipschitz-type condition on the joint density of $(\varepsilon_{i,t}^{2},\varepsilon_{i^{\prime},t}^{2})^{\top}$.

1. (A5)

The joint distribution of $\varepsilon_{i,t}^{2}$ and $\varepsilon_{i^{\prime},t}^{2}$, denoted by $f_{i,i^{\prime}}(u,v)$, satisfies the following: for any $a>0$, there exists fixed $K>0$ independent of $a$ such that

[TABLE]

uniformly over $i,i^{\prime}=1,\ldots,N,\,i\neq i^{\prime}$.

Proposition 1**.**

Under Assumptions (A1) and (A3)–(A5), there exists some $\alpha\in(0,1)$ such that

[TABLE]

where $M$ is a finite constant independent of $t$ and $k$.

See Appendix A.2 for the proof, which relies on the Markovian nature of the time-varying bivariate GARCH processes. The strong mixing property of the tv-MGARCH process derived in Proposition 1 proves useful in controlling the behaviour of empirical residuals from the data transformation proposed in Section 3.1 below, see Proposition 2 (ii), and this in turn allows us to rigorously establish the consistency of the proposed change point detection methodology.

3 Two-stage change point detection methodology

In this section, we describe the proposed two-stage change point detection methodology under the tv-MGARCH model in (1). More specifically, Section 3.1 introduces a data transformation of the $N$-dimensional tv-MGARCH process into a panel of $d\equiv d_{N}=N(N+1)/2$ series, which makes any change in its conditional variance and correlation structure as defined in (A1), detectable as that in the mean of the transformed data. Then in Section 3.2, we outline the application of the double CUSUM binary segmentation algorithm (Cho, , 2016) to the transformed panel data for the detection and localisation of the multiple change points. Section 3.3 establishes the consistency of the combined two-stage procedure in estimating both the total number and the locations of the change points $\eta_{b},\,b=1,\ldots,B$, and Section 3.4 discusses the selection of any tuning parameters involved in the methodology.

3.1 Stage 1: Transformation of tv-MGARCH processes

We first describe a data transformation for detecting change points in $\bm{\Omega}_{i}(t)$ for each $i=1,\ldots,N$, which involves a function $g_{0}:\mathbb{R}^{1+p+q}\to\mathbb{R}$ that takes $\mathbf{r}_{i,t}^{t-p}=(r_{i,t},\ldots,r_{i,t-p})^{\top}$ and $\mathbf{h}_{i,t-1}^{t-q}=(h_{i,t-1},\ldots,h_{i,t-q})^{\top}$ as an input, and generates:

[TABLE]

where $C_{i,0}\in\mathbb{R}_{+}$, $C_{i,j}\in\mathbb{R}_{+}\cup\{0\},\,j=1,\ldots,p+q$, and $\epsilon$ denotes a small positive constant. For now, we ignore that the conditional variance $\mathbf{h}_{i,t-1}^{t-q}$ is unobservable and that the GARCH orders $p$ and $q$ are unknown; we address their selection and the choice of $C_{i,j}$ in Section 3.4.

While it is possible to adopt an alternative transformation as $g_{0}$, the proposed transformation in (2) is closely related to empirical residuals which have been popularly adopted for change point analysis in time-varying and conditionally heteroscedastic univariate processes, see Kokoszka and Teyssière, (2002), Lee et al., (2003), De Pooter and Van Dijk, (2004), Fryzlewicz and Subba Rao, (2014) and Barassi et al., (2020). For stationary processes, empirical residuals approximates a series of i.i.d. innovations and, even in the presence of change points, such transformation tends to reduce the autocorrelations. Moreover, $U_{i,t}$ ‘encodes’ the presence and the locations of any change point in $\bm{\Omega}_{i}(t)$. More specifically,

[TABLE]

so that $U^{2}_{i,t}$ contains any change in $\bm{\Omega}_{i}(t)$ as a change in its ‘level’, recalling that $\mathbb{E}(\varepsilon_{i,t}^{2})=1$. For notational convenience, we define

[TABLE]

To detect any changes in the cross-sectional dependence structure of $\mathbf{r}_{t}$, we adopt the transformation $g_{2}:\mathbb{R}^{2+2p+2q}\to\mathbb{R}$:

[TABLE]

where $s_{i,i^{\prime}}\in\{1,-1\}$. Similarly to (3), $U_{ii^{\prime},t}$ admits the decomposition

[TABLE]

from which we conclude that a change in $\sigma_{i,i^{\prime}}(t)$ are detectable from $U_{ii^{\prime},t}$ as that in its level.

Regarding $U_{i,t}$ as empirical residuals obtained by applying volatility filters, Andreou and Ghysels, (2003) propose to examine $U_{1,t}U_{2,t}$, $U_{1,t}^{2}U_{2,t}^{2}$ or $|U_{1,t}U_{2,t}|$ for detecting structural changes in the co-movement of a pair of series $(r_{1,t},r_{2,t})^{\top}$. Instead, we adopt $U_{ii^{\prime},t}$, whose formulation is motivated by the observation made in Cho and Fryzlewicz, (2015): For given $(a_{t},b_{t})$, any changes in the second-order dependence structure, $\mathbb{E}(a_{t}^{2})$, $\mathbb{E}(b_{t}^{2})$ and $\mathbb{E}(a_{t}b_{t})$, are detectable by jointly examining $\mathbb{E}(a_{t}^{2})$, $\mathbb{E}(b_{t}^{2})$ and $\mathbb{E}\{(a_{t}+s_{a,b}b_{t})^{2}\}$ for any $s_{a,b}\in\{1,-1\}$. Moreover, $U_{i,t}^{2}$, $U_{i^{\prime},t}^{2}$ and $U_{ii^{\prime},t}$ are in the same scale in the sense that if $U_{i,t}$ and $U_{i^{\prime},t}$ are Gaussian random variables, $U_{i,t}^{2}$, $U_{i^{\prime},t}^{2}$ and $U_{ii^{\prime},t}$ are distributed as random variables following scaled $\chi^{2}_{1}$ distributions. This enables us to regard $U_{ii^{\prime},t}$ on an equal footing with $U_{i,t}^{2}$ and $U_{i^{\prime},t}^{2}$, which is essential in simultaneous segmentation of $\mathbf{r}_{t}$ via joint consideration of all the coordinate series of the panel data

[TABLE]

Remark 1* (Choice of $s_{i,i^{\prime}}$).*

Our theoretical results do not depend on the choice of $s_{i,i^{\prime}}$. In practice, we set $s_{i,i^{\prime}}=-\text{sign}(\widehat{\text{cor}}(U_{i,t},U_{i^{\prime},t}))$ where $\widehat{\text{cor}}(U_{i,t},U_{i^{\prime},t})$ denotes the sample correlation between $U_{i,t}$ and $U_{i^{\prime},t}$ over $t=1,\ldots,T$, in an attempt to better bring out any change in $\sigma_{i,i^{\prime}}(t)$ as that in the level of $U_{ii^{\prime},t}$; a similar choice has also been considered in Cho and Fryzlewicz, (2015) and Barigozzi et al., (2018).

In summary, we transform the $N$-dimensional time series $\mathbf{r}_{t}$ to the $d$-dimensional panel data with $d=N(N+1)/2$ in (5), where any changes in $\bm{\Omega}_{i}(t)$ and $\bm{\Theta}_{i}(t)$, $i=1,\ldots,N$ are detectable as changes in the level of the coordinate time series. In doing so, the problem of detecting multiple change points in the more complex dependence structure of $\mathbf{r}_{t}$, is transformed into a relatively easier problem of detecting change points in the level of the panel data in (5). While this brings in the increase of dimensionality, the panel data segmentation algorithm adopted in the second stage (detailed in Section 3.2) handles the high dimensionality well under a mild assumption on the dimensionality in Assumption (A2).

3.2 Stage 2: Double CUSUM binary segmentation

The double CUSUM binary segmentation (DCBS) algorithm. first introduced in Cho, (2016), is applicable to high-dimensional panel data for the detection and estimation of multiple change points in its level. The $d$-dimensional panel in (5) obtained from transforming the original time series $\mathbf{r}_{t}$, encodes the change points under (1) in its level as described in Section 3.1, and thus can serve as an input to the DCBS algorithm. For ease of notation, we introduce the DCBS algorithm with each element time series of the panel data in (5) denoted by

[TABLE]

(where $U_{ii,t}\equiv U_{i,t}^{2}$) such that $j\in\{1,\ldots,d\}$.

We first describe the double cumulative sum (CUSUM) statistics computed on a generic segment $[s,e]$ for some $1\leq s<e\leq T$, and then provide a full description of the DCBS algorithm. CUSUM statistics have been widely adopted for change point detection in both univariate and multivariate data. A series of (weighted) CUSUM statistics, calculated on $x_{j,t},\,s\leq t\leq e$, is given by

[TABLE]

In the context of testing the null hypothesis of no change point against the at-most-one-change alternative, asymptotic properties of the CUSUM-based test have been extensively studied, see e.g. Csörgö and Horváth, (1997) and Aue and Horváth, (2013). A large value of the CUSUM statistic $|\mathcal{X}^{j}_{s,c,e}|$ indicates the presence of a change point in the level of $x_{j,t}$ in the vicinity of $t=c$. Combined with the binary segmentation algorithm, the CUSUM statistic has frequently been adopted for univariate time series segmentation; see Vostrikova, (1981) and Venkatraman, (1992) for the theoretical treatment of its application to the canonical additive model with independent noise, and Fryzlewicz and Subba Rao, (2014) the case of piecewise stationary ARCH processes.

Cho, (2016) propose the DCBS algorithm for simultaneous segmentation of multivariate, possibly high-dimensional time series panel data. It guarantees consistency in estimating both the total number and locations of multiple change points, while permitting both serial and cross-sectional correlations in $x_{j,t}$, which is highly relevant to the time series setting considered in this paper. Over a given segment $[s,e]$, we aggregate the CUSUM statistics $\mathcal{X}^{j}_{s,c,e}$ over $1\leq j\leq d$, and generates a two-dimensional array of double CUSUM (DC) statistics

[TABLE]

for $s\leq c<e$ and $1\leq m\leq d$, where $|\mathcal{X}^{(j)}_{s,c,e}|$ denote the ordered CUSUM statistics at each $c$ such that $|\mathcal{X}^{(1)}_{s,c,e}|\geq|\mathcal{X}^{(2)}_{s,c,e}|\geq\ldots\geq|\mathcal{X}^{(d)}_{s,c,e}|$. Then, the test statistic is derived by maximising the two-dimensional array over both time and cross-sectional indices, as

[TABLE]

which is compared against a threshold, $\pi_{d,T}$, for determining the presence of a change point over the interval $[s,e]$. If $\mathcal{T}_{s,e}>\pi_{d,T}$, we regard that there exists at least one change point in $[s,e]$ and identify the location of a change point as

[TABLE]

Remark 2*.*

We briefly remark upon the choice of scaling for the DC statistics in (8), namely $\sqrt{m(2d-m)/(2d)}$ in place of $\sqrt{m(d-m)/d}$. While the latter is in line with the scaling adopted for the weighted CUSUM statistic in (7), when applied to the cross-sectional dimension, it does not favour a change point shared by more than $d/2$ rows of the panel data and in fact, acts as a penalty when all $d$ rows share a change point, which is counter-intuitive. The choice of the former scaling resolves this issue, and favourably regards the ‘density’ of a change point in considering its cross-sectional magnitude.

The DCBS algorithm recursively applies the above steps of testing and locating a single change point over the segments determined by the previously detected change points, until they are no longer partitioned according to a threshold $\pi_{d,T}$. A pseudocode of the DCBS algorithm is provided in Algorithm 1. The call of DCBinSeg with $s=1$, $e=T$ and $\widehat{\mathcal{B}}=\emptyset$ as the initial arguments returns the set of change point estimators $\widehat{\mathcal{B}}=\{\widehat{\eta}_{b}:\,b=1,\ldots,\widehat{B}\}$. We discuss the choice of $\pi_{d,T}$ in Section 3.4.

3.3 Theoretical properties

Theoretical investigation into the properties of the two-stage procedure is divided into two parts; firstly, we show that the transformed panel data in (5) admits a decomposition into piecewise constant signals and the residuals meeting some desirable properties. Then, we establish the consistency of the DCBS algorithm applied with such panel data as an input.

We introduce the following notations for presenting the first part of our results. For each $b=1,\ldots,B$, let $\{\mathbf{r}^{b}_{t}\}$ denote a stationary multivariate GARCH($p,q$) process defined with the parameters $\bm{\Omega}_{i}(\eta_{b}+1)$ and $\bm{\Theta}_{i}(\eta_{b}+1)$, and the innovations coinciding with $\bm{\varepsilon}_{t}$ over the associated segment $[\eta_{b}+1,\eta_{b+1}]$; we denote the corresponding conditional variance by $\{\mathbf{h}^{b}_{t}\}$. Then, let $\widetilde{U}^{b}_{i,t}=g_{0}(\mathbf{r}_{i,t}^{b,t-p},\mathbf{h}_{i,t-1}^{b,t-q})$ and $\widetilde{U}^{b}_{ii^{\prime},t}=g_{2}(\mathbf{r}_{i,t}^{b,t-p},\mathbf{h}_{i,t-1}^{b,t-q},\mathbf{r}_{i^{\prime},t}^{b,t-p},\mathbf{h}_{i^{\prime},t-1}^{b,t-q})$, which are constructed analogously as $U_{i,t}$ and $U_{ii^{\prime},t}$ with stationary $\mathbf{r}^{b}_{t}$ and $\mathbf{h}^{b}_{t}$ replacing their nonstationary counterparts in (2) and (4), respectively. Finally, we denote the index of the change point strictly to the left of and nearest to $t$ by $v(t)=\max\{0\leq b\leq B:\,\eta_{b}<t\}$, with which piecewise stationary processes $\mathbf{r}^{v(t)}_{t}$, $\mathbf{h}^{v(t)}_{t}$, $\widetilde{U}^{v(t)}_{i,t}$ and $\widetilde{U}^{v(t)}_{ii^{\prime},t}$ are defined.

Proposition 2**.**

Suppose that (A1)–(A5) hold, and recall the definition of $\{x_{j,t},\,1\leq j\leq d,\,1\leq t\leq T\}$ from (6). Then, the following decomposition holds:

[TABLE]

1. (i)

Setting $f_{j,t}=\widetilde{g}_{ii^{\prime},t}$ (with $\widetilde{g}_{ii,t}\equiv\widetilde{g}_{i,t}$) where $\widetilde{g}_{i,t}=\mathbb{E}\{(\widetilde{U}^{v(t)}_{i,t})^{2}\}$ and $\widetilde{g}_{ii^{\prime},t}=\mathbb{E}(\widetilde{U}^{v(t)}_{ii^{\prime},t})$, we have $f_{j,t}$ piecewise constant with all change points belong to $\mathcal{B}=\{\eta_{1},\ldots,\eta_{B}\}$; conversely, at each $\eta_{b}\in\mathcal{B}$, there exists a non-empty set $\widetilde{\Pi}_{b}\subset\{1,\ldots,d\}$ where

[TABLE] 2. (ii)

The residuals $z_{j,t}=U_{ii^{\prime},t}-\mathbb{E}(\widetilde{U}^{v(t)}_{ii^{\prime},t})$ satisfy

[TABLE]

The proof of Proposition 2 can be found in Appendix A.3. Unlike $\mathbb{E}(U_{i,t}^{2})$ or $\mathbb{E}(U_{ii^{\prime},t})$, we have $\widetilde{g}_{i,t}$ and $\widetilde{g}_{ii^{\prime},t}$ exactly piecewise constant without any boundary effects. By its construction, $z_{j,t}$ do not necessarily satisfy $\mathbb{E}(z_{j,t})=0$. However, due to the mixing properties of the tv-MGARCH processes (Proposition 1), scaled partial sums of $z_{j,t}$ can be appropriately bounded.

In order to establish the theoretical consistency of the DCBS algorithm, we impose the following conditions on the detectability of each change point $\eta_{b}$.

1. (B1)

There exists a fixed constant $c>0$ such that $\min_{0\leq b\leq B}(\eta_{b+1}-\eta_{b})\geq cT^{\gamma}$ for some $\gamma\in(6/7,1]$ (recalling that $\eta_{0}=0$ and $\eta_{B+1}=T$). 2. (B2)

The number of change points, $B\equiv B_{T}$, satisfies $B=o(\log(T))$. 3. (B3)

Recall the definitions of $\widetilde{\Pi}_{b}$ and $\Delta_{j,b}$ given in (11). Then,

[TABLE]

In conjunction with (B1), Assumption (B2) imposes a bound on the total number of change points $B$, which is permitted to grow slowly with $T$. Assumption (B3) specifies the minimum requirement on the cross-sectional size of the change, quantified by $\underline{\Delta}_{d,T}$, for all change points to be detected as well as being located with accuracy. The quantity $\underline{\Delta}_{d,T}$ combines both the cross-sectional ‘density’ of each change point, $|\widetilde{\Pi}_{b}|$, and the magnitude of jumps, $\sum_{j\in\widetilde{\Pi}_{b}}\Delta_{j,b}$, over all $b=1,\ldots,B$; for example, if $\Delta_{j,b}\equiv\Delta_{b}$ for all $j\in\widetilde{\Pi}_{b}$, we have $|\widetilde{\Pi}_{b}|^{-1/2}\sum_{j\in\widetilde{\Pi}_{b}}\Delta_{j,b}=|\widetilde{\Pi}_{b}|^{1/2}\Delta_{b}$, which increases with both $|\widetilde{\Pi}_{b}|$ and $\Delta_{b}$. Typically, change point detection becomes more challenging as the distance between two adjacent change points decreases (with decreasing $\gamma$), and also as the dimensionality $d$ increases relative to the density of the change point $|\widetilde{\Pi}_{b}|$, which is reflected in (B3). We highlight that our methodology does not require each change point to be common to all the cross-sections of $\{x_{j,t}\}$ (and, consequently, in all $\bm{\Omega}_{i}(t)$ and $\bm{\Theta}_{i}(t)$) provided that Assumption (B3) is met. In other words, a relatively sparse change point $\eta_{b}$ (due to small $d^{-1}|\widetilde{\Pi}_{b}|$) is detectable with accuracy as long as $\Delta_{j,b}$ are sufficiently large.

It is not trivial to relate $\Delta_{j,b}$, the magnitude of a jump in $f_{j,t}=\widetilde{g}_{ii^{\prime},t}=\mathbb{E}\{(U^{v(t)}_{ii^{\prime}})^{2}\}$ (with $j=j(i,i^{\prime})$ as defined in (6)) and the changes in $\bm{\Omega}_{i}(t)$ or $\bm{\Theta}_{i}(t)$ due to the presence of the complex transformation. However, the density or sparsity of a change point, measured by $|\Pi_{b}|$ (defined in Assumption (A1) ), is preserved by $|\widetilde{\Pi}_{b}|$ and, in fact, $(i,i^{\prime})\in\Pi_{b}$ iff $j=j(i,i^{\prime})\in\widetilde{\Pi}_{b}$. For further discussion on the high-dimensional efficiency of DC test statistic, we refer to Remark 3.1 of Cho, (2016).

Theorem 1**.**

Suppose that Assumptions (A1)–(A5) and (B1)–(B3) hold. Let $\widehat{\eta}_{b},\ b=1,\ldots,\widehat{B}$ (with $1<\widehat{\eta}_{1}<\ldots<\widehat{\eta}_{\widehat{B}}<T$), denote the change point estimators returned by the DCBS algorithm with a threshold $\pi_{d,T}$ satisfying $C^{\prime}d\underline{\Delta}_{d,T}^{-1}T^{5(1-\gamma)/2}\sqrt{\log(T)}<\pi_{d,T}<C^{\prime\prime}\underline{\Delta}_{d,T}T^{\gamma-1/2}$ for some constants $C^{\prime},C^{\prime\prime}>0$. Then, there exists $c_{0}>0$ such that

[TABLE]

as $T\to\infty$, where $\rho_{d,T}=d\underline{\Delta}_{d,T}^{-2}T^{5(1-\gamma)}\log(T)$.

For the proof of Theorem 1, see Appendix A.4. From the condition imposed on the rate of $\underline{\Delta}_{d,T}$ in (B3), it is easily seen that $\rho_{d,T}/T^{\gamma}\to 0$ as $T\to\infty$. That is, in the re-scaled time interval $[0,1]$, the change point estimators satisfy $T^{-1}|\widehat{\eta}_{b}-\eta_{b}|\leq T^{-\gamma}|\widehat{\eta}_{b}-\eta_{b}|\to 0$ for all $b=1,\ldots,B$. Defining the optimality in change point detection as when each of the true change points and the corresponding estimated change point are within the distance of $O_{p}(1)$ (see e.g. Korostelev, (1987)), it is attained up to a logarithmic factor when the change points are maximally spread ($\gamma=1$), and the jumps are dense ($|\widetilde{\Pi}_{b}|\asymp d$) and of large magnitude ($\sum_{j\in\widetilde{\Pi}_{b}}\Delta_{j,b}\asymp d$).

3.4 Selection of tuning parameters

3.4.1 Choice of parameters for transformation

Empirical performance of the two-stage methodology, its power in particular, is influenced by the choice of the transformation function $g_{0}$ determined by the coefficients $C_{i,j},\,j=0,\ldots,p+q$, in (2). In most references given at the beginning of Section 3.1, $C_{i,j}$ are set as the maximum likelihood estimates (MLEs) of the GARCH parameters obtained from the whole sample assuming the stationarity, say $\widehat{\omega}_{i},\widehat{\alpha}_{i,j},\,1\leq j\leq p$ and $\widehat{\beta}_{i,k},\,1\leq k\leq q$.

The BASTA–res algorithm proposed by Fryzlewicz and Subba Rao, (2014) performs change point detection in the univariate GARCH process by analysing the transformation of the input time series obtained similarly to $U_{i,t}^{2}$. They recommend the use of ‘dampened’ versions of the GARCH parameter estimates. In our setting, this leads to the choice of $C_{i,0}=\widehat{\omega}_{i}$, $C_{i,j}=\widehat{\alpha}_{i,j}/F_{i},\,1\leq j\leq p$, and $C_{i,p+k}=\widehat{\beta}_{i,k}/F_{i},\,1\leq k\leq q$, with within-series dampening parameters $F_{i}\geq 1$.

Empirically, the motivation behind the introduction of $F_{i}$ is as follows. For $r_{i,t}$ with time-varying parameters, we often observe that the GARCH parameters are over-estimated such that $\sum_{j=1}^{p}\widehat{\alpha}_{i,j}+\sum_{k=1}^{q}\widehat{\beta}_{i,k}$ is close to, or even exceeds, one. Therefore, using the raw estimates in place of $C_{i,j}$’s in (2) leads to change points not well-detectable in the resultant transformed panel data. Thus we adopt the dampening parameter $F_{i}$ and select it as

[TABLE]

By construction, $F_{i}$ is bounded as $F_{i}\in[1,99]$ and approximately brings $\widehat{\omega}_{i}$ and $\sum_{j=1}^{p}\widehat{\alpha}_{i,j}+\sum_{k=1}^{q}\widehat{\beta}_{i,k}$ to the same scale.

The transformation $g_{0}$ also involves the unobservable conditional variance $h_{i,t}$, which we propose to replace with the empirical estimates

[TABLE]

obtained with the MLEs of the GARCH parameters.

Typically, the GARCH orders $p$ and $q$ are unknown and may even vary over time. We propose to use $(p,q)=(1,1)$. The GARCH($1$, $1$) model is simple yet known to provide a good fit to a wide range of datasets (see e.g. Hansen and Lunde, (2005)). Besides, the model (1) is adopted for the purpose of change point analysis rather than for describing the time series themselves. In simulation studies reported in Section 4, we study the effect of mis-specifying the GARCH orders on change point analysis, which shows that the choice of $(p,q)=(1,1)$ works reasonably well even when it under-specifies the true GARCH orders.

3.4.2 Choice of threshold for DCBS algorithm

Theorem 1 provides a range of the threshold $\pi_{d,T}$ that guarantees the consistency of the proposed methodology. However, the theoretical range involves typically unattainable knowledge on the quantities such as $\gamma$ or $\underline{\Delta}_{d,T}$. Moreover, even when such knowledge is available, finite sample performance may be affected by the choice of the multiplicative constant to a given rate. Instead, we propose a parametric resampling procedure, which enables us to approximate the distribution of the DC test statistic in the absence of any change point. A similar approach has been widely adopted in the change point literature including Kokoszka and Teyssière, (2002) in the context of testing the presence of a change point in univariate GARCH processes.

Algorithm 2 outlines the proposed resampling scheme, where we derive the segment-dependent threshold $\pi^{(s,e)}_{d,T}$ for each segment $[s,e]$ considered at some iteration of the DCBS algorithm (see Algorithm 1). It takes the empirical residuals $\widehat{\varepsilon}_{i,t}=\widehat{h}_{i,t}^{-1/2}r_{i,t}$ where the conditional variance is estimated as in (12) with the MLE of the GARCH parameters $\widehat{\omega}_{i}$, $\widehat{\alpha}_{i,j}$ and $\widehat{\beta}_{i,k}$. While theoretical investigation into the validity of the bootstrap procedure is beyond the scope of this paper, we verify its good performance on simulated datasets in Section 4. In all our numerical studies, we set $\alpha=0.05$ and $R=100$.

4 Simulation study

4.1 Models

We investigate the numerical performance of the two-stage data segmentation methodology proposed in Section 3 on the datasets simulated from the models below. We are not aware of another change point detection methodology applicable to the segmentation of high-dimensional GARCH processes, which makes a comparative numerical study challenging. Instead, we adopt the simulation models from the literature on change point analysis in univariate GARCH processes; Model (M0) is motivated by the simulation models considered in Fryzlewicz and Subba Rao, (2014), Models (M1)–(M2) by those in Kokoszka and Teyssière, (2002). We note that under (M0), there does not exist any change point, i.e. $B=0$.

1. (M0)

Stationary MGARCH $(1,1)$ processes. Let $\omega_{i}=\omega+\delta_{\omega,i}$, $\alpha_{i,1}=\alpha+\delta_{\alpha,i}$ and $\beta_{i,1}=\beta+\delta_{\beta,i}$, where $\delta_{\cdot,1}\sim_{\text{\scriptsize{i.i.d}}}\mathcal{U}(-\Delta,\Delta)$ for some small $\Delta>0$ is added to each GARCH parameter so that every $r_{i,t}$ has a slightly different set of GARCH parameters. The innovations are generated from two different distributions, namely (i) $\bm{\varepsilon}_{t}\sim_{\text{\scriptsize{i.i.d}}}\mathcal{N}(\mathbf{0},\bm{\Sigma}_{\varepsilon})$ where $\sigma_{i,i^{\prime}}=\varrho^{|i-i^{\prime}|}$ with $\varrho=-0.75$ and (ii) $\varepsilon_{i,t}\sim_{\text{\scriptsize{i.i.d}}}t_{10}$ for each $i$ and $t$. We consider $T=1000$ and $N\in\{50,100\}$.

1. (M0.1)

$(\omega,\alpha_{1},\beta_{1})=(0.4,0.1,0.5)$. 2. (M0.2)

$(\omega,\alpha_{1},\beta_{1})=(0.1,0.1,0.8)$. 2. (M1)

tv-MGARCH $(1,1)$ processes with two change points. We introduce the first change point $\eta_{1}=[T/4]$ to the GARCH parameters such that for a randomly chosen $\mathcal{S}_{1}\subset\{1,\ldots,N\}$, GARCH parameters $\omega_{i}(t),\alpha_{i,1}(t)$ and $\beta_{i,1}(t)$ for $i\in\mathcal{S}_{1}$ change at $t=\eta_{1}$ as $\omega_{i}(t)=\omega^{(1)}\mathbb{I}(t\leq\eta_{1})+\omega^{(2)}\mathbb{I}(t>\eta_{1})+\delta_{\omega,i}$, $\alpha_{i,1}(t)=\alpha^{(1)}_{1}\mathbb{I}(t\leq\eta_{1})+\alpha^{(2)}_{1}\mathbb{I}(t>\eta_{1})+\delta_{\alpha,i}$ and $\beta_{i,1}(t)=\beta^{(1)}_{1}\mathbb{I}(t\leq\eta_{1})+\beta^{(2)}_{1}\mathbb{I}(t>\eta_{1})+\delta_{\beta,i}$, where $\delta_{\cdot,i}$ is as in (M0), and the parameters before and after the change point are chosen as below.

1. (M1.1)

$(\omega,\alpha_{1},\beta_{1}):$ $(0.1,0.3,0.3)$ $\to$ $(0.15,0.25,0.65)$. 2. (M1.2)

$(\omega,\alpha_{1},\beta_{1}):$ $(0.1,0.3,0.3)$ $\to$ $(0.125,0.1,0.6)$. 3. (M1.3)

$(\omega,\alpha_{1},\beta_{1}):$ $(0.1,0.3,0.3)$ $\to$ $(0.15,0.15,0.25)$.

The second change point $\eta_{2}=[3T/5]$ is introduced to the cross-correlation structure. Initially, $\bm{\Sigma}_{\varepsilon}(t)=\bm{\Sigma}_{\varepsilon}$ defined in (M0) for $t\leq\eta_{2}$. Then, for a randomly chosen $\mathcal{S}_{2}\subset\{1,\ldots,N\}$, the rows of $\bm{\Sigma}_{\varepsilon}(t)$ corresponding to $\varepsilon_{i,t},\,i\in\mathcal{S}_{2}$ swap their locations arbitrarily. For the generation of $\bm{\varepsilon}_{t}=(\bm{\Sigma}_{\varepsilon}(t))^{1/2}\mathbf{v}_{t}$, we consider two settings: (i) $v_{i,t}\sim_{\text{\scriptsize{i.i.d}}}\mathcal{N}(0,1)$ for all $i$ and $t$, and (ii) $v_{i,t}\sim_{\text{\scriptsize{i.i.d}}}t_{10}$ for all $i$ and $t$. We set $|\mathcal{S}_{1}|=|\mathcal{S}_{2}|=[\varrho N]$ with $\varrho\in\{1,0.75,0.5,0.25\}$ controlling the ‘sparsity’ of the change points, and consider $N\in\{50,100\}$ and the sample size $T=500$ for (i) and $T=2500$ for (ii); the large sample size under (ii) is to ensure the stability of the quasi-MLE of the GARCH coefficients for the transformation. 3. (M2)

Mis-specification of the orders $p$ and $q$.

1. (M2.1)

Over-specification. The two change points are introduced to tv-MGARCH$(1,1)$ processes with Gaussian innovations as in (M1)(M1.1)–(M1)(M1.2) (referred to as (M2.1.1)–(M2.1.2)), but the GARCH orders are mis-specified in the transformation $g_{0}$ (see (2)) as $(p,q)=(2,2)$. 2. (M2.2)

Under-specification. The two change points are introduced to tv-MGARCH$(2,2)$ processes with Gaussian innovations as in (M1), with the GARCH parameters change at $\eta_{1}$ as:

1. (M2.2.1)

$(\omega,\alpha_{1},\alpha_{2},\beta_{1},\beta_{2})$: $(0.1,0.1,0.2,0.1,0.2)$ $\to$ $(0.15,0.15,0.1,0.35,0.3)$; or 2. (M2.2.2)

$(\omega,\alpha_{1},\alpha_{2},\beta_{1},\beta_{2})$: $(0.1,0.1,0.2,0.1,0.2)$ $\to$ $(0.125,0.1,0,0.3,0.3)$.

The covariance matrix of the innovations change at $\eta_{2}$ as in (M1). The GARCH orders are mis-specified in the transformation $g_{0}$ as $(p,q)=(1,1)$. 4. (M3)

Full-factor multivariate GARCH($1$, $1$) model with time-varying factors and loadings. Proposed in Vrontos et al., (2003), each $r_{i,t}$ is generated as a linear combination of the independent factors $f_{j,t}$, $j=1,\ldots,N$ which are GARCH$(1,1)$ processes.

[TABLE]

with $w_{i,i^{\prime}}\sim_{\text{\scriptsize{i.i.d}}}\mathcal{N}(1,1)$ for the loading matrix $\mathbf{W}$. For a randomly chosen $\mathcal{S}_{1}\subset\{1,\ldots,N\}$, GARCH parameters of $\mathbf{f}$ change at $\eta_{1}=[T/4]$ as

1. (M3.1)

$(\omega,\alpha_{1},\beta_{1}):$ $(0.1,0.3,0.3)$ $\to$ $(0.15,0.25,0.65)$. 2. (M3.2)

$(\omega,\alpha_{1},\beta_{1}):$ $(0.1,0.3,0.3)$ $\to$ $(0.125,0.1,0.6)$.

Another change point is introduced to the loading matrix at $\eta_{2}=[3T/5]$, by swapping the rows of $\mathbf{W}$ corresponding to a randomly chosen $\mathcal{S}_{2}\subset\{1,\ldots,N\}$, which brings in a change in the conditional cross-correlations as well as within-series conditional variance. The cardinality of $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$ is controlled as in (M1) with $\varrho\in\{1,0.75,0.5,0.25\}$, and we consider $T=500$ and $N\in\{50,100\}$.

In all follows, we set $p=q=1$ as suggested in Section 3.4.1 for the transformation $g_{0}$ unless specifically chosen otherwise (see (M2)(M2.2)), and set $\alpha=0.05$ and $R=100$ for threshold selection as described in Section 3.4.2. All simulation results reported are based on $100$ realisations.

4.2 Results

Firstly, we perform at-most-one-change test on the data simulated from (M0) by conducting only a single iteration of the DCBS algorithm in Algorithm 1, which enables us to investigate its size behaviour, see Tables 1. In Appendix B, we investigate the power and localisation accuracy of such a test in various single change point scenarios motivated by Fryzlewicz and Subba Rao, (2014).

Overall, we observe that the DC-based test performs well in size control. When the innovations are generated from a Gaussian distribution with cross-correlations, the test manages to keep the size below the significance level $\alpha=0.05$ when $N=50$, while more spurious false alarms are observed as the dimensionality grows for (M0)(M0.2). We note that, although not directly comparable, Table 1 of Fryzlewicz and Subba Rao, (2014) observed similar size behaviour from their procedure as well as the change point test from Andreou and Ghysels, (2002) applied to univariate GARCH processes generated with the same GARCH parameters as in (M0)(M0.2). When the innovations are drawn from a $t_{10}$-distribution, the parameter configuration of (M0)(M0.2) brings in greater size distortion.

Table 2–3 report the results from applying the proposed methodology to multiple change point detection from tv-MGARCH processes generated as described in (M1), with Gaussian and heavy-tailed $\bm{\varepsilon}_{t}$; Figures 1–2 in Appendix illustrate the locations of the estimated change points. For the change point $\eta_{1}$ attributed solely to the change in GARCH parameters, when it results in time-varying unconditional variance $\text{var}(r_{i,t}),\,i\in\mathcal{S}_{1}$, as in (M1)(M1.1)–(M1)(M1.2) (with (M1)(M1.1) bringing a larger jump in the variance), the change point is more easily detected than when the unconditional variance is kept approximately constant as in (M1)(M1.3). Between the two types of change points considered in (M1), the detection of $\eta_{1}$ is more challenging as it becomes sparser and $N$ increases, compared to $\eta_{2}$ where the cross-sectional correlations undergo a change. This is explained by the fact that, the sparsity of $\eta_{1}$, measured by the number of series containing $\eta_{1}$ scaled by $d$ ($\approx N^{2}$), is in the order of $\varrho N^{-1}$ whereas that of $\eta_{2}$ is in the order of $\varrho^{2}$. With heavier-tailed innovations, both the detection and the localisation accuracy deteriorates, possibly as the quasi-MLE procedure for the first-stage transformation suffers due to the heavy-tailedness and in particular, the detection of $\eta_{1}$ is the most affected between the two.

Comparing Tables 4–5 to Table 2, we observe that mis-specifying the GARCH orders $(p,q)$ does not noticeably worsen the performance of our methodology, confirming the robustness of our methodology to the choice of $p$ and $q$.

For the full-factor MGARCH model in (M3), Table 6 shows that our methodology manages to estimate the two change points with high accuracy, but it also tends to return spurious estimators (see also Figures 7–8 in Appendix). In this scenario, the cross-sectional dependence is strong due to the presence of common factors, which renders the re-sampling procedure for threshold selection discussed in Section 3.4.2 less reliable.

5 Application to financial risk management

In this section, we consider the application of the methodology proposed in Section 3 to risk management using the Value-at-Risk (VaR) of a portfolio, a widely used measure of market risk embraced by financial institutions for regulatory and other internal purposes. In this exercise, we demonstrate the peril of ignoring the change points in the volatility and the cross-sectional correlations of a multi-asset portfolio, and show that the stress period identified by our proposed methodology is effective in providing robust risk management.

5.1 Background: VaR, stressed VaR and its backtests

VaR measures the extreme loss (change in value) of an asset or a portfolio of assets with a prescribed probability level during a given holding period. It has been criticised due to unrealistic assumptions (linearity and normality), parameter sensitivity (to estimation and holding periods) and its inadequacy during crises especially when correlations between assets are observed to vary over time (Persaud, , 2000; Danielsson, , 2002). The last point is of particular interest since, compared to the period of market stability, correlations are observed to be significantly higher when markets are falling (Li et al., , 2017). Jäckel and Rebonato, (2001) noted that a risk manager ‘would greatly over-estimate the degree of diversification in his portfolio in the event of a crash if he used the [correlation] matrix estimated during normal periods’. Works to address the criticism on VaR exist: Valentinyi-Endrész, (2004) examined whether detecting and taking into account change points improves upon VaR forecast. Similarly, Spokoiny, (2009) proposed to perform local change point analysis to detect regions of volatility homogeneity as an alternative to stationary GARCH modelling, and found that the local volatility estimator performed well in the application to VaR.

The Basel Accord (1996 Amendment) requires the use of stressed VaR (sVaR) which is based on a covariance matrix from a crisis period in the past. The accord does not specify the exact time period to be used but instead, proposes the judgement-based and the formulaic approaches (European Banking Authority, , 2012). The former relies on a high-level analysis of the risks related to the holding portfolio, while the latter is a more systematic, quantitative approach where our proposed methodology can contribute by supplying the information about the latest stress periods.

The VaR metric is defined as follows:

[TABLE]

where $F^{-1}(\cdot|\mathcal{G}_{t})$ is the quantile function of the loss and profit distribution with market- and portfolio-specific conditions contained in $\mathcal{G}_{t}$. To backtest VaR by means of statistical tests, we adopt two tests from Kupiec, (1995) using the Proportion of Failure (PoF) and the Time until First Failure (TFF). For a sample of $T$ observations, the test statistics take the form of a likelihood ratio test statistic:

[TABLE]

for some $\alpha\in(0,1)$, where $x_{f}$ denotes the number of failures (losses in excess of the reported VaR) occurred and $t_{f}$ the number of days until the first failure within the $T$ observations occurs. Under suitable assumptions, both $LR_{\text{PoF}}$ and $LR_{\text{TFF}}$ follow a $\chi^{2}_{1}$-distribution asymptotically, which can be used to test whether the VaR model is adequate. These tests do not control for the dependence in the failures, i.e. excess losses beyond the reported VaR may cluster while the overall (unconditional) number of failures is not significantly different from $\alpha T$. To address this limitation, we adopt the dynamic conditional quantile (DQ) test by Engle and Manganelli, (2004) whose test statistic is defined as

[TABLE]

where $\mathbf{Z}$ is the matrix of explanatory variables (e.g. raw and squared past returns) with $\bar{q}=\text{rank}(\mathbf{Z})$, and $\mathbf{H}$ the vector collecting $\text{Hit}_{t}(\alpha):=\mathbb{I}(r_{t}<-\text{VaR}_{t}(\alpha))-\alpha$ over time with $r_{t}$ being the time series of portfolio returns. Later in Sections 5.3.2–5.3.3, we use the three backtests in (13) and (14) when assessing the performance of our proposed methodology in the applications to stressed VaR.

5.2 Change point detection in a multi-asset portfolio

We collect the daily log-returns of the stocks composing the S&P100 index from 1 January 2007 to 31 December 2020 (available from Yahoo! Finance). We winsorize the time series in order to restrict the influence of outliers on the MLE of the GARCH parameters for individual log-return series. In addition, identifying the pairs of stocks whose unconditional correlation is above $0.8$, one of each of the pair is randomly removed, which results in $N=79$ time series ($d_{N}=3160$). We take this step solely for the purpose of ensuring that the estimation of multivariate GARCH processes we consider later in Section 5.3 for out-of-sample validation of the sVaR, is feasible and free of convergence issues; such high correlations are often attributed to the shares of the same company in different classes (e.g. GOOG and GOOGL). We set 75% of the total sample (from 1 January 2007 to 28 April 2016, $T=2347$) as the training data that serves as an input to the proposed change point detection methodology. The remaining observations are used for out-of-sample validation of the sVaR from different stress periods, which is discussed in the next section.

Our method, using the default parameters described in Section 3.4, detects three change points from the training data, which form the four periods reported in Table 7 (Periods 1–4). The table also reports Period 5, which covers the 12-month period following the bankruptcy of Lehman Brothers in September 2008 as the crisis period, with its length chosen in accordance with Basel 2.5 (Basle Committee on Banking Supervision, , 2010). The results indicate that the highest VaR (obtained as quantiles of an equally weighted portfolio of the 79 stocks in modulus) was obtained from Period 2 which spans from 8 September 2008 to 5 May 2009, a rather expected outcome given the high volatility the market experienced after the bankruptcy of Lehman Brothers. We note that this period spans less than 12 months, and is shorter than the maximum allowed duration by Basel 2.5. In addition, it coincides with Bank of England’s view on the historical periods per region with the worst market moves (Bank of England, , 2018). Period 4 also exhibits high stress characteristics but less severe than those observed in the aftermath of the Lehman Brothers bankruptcy. If a bank selects Period 5 (the 95% VaR of which is 4.41%) to calibrate its sVaR model, it will likely under-estimate the measure compared to that calibrated with Period 2 (the 95% VaR of which is 5.32%). In Section 5.3, we explore this argument in further detail.

Out of independent interest, we applied the proposed change point detection methodology to the same dataset covering the recent pandemic crisis, i.e. from 6 August 2009 (after the latest change point identified from the training data) to 31 December 2020. It detects a change point at the beginning of the pandemic located around the end of February 2020, a few days before the Federal Reserve System decided to cut interest rates in order to reduce the economic impacts from the Covid-19 outbreak, see Period 6 and 7 in Table 7. From the same table, we observe that the current pandemic crisis is not, measured by either the 95% VaR or the 95% VaR, as ‘stressed’ as Period 2 coinciding with global financial crisis. One reason might be the extremely fast recovery of the capital markets which likely offset the initial shocks. However, without further data it is not trivial to make a conclusive assessment.

5.3 Out-of-sample performance of sVaR

Our proposed methodology supplies a means to segment the data and identify the period of most volatility, but does not automatically provide an sVaR model. After the change points are detected from a given dataset, a user should decide on the model according to regulatory or internal requirements. In this section, we form optimal portfolios of the S&P100 constituent stocks (Section 5.3.1) and assess which period from Table 7 is the most appropriate for an sVaR purpose using two sVaR models; in order to ensure that we have enough test data for the out-of-sample exercises, only Periods 1–5 have been used for model calibration. One based on the historical simulation (unconditional covariance modelling, Section 5.3.2) and the other based on the DCC model of Engle, (2002) (conditional, Section 5.3.3). The former is the most popular VaR model used by banks and the latter is a dynamic extension of the CCC model (Bollerslev, , 1990), which in turn is closely related to the time-varying model (1) adopted in this paper.

5.3.1 Optimal portfolio formulation

We form portfolios of the $N$ stocks over a given period for testing the out-of-sample performance of the sVaR. In particular, we calculate the portfolio return from the vector of the returns $\mathbf{r}_{t}$ as $r^{\text{opt}}_{t}=\mathbf{w}_{t}^{\top}\bm{\mu}_{t}$, where $\mathbf{w}_{t}\in\mathbb{R}^{N}$ denotes a vector of the portfolio weights and $\bm{\mu}_{t}\in\mathbb{R}^{N}$ a vector of the mean returns of the assets at time $t$. There are many different methods for constructing $\bm{\mu}_{t}$, but we simply calculate the sample mean return for each stock using the observations from the most recent $T_{\text{sv}}\in\{250,500,650\}$ days, which forms a rolling window that moves forward one day at a time. To obtain the optimal weights $\mathbf{w}_{t}$, we solve the following Markowitz optimisation problem

[TABLE]

where $\widehat{\bm{\Sigma}}_{r}(t)$ is an estimator of the covariance matrix of $\mathbf{r}_{t}$; to ensure that it is well-conditioned, we obtain $\widehat{\bm{\Sigma}}_{r}(t)$ via a non-linear shrinkage estimation technique proposed by Ledoit and Wolf, (2012) using the same $T_{\text{sv}}$-day rolling window. The quantity $\mu_{\text{target}}$ denotes a target expected portfolio return that takes its value from $\{0.25\%,0.50\%,0.75\%\}$. We repeatedly solve the optimisation problem every $h_{\text{reb}}\in\{5,20\}$ days keeping the portfolio weights $\mathbf{w}_{t}$ unchanged in between, which implies a weekly or monthly portfolio rebalancing. Following the above steps, we create an optimal portfolio for every triplet $(T_{\text{sv}},\mu_{\text{target}},h_{\text{reb}})$ and use it to evaluate the adequacy of the sVaR as described in Sections 5.3.2 and 5.3.3.

5.3.2 Historical simulation

Let the five periods identified in Table 7 be indexed with $b=1,\ldots,5$. A financial institution typically uses historical simulation over one or two years ($T_{\text{sv}}\in\{250,500\}$), to produce a one-day ahead forecast of sVaR as

[TABLE]

with $\mathbf{L}^{(b)}$ being the Cholesky decomposition of $\widehat{\bm{\Sigma}}^{(b)}_{r}$, an estimator of the (unconditional) covariance matrix for Period $b$ (also obtained using the non-linear shrinkage estimation method of Ledoit and Wolf, (2012), see Section 5.3.1), $\mathbf{L}(t)$ the Cholesky decomposition of $\widehat{\bm{\Sigma}}_{r}(t)$ and $\mathbf{w}_{t}$ the optimal weights obtained from (15). Equation (17) transforms the return vector $\mathbf{r}_{t}$ with covariance matrix $\widehat{\bm{\Sigma}}_{r}(t)$, into the portfolio return $R^{(b)}_{t}$ of a return vector with a (stressed) covariance $\widehat{\bm{\Sigma}}^{(b)}_{r}$, for given $b=1,\ldots,5$ (Duffie and Pan, , 1997). For each $b$, we repeatedly forecast the one-day ahead sVaR of an optimal portfolio using $\widehat{\text{sVaR}}_{t+1}^{(b)}$, until we reach the end of the test data on 31 December 2020. We repeat this exercise for each of the optimal portfolios formed with the triplets $(T_{\text{sv}},\mu_{\text{target}},h_{\text{reb}})\in\{250,500\}\times\{0.25\%,0.50\%,0.75\%\}\times\{5,20\}$.

The results from sVaR backtests are provided in Table 8. In particular, the table indicates that Period 2 can safely be used to calibrate sVaR: The corresponding number of failures is below the expected number of violations for the out-of-sample period at 99% levels (approximately 5 and 7 for $500$ and $750$ days, respectively), and the first failure occurred only after at least 173 days at the 99% level depending on $\mu_{\text{target}}$. Note that the alternative hypothesis in the likelihood ratio tests is two-side and a small number of failures also reject the adequacy of an sVaR model. The sVaR model calibrated using Period 2 yields a favourable, albeit conservative, result and is more likely to be accepted by risk managers. Besides, it passes the traffic light test (Basle Committee on Banking Supervision, , 1996) in all cases, whereby a VaR model is deemed valid (green zone) if the probability of observing up to $x_{f}$ failures is less than $0.95$ under the binomial distribution with $T_{\text{sv}}$ and $\alpha$ as the parameters. Assessing the adequacy of stress periods using the DQ test, Period 2 is the best candidate for the sVaR calibration with $p$-value $>0.1$ in most scenarios except when $\mu_{\text{target}}$ is large and $T_{\text{sv}}=500$.

We mention that Period 5 also exhibits a small number of failures, rather expected behaviour considering that it overlaps with Period 2 for eight months. It is a good candidate for a stress period as it passes almost all tests (except for the DQ test). However, the results could be significantly worse if the actual stress period were shorter, in which case a twelve-month period would be an equal mix of a stress period and a much less stressed period.

5.3.3 Dynamic conditional correlation model

In the second exercise, we adopt the DCC model of Engle, (2002) to model each stationary segment identified by the proposed change point detection methodology. We restrict our analysis to $N=30$ stocks randomly selected from the $79$ stocks in order to avoid convergence issues with the estimation of the DCC model.

Under the DCC, the conditional covariance matrix is factorised as $\mathbf{H}_{t}=\mathbf{D}_{t}\mathbf{R}_{t}\mathbf{D}_{t}$, where $\mathbf{D}_{t}=\text{diag}(h_{1,t}^{1/2},\ldots,h_{N,t}^{1/2})$ and the conditional variance $h_{i,t}$ follows the GARCH($p,q$) model (1) with time-invariant parameters. We consider the GARCH models with $p=q=1$ and normally distributed innovations for their estimation; the results were similar with other model orders $p,q$ or innovation distributions. Engle, (2002) propose the following dynamic correlation structure for $\mathbf{R}_{t}$:

[TABLE]

where $\mathbf{v}_{t}=\mathbf{D}_{t}^{-1}\mathbf{r}_{t}$ is the standardised residual vector and $\bar{\bm{\Sigma}}=[\widehat{\varrho}_{i,i^{\prime}}]_{i,i^{\prime}=1}^{N}$ with $\widehat{\varrho}_{i,i^{\prime}}$ denoting the unconditional sample correlations between $v_{i,t}$ and $v_{i^{\prime},t}$.

Similarly to the approach taken in Section 5.3.2, we consider a portfolio of assets with the return vector $\mathbf{r}_{t}$ and the vector of portfolio weights $\mathbf{w}_{t}$ obtained as in (15) with the triplets $(T_{\text{sv}},\mu_{\text{target}},h_{\text{reb}})\in\{650\}\times\{0.25\%,0.50\%,0.75\%\}\times\{5,20\}$. For each Period $b$, we estimate the parameters $a_{\text{dcc}},b_{\text{dcc}}$ along with $\omega_{i},\alpha_{i,1}$ and $\beta_{i,1}$ (their dependence on $b$ suppressed). Then from the estimated DCC model for the corresponding period, we obtain the one-day ahead forecast of the conditional covariance matrix, $\widehat{\mathbf{H}}^{(b)}_{t+1}$, via (18), produce a forecast of the portfolio volatility as $\widehat{\sigma}^{(b)}_{t+1}=\sqrt{\mathbf{w}_{t}^{\top}\widehat{\mathbf{H}}^{(b)}_{t+1}\mathbf{w}_{t}}$, and compare $-2.33\times\widehat{\sigma}^{(b)}_{t+1}$ (99% sVaR up to its sign) to the actual portfolio returns, until we reach the end of the test data on 31 December 2020. The thus-obtained backtesting results are given in Table 9.

Using the estimated parameters from Period 2 yields the smallest number of failures compared with the rest: At the 99% level, the first failure occurred after 295 to 383 days, compared with 12 days taken for the rest of the stress periods. The DQ test also indicates that Period 2 is the most suitable for calibrating the 99% sVaR metric with $p$-values $>0.10$ in all scenarios. Period 5 fails to pass all the tests as the corresponding numbers of failures exceeded the expected ones. This reinforces our argument that, since Period 5 includes the three month period (May to August 2009) of low stress characteristics, overlapping with Period 3 in Table 7, it is not suitable for calibrating the 99% sVaR metric, see also that Period 3 yields a large number of failures compared to Period 2.

6 Conclusions

In this paper, we propose a two-stage methodology for detecting multiple change points in both within-series and cross-correlation structures of multivariate volatility processes. It first transforms the $N$-dimensional series so that complex structural change points are detectable as change points in the level of $N(N+1)/2$-dimensional transformed data, which provides an input to the multiple change point detection algorithm in the second stage. We show the theoretical consistency of the combined methodology in terms of the total number and locations of estimated change points, and verify its good performance on simulated datasets. Also, we demonstrate the efficacy of the proposed methodology in financial risk management. This exercise shows that, by identifying the period of stress from the dataset of the S&P100 constituent stocks, our method can serve as a formulaic approach to calculating the (stressed) Value-at-Risk of a portfolio of risky assets.

Appendix A Proofs of theoretical results

Notation

In addition to the notations introduced in the main text, $C_{i},\ i=1,2,\ldots$ are adopted to denote fixed positive constants throughout the proofs, which may represent different values from one usage to another. We write $[\mathbf{a}]_{i}$ to denote the $i$-th element of a vector $\mathbf{a}$. Also, $\mathbf{a}\leq\mathbf{b}$ for two vectors $\mathbf{a},\mathbf{b}\in\mathbb{R}^{n}$ indicates that $[\mathbf{a}]_{i}\leq[\mathbf{b}]_{i}$ for all $i=1,\ldots,n$. We write $\mathbf{1}$ to denote a vector of ones, whose dimension should be clear from the context.

A.1 Preliminary results

Recall that $v(t)$ denotes the index of the nearest change point among those satisfying $\eta_{b}<t$, and we denote the location of such a change point by $\underline{\eta}(t)$, i.e., $\underline{\eta}(t)=\max\{\eta_{b}:\,1\leq b\leq B,\,\eta_{b}<t\}=\eta_{v(t)}$. Similarly defined is $\bar{\eta}(t)=\min\{\eta_{b}:\,1\leq b\leq B,\,\eta_{b}\geq t\}$. Also, recall the definitions of $\{\mathbf{r}^{b}_{t}\}$ from above Proposition 2.

Lemma 1**.**

Under Assumptions (A1) and (A3)–(A5), it holds that for each $i$,

[TABLE]

where $\max_{1\leq i\leq N}\mathbb{E}(V_{i,t})\leq c_{1}\rho_{1}^{t-\underline{\eta}(t)-1}$ with some fixed $\rho_{1}\in(0,1)$ and $c_{1}>0$. Similarly, $|h_{i,t}^{2}-(h^{v(t)}_{i,t})^{2}|\leq W_{i,t}$ where $\max_{1\leq i\leq N}\mathbb{E}(W_{i,t})\leq c_{2}\rho_{2}^{t-\underline{\eta}(t)-1}$ with some fixed $\rho_{2}\in(0,1)$ and $c_{2}>0$.

Proof.

The proof below is an extension of the proof given for Lemma 1 of Fryzlewicz and Subba Rao, (2014). For brevity, we omit the subscript $i$ from $r_{i,t}$ and other relevant notations when there is no confusion.

The random recurrence equation for a GARCH($p,q$) process $\mathbf{r}_{t}$ in (1) is defined as below (Chen and An, , 1998):

[TABLE]

where both $\mathcal{A}_{t}$ and $\mathbf{B}_{t}$ are independent over $t$ and have all their elements non-negative. After iterating the equation $k$ steps backwards, we obtain

[TABLE]

(where $\prod_{l=1}^{0}\mathcal{A}_{t-l+1}=1$). As for $r^{v(t)}_{t}$,

[TABLE]

where $\widetilde{\mathbf{Y}}_{t}=((r_{t}^{v(t)})^{2},(r_{t-1}^{v(t)})^{2},\ldots,(r_{t-p+1}^{v(t)})^{2},h^{v(t)}_{t},h^{v(t)}_{t-1},\ldots,h^{v(t)}_{t-q+1})^{\top}$ and $\widetilde{\mathcal{A}}_{t}$ and $\widetilde{\mathbf{B}}_{t}$ are defined analogously as $\mathcal{A}_{t}$ and $\mathbf{B}_{t}$, respectively, with the GARCH parameters associated with the stationary segment $[\underline{\eta}(t)+1,\bar{\eta}(t)]$. By its construction, $r^{v(t)}_{t}$ shares the same GARCH parameters and innovation sequence as $r_{t}$ over $[\underline{\eta}(t)+1,\bar{\eta}(t)]$. Therefore setting $k=t-\underline{\eta}(t)-1$ leads to $II=IV$. As $r_{t}^{2}-(r^{v(t)}_{t})^{2}$ is equal to the first element of $\mathbf{Y}_{t}-\widetilde{\mathbf{Y}}_{t}$, we have $|r_{t}^{2}-(r^{v(t)}_{t})^{2}|\leq V_{t}=[\prod_{j=0}^{k-1}\mathcal{A}_{t-j}\mathbf{Y}_{t-k}]_{1}+[\prod_{j=0}^{k-1}\widetilde{\mathcal{A}}_{t-j}\widetilde{\mathbf{Y}}_{t-k}]_{1}$. Since $\mathcal{A}_{t-j},\,j=0,\ldots,k-1$ and $\mathbf{Y}_{t-k}$ ($\widetilde{\mathcal{A}}_{t-j},\,j=0,\ldots,k-1$ and $\widetilde{\mathbf{Y}}_{t-k}$) are independent, and $\mathcal{A}^{*}=\mathbb{E}(\mathcal{A}_{t})=\mathbb{E}(\widetilde{\mathcal{A}}_{t})$ remains the same over $t\in[\underline{\eta}(t)+1,\bar{\eta}(t)]$, $\mathbb{E}(V_{t})=\left[(\mathcal{A}^{*})^{k}\mathbb{E}(\mathbf{Y}_{t-k})\right]_{1}+\left[(\mathcal{A}^{*})^{k}\mathbb{E}(\widetilde{\mathbf{Y}}_{t-k})\right]_{1}$. Further, as the entries of $\mathcal{A}^{*}$ are non-negative, if $\mathbf{x}$ and $\mathbf{y}$ are $(p+q)$-vectors satisfying $\mathbf{0}\leq\mathbf{x}\leq\mathbf{y}$, we have $[(\mathcal{A}^{*})^{k}\mathbf{x}]_{1}\leq[(\mathcal{A}^{*})^{k}\mathbf{y}]_{1}$. Under (A3)–(A4), both $\mathbb{E}(\mathbf{Y}_{t-k})$ and $\mathbb{E}(\widetilde{\mathbf{Y}}_{t-k})$ are bounded from above by $\epsilon_{2}^{-1}\Xi_{1}\mathbf{1}$ and hence $\mathbb{E}(V_{t})\leq 2\epsilon_{2}^{-1}\Xi_{1}[(\mathcal{A}^{*})^{k}\mathbf{1}]_{1}$. Under (A4), it is clear that $(\mathcal{A}^{*})^{p+q}\mathbf{1}\leq C_{1}(1-\epsilon_{2})\mathbf{1}$ for some $C_{1}$. Therefore, we can find for some fixed constant $c_{1}>0$ and $\rho_{1}\in(0,1)$ such that

[TABLE]

Since $\rho_{1}$ depends on $p$, $q$, $\epsilon_{2}$ and $\Xi_{1}$ only, the same arguments hold for all $i=1,\ldots,N$ and hence the bound on $\mathbb{E}(V_{i,t})$ holds uniformly in $i$ for an appropriately chosen $\rho_{1}$. As for $W_{t}$, note that $h_{t}^{2}-(h^{v(t)}_{t})^{2}$ is equal to the $(p+1)$-th element of $\mathbf{Y}_{t}-\widetilde{\mathbf{Y}}_{t}$. Analogous arguments show that there exists $c_{2}>0$, $\rho_{2}\in(0,1)$ such that $\mathbb{E}(W_{i,t})\leq c_{2}\rho_{2}^{t-\underline{\eta}(t)-1}$. ∎

A.2 Proof of Proposition 1

Fryzlewicz and Subba Rao, (2011) showed that time-varying, univariate ARCH processes are strong mixing at a geometric rate. As the techniques used by the authors for the proof of their Theorem 3.1 are applicable to other Markov processes, we extend their results by deriving the mixing rate for time-varying bivariate GARCH processes. The arguments adopted below for showing the mixingness of $\mathbf{r}_{12,t}=(r_{1,t},r_{2,t})^{\top}$, are analogously applied to showing that of any pair of $r_{i,t}$ and $r_{i^{\prime},t},\,1\leq i<i^{\prime}\leq N$.

Let $\mathbf{X}_{t}=(X_{1,t},X_{2,t})^{\top}$ with $X_{i,t}=r_{i,t}^{2}$ and $\mathbf{Z}_{t}=(Z_{1,t},Z_{2,t})^{\top}$ with $Z_{i,t}=\varepsilon_{i,t}^{2}$ for $i=1,2$, and $R=p\vee q$. Recall (A5), which is met by the density function of $\mathbf{Z}_{t}$. Denote the $\sigma$-algebra generated by $\mathbf{X}_{t},\ldots,\mathbf{X}_{t+k}$ by $\mathcal{F}_{t+k}^{t}=\sigma(\mathbf{X}_{t},\ldots,\mathbf{X}_{t+k})$. Also, for any $k>0$, let $\underline{\mathbf{X}}_{t}^{t-k}=(\mathbf{X}_{t}^{\top},\ldots,\mathbf{X}_{t-k}^{\top})^{\top}$ and $\underline{\mathbf{X}}_{t}^{t-k}=0$ for $k\leq 0$. We similarly define $\underline{\mathbf{Z}}_{t}^{t-k}=(\mathbf{Z}_{t}^{\top},\ldots,\mathbf{Z}_{t-k}^{\top})^{\top}$, $\underline{\mathbf{X}}_{i,t}^{t-k}=(X_{i,t},\ldots,X_{i,t-k})^{\top}$, $\underline{\mathbf{h}}_{i,t}^{t-k}=(h_{i,t},\ldots,h_{i,t-k})^{\top}$ and $\underline{\mathbf{Z}}_{i,t}^{t-k}=(Z_{i,t},\ldots,Z_{i,t-k})^{\top}$.

Recall the random recurrence equation in the proof of Lemma 1. We define $\mathbf{Y}_{i,t}$, $\mathcal{A}_{i,t}$ and $\mathbf{B}_{i,t}$ similarly as $\mathbf{Y}_{t}$, $\mathcal{A}_{t}$ and $\mathbf{B}_{t}$ of (19) with respect to the innovations $Z_{i,t}$ and the time-varying GARCH parameters associated with $r_{i,t}$ for $i=1,2$, except that they are now defined with $\alpha_{i,j}(t),\beta_{i,k}(t),\,j,k=1,\ldots,R$ (with $\alpha_{i,j}(t)=0,\,j>p$ and $\beta_{i,k}(t)=0,\,k>q$ whenever relevant) such that their dimensions are $2R$ ($\mathbf{Y}_{i,t}$, $\mathbf{B}_{i,t}$) and $(2R)\times(2R)$ ($\mathcal{A}_{i,t}$) with appropriately padded zeros. Finally, let $\mathbf{Y}_{t}=(\mathbf{Y}_{1,t}^{\top},\mathbf{Y}_{2,t}^{\top})^{\top}$, $\mathbf{B}_{t}=(\mathbf{B}_{1,t}^{\top},\mathbf{B}_{2,t}^{\top})^{\top}$ and $\mathcal{A}_{t}$ be the block diagonal matrix with $\mathcal{A}_{1,t}$ and $\mathcal{A}_{2,t}$ in its diagonal. Then it can be shown that

[TABLE]

By definition, $\mathcal{A}_{t}$ and $\mathbf{B}_{t}$ are independent random matrices and vectors over $t$, respectively. Let $\mathcal{A}^{*}_{t}=\mathbb{E}(\mathcal{A}_{t})$ and $\mathcal{A}^{*}_{i,t}=\mathbb{E}(\mathcal{A}_{i,t})$. For the proof of Proposition 1, we need the following Lemmas.

Lemma 2**.**

For $s\geq 0$ and $k\geq R$, we have

[TABLE]

where $\underline{\mathbf{V}}_{i,t}^{t-l}=((\underline{\mathbf{X}}_{i,t}^{t-l})^{\top},(\underline{\mathbf{h}}_{i,t}^{t-l})^{\top})^{\top}$ for $l>0$, and $\underline{\mathbf{V}}_{i,t}^{t-l}=0$ otherwise.

(i) When $s=0$: we have $\underline{\mathbf{V}}_{i,t+k+s-1}^{t+k}=0$ and

[TABLE]

(ii) When $1\leq s<R$:

[TABLE]

(iii) When $s\geq R$:

[TABLE]

Proof.

The proof of lemma follows trivially from the representation in (29). ∎

From Lemma 2, the conditional density of $\mathbf{X}_{t+k+s}$ given $\underline{\mathbf{V}}_{t+k+s-1}^{t+k}=((\underline{\mathbf{V}}_{1,t+k+s-1}^{t+k})^{\top},(\underline{\mathbf{V}}_{2,t+k+s-1}^{t+k})^{\top})^{\top}$, $\underline{\mathbf{Z}}_{t+k-1}^{t+1}=((\mathbf{Z}_{1,t+k-1}^{t+1})^{\top},(\mathbf{Z}_{2,t+k-1}^{t+1})^{\top})^{\top}$ and $\underline{\mathbf{V}}_{t}^{t-R+1}=((\underline{\mathbf{V}}_{1,t}^{t-R+1})^{\top},(\underline{\mathbf{V}}_{2,t}^{t-R+1})^{\top})^{\top}$, is a function of

[TABLE]

That is, denoting the conditional density by $f_{s,k,t}$,

[TABLE]

The following lemma is adapted from Proposition 2.1 of Fryzlewicz and Subba Rao, (2011).

Lemma 3**.**

Provided its existence, denote the conditional density of $\underline{\mathbf{X}}_{t+k+R_{2}}^{t+k}$ given $\underline{\mathbf{V}}_{t}^{t-R_{1}}$, by $f_{\underline{\mathbf{X}}_{t+k+R_{2}}^{t+k}|\underline{\mathbf{V}}_{t}^{t-R_{1}}}$. For any $\bm{\gamma}=(\gamma_{1},\ldots,\gamma_{4R_{1}+4})^{\top}\in(\mathbb{R}^{+})^{4R_{1}+4}$, define the set $E=\{\omega:\,\underline{\mathbf{V}}_{t}^{t-R_{1}}(\omega)\in\mathcal{E}\}$ where $\mathcal{E}=\{(\nu_{1},\ldots,\nu_{4R_{1}+4})^{\top}:\,|\nu_{l}|\leq\gamma_{l}\mbox{ for all }l=1,\ldots,4R_{1}+4\}$. Let $\underline{\mathbf{W}}$ denote a random vector independent of $\underline{\mathbf{V}}_{t}^{t-R_{1}}$ and denote its density function by $f_{\underline{\mathbf{W}}}$. Recalling the definition of $f_{s,k,t}$ in (30), define

[TABLE]

where $\mathbf{v}_{s-1}$, $\mathbf{z}$ and $\mathbf{v}$ are vectors of appropriate dimensions.

Then, for any $R_{1},R_{2}\geq 0$, the followings hold:

(i)

[TABLE]

(ii)

[TABLE]

Proof.

The proof given in Fryzlewicz and Subba Rao, (2011) for their Proposition 2.1 is readily applicable to show the above, since $\underline{\mathbf{h}}^{t-R_{1}}_{t}\in\sigma(\underline{\mathbf{X}}^{t-R_{1}}_{t})$ and $\underline{\mathbf{h}}^{t+k}_{t+k+s-1}\in\sigma(\underline{\mathbf{X}}^{t+k}_{t+k+s-1})$. ∎

Lemma 4**.**

For any $\bm{\gamma}\in(\mathbb{R}^{+})^{4R}$, we have

[TABLE]

Proof.

It holds that

[TABLE]

for any $R^{\prime}\geq R=p\vee q$. Hence

[TABLE]

and applying Lemma 3 with $R_{1}=R_{2}=R$ completes the proof. ∎

Lemma 5**.**

If (A5) holds, then for any positive $A_{i}$ and $B_{i}$, $i=1,2$, we have

[TABLE]

Proof.

Observe that the LHS of the above inequality is bounded from the above by

[TABLE]

To bound $I$, changing the variables with $u/(A_{1}+B_{1})=u^{\prime}$ and $v/(A_{2}+B_{2})=v^{\prime}$,

[TABLE]

It is trivial to show that $II\leq B_{1}/(A_{1}+B_{1})$, and the same arguments apply to bound $III$ and $IV$, which completes the proof. ∎

Lemma 6**.**

Under (A3)–(A5), for all $\mathbf{v}\in\mathcal{E}$, we have

[TABLE]

where $\mathbf{v}_{i},\,i=1,2$ denote the sub-vectors of $\mathbf{v}$ of the equal length, and $\epsilon_{3}\in(0,\epsilon_{2})$ for $\epsilon_{2}$ defined in (A4).

Proof.

When $\mathbf{v}_{i}=\mathbf{0}$, we have $\mathcal{Q}^{i}_{s,k,t}(\mathbf{0},\underline{\mathbf{Z}}_{i,t+k-1}^{t+1})=0$ and

[TABLE]

Recalling (30), we have

[TABLE]

Applying Lemma 5,

[TABLE]

and hence

[TABLE]

Let $\bar{\gamma}=\max_{1\leq l\leq 4R}\gamma_{l}$. From Lemmas 1 and 2, the following holds for any $\mathbf{v}\in\mathcal{E}$:

[TABLE]

for some $C_{4}>0$ and $\epsilon_{3}\in(0,\epsilon_{2})$ which, when plugged into (33), yields (32). ∎

Now we are fully equipped to prove Proposition 1. Lemma 6 provides the upper bound on the first term in (31). Therefore it only remains to bound the probability terms in (31). Using Markov’s inequality,

[TABLE]

for some fixed $C_{5}>0$, and $\mathbb{P}(h_{i,t-l}\geq\gamma)$ is similarly bounded. Thus,

[TABLE]

Setting $\gamma_{l}=(1-\epsilon_{3})^{-k/2}$ for all $l=1,\ldots,4R$, we obtain

[TABLE]

for some $\alpha\in[\sqrt{1-\epsilon_{3}},1)$ and fixed $M>0$. Since $\alpha$ does not depend on the index of $\mathbf{r}_{12}$, the same arguments apply to any $\mathbf{r}_{ii^{\prime},t},\,1\leq i<i^{\prime}\leq N$ and hence the proposition is proved. $\square$

A.3 Proof of Proposition 2

By its definition, $g_{0}$ is bounded and and is Lipschitz continuous in its squared arguments, i.e. there exist $\bar{g},C_{g}\in(0,\infty)$ such that

[TABLE]

Part (i) follows trivially from the definition of $\widetilde{U}^{v(t)}_{i,t}$ and $\widetilde{U}^{v(t)}_{ii^{\prime},t}$. More specifically, since $r^{v(t)}_{i,t}$ ‘jumps’ from one stationary process $r^{b}_{i,t}$ to another $r^{b+1}_{i,t}$ without any boundary effect, we have the expectations of $(\widetilde{U}^{v(t)}_{i,t})^{2}$ and $\widetilde{U}^{v(t)}_{ii^{\prime},t}$ exactly piecewise constant with their change points coinciding with $\mathcal{B}$. For the proof of (ii), we adopt the arguments similar to those used in the proof of Lemma 4 in Fryzlewicz and Subba Rao, (2014).

Recall that $v(t)$ denotes the index of the nearest change point among those satisfying $\eta_{b}<t$, and the definition of $z_{j,t}$ in the context of time series segmentation, namely $z_{j,t}=U_{i,t}^{2}-\widetilde{g}_{i,t}$ or $z_{j,t}=U_{ii^{\prime},t}-\widetilde{g}_{ii^{\prime},t}$ for some $i,i^{\prime}\in\{1,\ldots,N\}$. Note that the claim of the proposition follows when $\mathbb{P}(\mathcal{E}_{1}\cap\mathcal{E}_{2})\to 1$, where

[TABLE]

for some fixed $c^{\prime},c^{\prime\prime}>0$.

We first focus on the case when $z_{j,t}=U_{i,t}^{2}-\widetilde{g}_{i,t}$. We investigate the probability of the following event:

[TABLE]

for $\lambda_{T}\asymp\sqrt{\log(T)}$. Let $\nu=e-s+1$.

(a) When $\nu<\lambda_{T}^{2}/(4\bar{g}^{4})$: Note that $|z_{j,t}|\leq|U_{i,t}^{2}|+|\widetilde{g}_{i,t}|\leq 2\bar{g}^{2}$ from (34). Hence $\nu^{-1/2}|\sum_{t=s}^{e}z_{t}|\leq\nu^{1/2}2\bar{g}^{2}<\lambda_{T}$.

(b) When $\nu>\lambda_{T}^{2}/(4\bar{g}^{4})$: Decompose $z_{t}$ as

[TABLE]

Also from (34),

[TABLE]

for $\rho=\rho_{1}\vee\rho_{2}\in(0,1)$ and fixed $C_{1}>0$ thanks to Lemma 1. Therefore

[TABLE]

for some $C_{2}>0$. Then, the probability of the event in (35) is bounded from the above by $\mathbb{P}\{\nu^{-1/2}|\sum_{t=s}^{e}\xi_{i,t}|\geq\widetilde{\lambda}_{T}\}$ with $\widetilde{\lambda}_{T}=\lambda_{T}-C_{2}(p+q)B\nu^{-1/2}$. From Proposition 1, $U_{i,t}$ is a bounded, strong mixing process with its $\alpha$-mixing coefficient $\alpha(k)\asymp\alpha^{k}$ for some $\alpha\in(0,1)$. Combining the mixing property of $\xi_{i,t}$ with the fact that $|\xi_{i,t}|\leq 2\bar{g}^{2}$, we apply Theorem 1 of Merlevède et al., (2009) and derive

[TABLE]

where $C_{2}$ depends only on $\alpha$ and hence not on $j$ ($i$), $s$ or $e$. Since $\widetilde{\lambda}_{T}\to\infty$ as $T\to\infty$, while $\widetilde{\lambda}_{T}<\lambda_{T}$ and $\nu\geq\lambda_{T}^{2}/(4\bar{g}^{4})$, we have the RHS of (36) bounded by $C_{3}\exp(-C_{4}\widetilde{\lambda}_{T}^{2})$ for some fixed $C_{3},C_{4}>0$. Summarising (a)–(b) above, it can be shown that

[TABLE]

as $T\to\infty$ under (A2) and (B2), for a sufficiently chosen $c^{\prime}$.

We now turn our attention to $\mathbb{P}(\mathcal{E}_{2}|\mathcal{E}_{1})$. Let $z_{j,t}=U_{ii^{\prime},t}-\widetilde{g}_{ii^{\prime},t}$. By construction,

[TABLE]

Given $\mathcal{E}_{1}$, we only need to investigate the boundedness of the last term. Observe that

[TABLE]

Since $g_{3}$ (defined in Appendix A.5) is also Lipschitz continuous in its squared arguments, we can bound $\nu^{-1/2}|\sum_{t=s}^{e}\zeta^{\prime}_{ii^{\prime},t}|$ similarly as $\nu^{-1/2}|\sum_{t=s}^{e}\zeta_{i^{\prime},t}|$. As for $\nu^{-1/2}|\sum_{t=s}^{e}\xi^{\prime}_{ii^{\prime},t}|$ since $\xi^{\prime}_{ii^{\prime},t}$ is a zero-mean, strong mixing process from Proposition 1, its boundedness analogously follows as that of $\nu^{-1/2}|\sum_{t=s}^{2}\xi_{i,t}|$, and $\mathbb{P}(\mathcal{E}_{2}|\mathcal{E}_{1})\to 1$, which completes the proof.

A.4 Proof of Theorem 1

Lemma 7**.**

Under the conditions of Proposition 2,

[TABLE]

where

[TABLE]

Proof.

Note that

[TABLE]

From Lemma 1,

[TABLE]

and similarly $VI=O_{p}(\sqrt{\log(T)}),$, both of which hold uniformly in $j$ and $(s,c,e)$. Therefore, for some fixed $c>0$,

[TABLE]

with probability tending to one as $T\to\infty$. ∎

By the definition of $g_{0}$ and the discussion in Appendix A.5, there exists a fixed constant $\bar{f}>0$ such that $\max_{1\leq j\leq d}\max_{1\leq t\leq T}|f_{j,t}|\leq\bar{f}$. Proposition 2 and Lemma 7 place a logarithmic bound on the partial sums and the cumulative sums of $z_{j,t}$. With these results in place of Lemmas 1–2 of Cho, (2016), we are able to apply the arguments identical to those employed there to prove Theorem 1 under (B1)–(B3).

A.5 Properties of $g_{2}$

We show that from (34), $g_{2}$ is bounded and Lipschitz continuous in its squared arguments. Firstly, note that

[TABLE]

and hence $|g_{2}(\mathbf{z},\mathbf{z}^{\prime})|\leq|g_{1}(\mathbf{z})|+|g_{1}(\mathbf{z}^{\prime})|+2|g_{0}(\mathbf{z})g_{0}(\mathbf{z}^{\prime})|\leq 4\bar{g}$ for any $\mathbf{z},\mathbf{z}^{\prime}\in\mathbb{R}^{p+q+1}$. From (37), we only need to establish the Lipschitz continuity of $g_{3}(\mathbf{z},\mathbf{z}^{\prime})=g_{0}(\mathbf{z})g_{0}(\mathbf{z}^{\prime})$. Further introducing the notations $\mathbf{w},\mathbf{w}^{\prime}\in\mathbb{R}^{p+q+1}$,

[TABLE]

and thus follows the Lipschitz continuity of $g_{2}$ in its squared arguments.

Appendix B Additional simulation results

B.1 Power of the test

In this section, we provide additional simulation results complementing the simulation studies conducted in Section 4 of the main text.

1. (M4)

tv-MGARCH ($1$, $1$) processes with a single change point. Change points are introduced to GARCH parameters in (1) at $\eta_{1}\in\{[T/2],[9T/10]\}$, where $T=1000$ and $N\in\{50,100\}$. For a randomly chosen $\mathcal{S}_{1}\subset\{1,\ldots,N\}$, GARCH parameters $\omega_{i}(t),\alpha_{i,1}(t)$ and $\beta_{i,1}(t)$ for $i\in\mathcal{S}_{1}$ change at $t=\eta_{1}$ as $\omega_{i}(t)=\omega^{(1)}\mathbb{I}(t\leq\eta_{1})+\omega^{(2)}\mathbb{I}(t>\eta_{1})+\delta_{\omega,i}$, $\alpha_{i,1}(t)=\alpha^{(1)}_{1}\mathbb{I}(t\leq\eta_{1})+\alpha^{(2)}_{1}\mathbb{I}(t>\eta_{1})+\delta_{\alpha,i}$ and $\beta_{i,1}(t)=\beta^{(1)}_{1}\mathbb{I}(t\leq\eta_{1})+\beta^{(2)}_{1}\mathbb{I}(t>\eta_{1})+\delta_{\beta,i}$, where $|\mathcal{S}_{1}|=[\varrho N]$ with $\varrho\in\{1,0.75,0.5,0.25\}$ controlling the ‘sparsity’ of the change point. We have $\delta_{\cdot,i}\sim_{\text{\scriptsize{i.i.d}}}\mathcal{U}(-\Delta,\Delta)$ is as in (M0), and $\bm{\varepsilon}_{t}\sim_{\text{\scriptsize{i.i.d}}}\mathcal{N}(\mathbf{0},\bm{\Sigma}_{\varepsilon}(t))$ with $\bm{\Sigma}_{\varepsilon}(t)=\bm{\Sigma}_{\varepsilon}$ defined in (M0).

1. (M4.1)

$(\omega,\alpha_{1},\beta_{1}):$ $(0.4,0.1,0.5)$ $\to$ $(0.4,0.1,0.6)$. 2. (M4.2)

$(\omega,\alpha_{1},\beta_{1}):$ $(0.4,0.1,0.5)$ $\to$ $(0.4,0.1,0.8)$. 3. (M4.3)

$(\omega,\alpha_{1},\beta_{1}):$ $(0.1,0.1,0.8)$ $\to$ $(0.1,0.1,0.7)$. 4. (M4.4)

$(\omega,\alpha_{1},\beta_{1}):$ $(0.1,0.1,0.8)$ $\to$ $(0.1,0.1,0.4)$. 5. (M4.5)

$(\omega,\alpha_{1},\beta_{1}):$ $(0.4,0.1,0.5)$ $\to$ $(0.5,0.1,0.5)$. 6. (M4.6)

$(\omega,\alpha_{1},\beta_{1}):$ $(0.4,0.1,0.5)$ $\to$ $(0.8,0.1,0.5)$. 7. (M4.7)

$(\omega,\alpha_{1},\beta_{1}):$ $(0.1,0.1,0.8)$ $\to$ $(0.3,0.1,0.8)$. 8. (M4.8)

$(\omega,\alpha_{1},\beta_{1}):$ $(0.1,0.1,0.8)$ $\to$ $(0.5,0.1,0.8)$.

Tables 10–13 present the results on the power of the test from conducting a single iteration of the DCBS algorithm. As expected, the test achieves higher power when the change point is cross-sectionally dense (with larger $\varrho$) and located centrally, and the same applies to the localisation accuracy of the change point estimator. For most GARCH parameter configurations, the test attains power above $0.9$ even when the change point is relatively sparse ($\varrho=0.25$) except for (M4)(M4.2), (M4)(M4.4) and (M4)(M4.8). Even when the location of the change point is skewed ($\eta_{1}=[9T/10]$), our method generally attains high power and localisation accuracy if the change point is not too sparse, in most settings where it shows good performance when the change point is centrally located ($\eta_{1}=[T/2]$).

B.2 Localisation accuracy

Figures 1–8 illustrate the locations of estimated change points from the simulation studies conducted in Section 4 of the main text.