Bayesian MIDAS Penalized Regressions: Estimation, Selection, and Prediction

TL;DR

This paper introduces a Bayesian Group Lasso approach with spike-and-slab priors for high-dimensional mixed-frequency regressions, enhancing prediction accuracy and variable selection, especially in economic forecasting contexts.

Contribution

It develops a novel Bayesian penalized regression method with adaptive hyper-parameter tuning, establishing theoretical properties and demonstrating superior empirical performance.

Findings

Good selection and forecasting in small samples

Outperforms competitors in U.S. GDP forecasting

Limited short-term predictive power of financial variables

Abstract

We propose a new approach to mixed-frequency regressions in a high-dimensional environment that resorts to Group Lasso penalization and Bayesian techniques for estimation and inference. In particular, to improve the prediction properties of the model and its sparse recovery ability, we consider a Group Lasso with a spike-and-slab prior. Penalty hyper-parameters governing the model shrinkage are automatically tuned via an adaptive MCMC algorithm. We establish good frequentist asymptotic properties of the posterior of the in-sample and out-of-sample prediction error, we recover the optimal posterior contraction rate, and we show optimality of the posterior predictive density. Simulations show that the proposed models have good selection and forecasting performance in small samples, even when the design matrix presents cross-correlation. When applied to forecasting U.S. GDP, our penalized regressions can outperform many strong competitors. Results suggest that financial variables may have some, although very limited, short-term predictive content.

Figures

Tables

Table 1. Table 1: Monte Carlo simulations: estimation and selection accuracy

$K$	${\sigma }_{\epsilon }$	$\overline{\sigma }$	MSE	MSE$\left(𝒜\right)$	MSE$\left({𝒜}^c\right)$	VAR	BIAS2	TPR	FPR	MCC
			DGP 1
			BMIDAS-AGL
30	0.50	1.3	4.4E-03	1.7E-02	1.9E-03	1.4E-03	3.0E-03	0.96	0.03	0.90
	0.95	2.1	6.5E-02	2.2E-01	3.3E-02	1.3E-02	5.2E-02	0.36	0.04	0.41
50	0.50	1.3	3.8E-03	1.9E-02	2.0E-03	1.1E-03	2.7E-03	0.94	0.03	0.85
	0.95	2.1	4.6E-02	2.3E-01	2.5E-02	9.2E-03	3.6E-02	0.35	0.03	0.40
			BMIDAS-AGL-SS
30	0.50	1.3	3.3E-03	1.6E-02	7.1E-04	1.3E-03	1.9E-03	0.96	0.01	0.95
	0.95	2.1	6.3E-02	2.7E-01	2.2E-02	1.6E-02	4.8E-02	0.36	0.02	0.49
50	0.50	1.3	2.2E-03	1.8E-02	4.7E-04	8.9E-04	1.3E-03	0.94	0.00	0.95
	0.95	2.1	3.8E-02	2.7E-01	1.2E-02	8.8E-03	2.9E-02	0.35	0.01	0.50
			DGP 2
			BMIDAS-AGL
30	0.50	1.1	2.3E-03	8.7E-03	1.1E-03	8.3E-04	1.5E-03	0.99	0.02	0.93
	0.95	1.8	4.5E-02	1.8E-01	1.9E-02	6.9E-03	3.9E-02	0.48	0.04	0.53
50	0.50	1.1	2.0E-03	1.0E-02	1.1E-03	6.3E-04	1.4E-03	0.98	0.02	0.90
	0.95	1.8	3.1E-02	1.8E-01	1.5E-02	5.1E-03	2.6E-02	0.46	0.03	0.51
			BMIDAS-AGL-SS
30	0.50	1.1	1.6E-03	7.7E-03	3.7E-04	7.8E-04	8.2E-04	0.99	0.01	0.98
	0.95	1.8	5.4E-02	2.4E-01	1.7E-02	1.2E-02	4.2E-02	0.47	0.02	0.58
50	0.50	1.1	1.2E-03	9.0E-03	2.9E-04	5.4E-04	6.2E-04	0.98	0.00	0.97
	0.95	1.8	3.3E-02	2.5E-01	9.5E-03	7.4E-03	2.6E-02	0.45	0.01	0.60
			DGP 3
			BMIDAS-AGL
30	0.50	1.0	1.2E-02	3.9E-02	6.5E-03	1.3E-03	1.1E-02	0.85	0.15	0.61
	0.95	1.6	9.5E-02	2.5E-01	6.5E-02	1.1E-02	8.4E-02	0.35	0.10	0.27
50	0.50	1.0	1.4E-02	6.0E-02	8.8E-03	1.2E-03	1.3E-02	0.77	0.16	0.46
	0.95	1.6	8.2E-02	2.7E-01	6.1E-02	9.4E-03	7.2E-02	0.32	0.10	0.22
			BMIDAS-AGL-SS
30	0.50	1.0	7.3E-03	3.1E-02	2.5E-03	1.3E-03	6.0E-03	0.87	0.06	0.78
	0.95	1.6	6.8E-02	2.7E-01	2.7E-02	1.1E-02	5.7E-02	0.33	0.03	0.41
50	0.50	1.0	6.6E-03	4.1E-02	2.8E-03	1.0E-03	5.6E-03	0.81	0.05	0.71
	0.95	1.6	4.3E-02	2.8E-01	1.7E-02	7.4E-03	3.6E-02	0.30	0.02	0.41
Notes: BMIDAS-AGL and BMIDAS-AGL-SS refer to the models described in Section 3. $\overline{\sigma }$ is the average variance of the error process across MC simulations. MSE, VAR, and BIAS2 denote the Mean Squared Error, the Variance, and the Squared Bias, respectively. $\text{MSE}\left(𝒜\right)$ and $\text{MSE}\left({𝒜}^c\right)$ denote the MSE for the active and inactive set, respectively. TPR, FPR, and MCC denote the True Positive Rate, the False Positive Rate, and the Matthews Correlation Coefficient, respectively.

Table 2. Table 2: Monte Carlo simulations: restricted vs unrestricted weighting schemes

		BMIDAS-AGL			BMIDAS-AGL-SS
$K$	${\sigma }_{\epsilon }$	TPR	FPR	MCC	TPR	FPR	MCC
		DGP 1
30	0.50	0.15	-0.02	0.15	0.14	-0.01	0.11
	0.95	0.06	0.01	0.03	0.06	-0.01	0.08
50	0.50	0.25	-0.04	0.27	0.15	-0.01	0.12
	0.95	0.08	0.00	0.09	0.06	0.00	0.09
		DGP 2
30	0.50	0.09	-0.02	0.09	0.08	-0.01	0.07
	0.95	0.10	-0.01	0.12	0.11	-0.01	0.13
50	0.50	0.18	-0.03	0.21	0.11	-0.01	0.11
	0.95	0.14	-0.01	0.15	0.09	0.00	0.11
		DGP 3
30	0.50	-0.10	0.11	-0.27	-0.08	0.04	-0.15
	0.95	-0.09	0.06	-0.20	-0.10	0.01	-0.12
50	0.50	-0.12	0.12	-0.33	-0.11	0.04	-0.19
	0.95	-0.08	0.07	-0.23	-0.10	0.01	-0.13
Notes: See Table 1. The reported values denote the difference in TPR, FPR, and MCC between models estimated with restricted $\left(r=2\right)$ and unrestricted $\left(r=0\right)$ Almon lag polynomials.

Table 3. Table 3: Monte Carlo simulations: forecasting performance

	RMSFE	LS	CRPS	RMSFE	LS	CRPS	RMSFE	LS	CRPS	RMSFE	LS	CRPS
$K$	$30$						$50$
${\sigma }_{\epsilon }$	0.50			0.95			0.50			0.95
	DGP 1
Oracle	$1.17$	-$1.61$	$0.68$	$1.87$	-$2.13$	$1.10$	$1.12$	-$1.59$	$0.65$	$1.88$	-$2.12$	$1.10$
BMIDAS-AGL	1.15	-0.11	1.14	1.11	-0.04	1.09	1.22	-0.15	1.20	1.18	-0.09	1.15
BMIDAS-AGL-SS	1.13	-0.09	1.12	1.10	-0.04	1.08	1.14	-0.10	1.13	1.15	-0.07	1.12
BMIDAS-AL	1.29	-0.22	1.28	1.19	-0.11	1.16	1.53	-0.38	1.50	1.31	-0.20	1.27
MIDAS-L	1.29	-0.22	1.28	1.18	-0.09	1.14	1.59	-0.42	1.58	1.35	-0.23	1.31
MIDAS-EN	1.32	-0.24	1.31	1.17	-0.09	1.14	1.64	-0.45	1.63	1.37	-0.25	1.33
MIDAS-SCAD	1.31	-0.25	1.31	1.20	-0.11	1.17	1.75	-0.61	1.77	1.39	-0.30	1.36
MIDAS-MC$+$	1.39	-0.31	1.40	1.18	-0.09	1.15	1.69	-0.53	1.69	1.41	-0.32	1.39
	DGP 2
Oracle	$0.97$	-$1.42$	$0.56$	$1.55$	-$1.91$	$0.90$	$0.93$	-$1.40$	$0.54$	$1.51$	-$1.88$	$0.87$
BMIDAS-AGL	1.08	-0.06	1.07	1.09	-0.05	1.08	1.14	-0.10	1.13	1.14	-0.07	1.11
BMIDAS-AGL-SS	1.07	-0.05	1.06	1.09	-0.06	1.08	1.11	-0.07	1.10	1.12	-0.07	1.10
BMIDAS-AL	1.13	-0.08	1.11	1.19	-0.12	1.16	1.23	-0.16	1.20	1.27	-0.17	1.23
MIDAS-L	1.10	-0.06	1.09	1.08	-0.04	1.07	1.14	-0.09	1.12	1.13	-0.07	1.10
MIDAS-EN	1.11	-0.07	1.10	1.08	-0.04	1.07	1.16	-0.11	1.14	1.13	-0.07	1.11
MIDAS-SCAD	1.08	-0.05	1.07	1.11	-0.06	1.10	1.11	-0.06	1.09	1.22	-0.14	1.20
MIDAS-MC$+$	1.08	-0.04	1.06	1.09	-0.04	1.08	1.11	-0.06	1.09	1.18	-0.11	1.16
	DGP 3
Oracle	$0.95$	-$1.41$	$0.55$	$1.44$	-$1.86$	$0.85$	$0.91$	-$1.39$	$0.53$	$1.52$	-$1.88$	$0.88$
BMIDAS-AGL	1.20	-0.13	1.18	1.17	-0.09	1.14	1.26	-0.18	1.24	1.21	-0.16	1.21
BMIDAS-AGL-SS	1.16	-0.12	1.15	1.16	-0.09	1.14	1.20	-0.13	1.17	1.16	-0.11	1.15
BMIDAS-AL	1.37	-0.26	1.34	1.26	-0.16	1.22	1.54	-0.43	1.51	1.34	-0.26	1.33
MIDAS-L	1.37	-0.28	1.36	1.33	-0.21	1.30	1.36	-0.26	1.34	1.30	-0.24	1.29
MIDAS-EN	1.38	-0.28	1.36	1.31	-0.19	1.27	1.40	-0.29	1.37	1.30	-0.22	1.28
MIDAS-SCAD	1.38	-0.30	1.36	1.30	-0.20	1.27	1.50	-0.43	1.50	1.45	-0.45	1.46
MIDAS-MC$+$	1.42	-0.32	1.41	1.27	-0.16	1.23	1.47	-0.39	1.46	1.37	-0.32	1.36
Notes: RMSFE, LS, and CRPS denote respectively the root mean squared forecast error, the log-score, and the continuously ranked probability score, in relative terms with respect to the Oracle benchmark. Values for the Oracle (in italic) denote absolute outcomes. MIDAS regressions with Lasso, Elastic-Net, SCAD, and MC$+$ penalties are estimated using the penalized MATLAB toolbox. Regularization parameters are tuned by optimizing the BIC over a fine grid of values.

Table 4. Table 4: Out-of-sample results: relative accuracy of point and density forecasts

	RMSFE	LS	CRPS	RMSFE	LS	CRPS	RMSFE	LS	CRPS	RMSFE	LS	CRPS
	$h=0$			$h=1/3$			$h=2/3$			$h=1$
RW	$0.028$	$2.138$	$0.016$	$0.028$	$2.138$	$0.016$	$0.028$	$2.138$	$0.016$	$0.029$	$2.069$	$0.016$
BMIDAS-AGL	0.60	0.49	0.60	0.53	0.59	0.53	0.65	0.45	0.64	0.78	0.31	0.76
BMIDAS-AGL-SS	0.59	0.55	0.58	0.51	0.70	0.50	0.55	0.62	0.54	0.72	0.40	0.71
AR(1)	0.85	0.17	0.81	0.85	0.17	0.81	0.85	0.17	0.81	0.84	0.19	0.80
BMIDAS-AL	0.59	0.50	0.59	0.56	0.55	0.57	0.70	0.40	0.68	0.82	0.28	0.79
MIDAS-L	0.80	0.23	0.77	0.75	0.29	0.74	0.83	0.19	0.80	0.80	0.24	0.77
MIDAS-EN	0.80	0.23	0.77	0.75	0.29	0.74	0.83	0.19	0.80	0.80	0.25	0.77
MIDAS-SCAD	0.63	0.48	0.62	0.76	0.34	0.73	0.74	0.30	0.73	0.80	0.24	0.77
MIDAS-MC$+$	0.68	0.36	0.69	0.67	0.44	0.66	0.68	0.40	0.67	0.81	0.30	0.77
BMIDAS-comb	0.64	0.45	0.64	0.66	0.45	0.65	0.68	0.44	0.65	0.78	0.33	0.75
BMA-MIDAS	0.61	0.52	0.60	0.59	0.56	0.58	0.55	0.61	0.54	0.81	0.25	0.81
BMS-MIDAS	0.62	0.51	0.60	0.59	0.54	0.59	0.59	0.57	0.58	0.82	0.23	0.83
Factor-BMIDAS PCA	0.67	0.43	0.67	0.67	0.46	0.66	0.66	0.48	0.65	0.77	0.35	0.75
Factor-BMIDAS QML	0.71	0.39	0.70	0.73	0.39	0.71	0.71	0.42	0.69	0.88	0.29	0.81
$t$-Factor-BMIDAS PCA	0.66	0.45	0.66	0.62	0.50	0.62	0.65	0.46	0.66	0.74	0.36	0.74
$t$-Factor-BMIDAS QML	0.63	0.49	0.62	0.61	0.54	0.60	0.64	0.50	0.63	0.74	0.33	0.76
	$h=4/3$			$h=5/3$			$h=2$			$h=4$
RW	$0.029$	$2.069$	$0.016$	$0.029$	$2.069$	$0.016$	$0.033$	$1.967$	$0.018$	$0.036$	$1.872$	$0.020$
BMIDAS-AGL	0.91	0.20	0.89	1.16	0.08	1.05	0.87	0.18	0.86	1.05	-0.03	1.05
BMIDAS-AGL-SS	0.83	0.24	0.82	0.82	0.24	0.80	0.76	0.28	0.74	0.77	0.29	0.76
AR(1)	0.84	0.19	0.80	0.84	0.19	0.80	0.79	0.24	0.76	0.73	0.31	0.70
BMIDAS-AL	1.01	0.23	0.90	1.32	0.10	1.06	0.84	0.23	0.82	1.01	-0.64	1.00
MIDAS-L	0.83	0.19	0.80	0.82	0.22	0.78	0.79	0.24	0.76	0.71	0.33	0.68
MIDAS-EN	0.83	0.19	0.80	0.82	0.22	0.78	0.79	0.25	0.76	0.71	0.32	0.69
MIDAS-SCAD	0.83	0.19	0.80	0.82	0.22	0.78	0.79	0.24	0.76	0.71	0.33	0.68
MIDAS-MC$+$	0.76	0.29	0.75	0.85	0.21	0.80	0.79	0.24	0.76	0.72	0.33	0.69
BMIDAS-comb	0.81	0.24	0.80	0.88	0.17	0.86	0.82	0.21	0.81	0.72	0.35	0.71
BMA-MIDAS	0.84	0.22	0.83	0.82	0.13	0.78	0.81	0.24	0.79	0.75	0.31	0.75
BMS-MIDAS	0.85	0.19	0.84	0.89	-0.02	0.86	0.81	0.22	0.79	0.75	0.30	0.74
Factor-BMIDAS PCA	0.86	0.27	0.81	0.96	0.21	0.87	0.92	0.18	0.85	0.83	0.20	0.81
Factor-BMIDAS QML	0.91	0.22	0.87	0.95	0.18	0.91	0.97	0.13	0.93	0.91	0.14	0.88
$t$-Factor-BMIDAS PCA	0.76	0.32	0.75	0.81	0.25	0.79	0.74	0.30	0.72	0.72	0.32	0.70
$t$-Factor-BMIDAS QML	0.75	0.32	0.76	0.78	0.29	0.78	0.74	0.31	0.74	0.73	0.32	0.70
Notes: RMSFE, LS, and CRPS denote respectively the root mean squared forecast error, the log-score, and the continuously ranked probability score, in relative terms with respect to the random-walk (RW) benchmark. Values for the RW (in italic) denote absolute outcomes. Bold values denote the best outcomes.

Table 5. Table 5: Out-of-sample results: DMW test for equal predictive accuracy with respect to the BMIDAS-AGL-SS model

	RMSFE	LS	CRPS	RMSFE	LS	CRPS	RMSFE	LS	CRPS	RMSFE	LS	CRPS
	$h=0$			$h=1/3$			$h=2/3$			$h=1$
RW	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
AR(1)	0.03	0.00	0.00	0.01	0.00	0.00	0.01	0.00	0.00	0.18	0.02	0.14
BMIDAS-AL	0.49	0.12	0.27	0.02	0.00	0.00	0.00	0.00	0.00	0.08	0.02	0.06
MIDAS-L	0.03	0.00	0.01	0.00	0.00	0.00	0.02	0.00	0.00	0.22	0.05	0.21
MIDAS-EN	0.03	0.00	0.01	0.00	0.00	0.00	0.02	0.00	0.00	0.23	0.05	0.21
MIDAS-SCAD	0.03	0.03	0.03	0.00	0.00	0.00	0.00	0.00	0.00	0.22	0.05	0.21
MIDAS-MC$+$	0.01	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.19	0.12	0.18
BMIDAS-comb	0.04	0.00	0.01	0.00	0.00	0.00	0.02	0.00	0.00	0.23	0.14	0.28
BMA-MIDAS	0.12	0.16	0.13	0.00	0.00	0.00	0.33	0.31	0.43	0.05	0.03	0.04
BMS-MIDAS	0.11	0.11	0.11	0.00	0.00	0.00	0.06	0.08	0.06	0.03	0.02	0.02
Factor-BMIDAS PCA	0.03	0.01	0.01	0.00	0.00	0.00	0.01	0.00	0.00	0.19	0.15	0.21
Factor-BMIDAS QML	0.00	0.00	0.00	0.01	0.00	0.00	0.00	0.00	0.00	0.09	0.11	0.10
$t$-Factor-BMIDAS PCA	0.08	0.06	0.03	0.01	0.00	0.00	0.01	0.01	0.00	0.36	0.27	0.31
$t$-Factor-BMIDAS QML	0.09	0.08	0.08	0.00	0.00	0.00	0.01	0.00	0.00	0.34	0.13	0.19
	$h=4/3$			$h=5/3$			$h=2$			$h=4$
RW	0.02	0.00	0.00	0.02	0.00	0.00	0.06	0.05	0.01	0.06	0.02	0.03
AR(1)	0.47	0.25	0.59	0.39	0.17	0.45	0.23	0.09	0.24	0.66	0.61	0.78
BMIDAS-AL	0.15	0.45	0.18	0.09	0.07	0.05	0.00	0.14	0.00	0.00	0.10	0.00
MIDAS-L	0.49	0.26	0.61	0.52	0.32	0.60	0.24	0.12	0.27	0.81	0.72	0.88
MIDAS-EN	0.49	0.26	0.61	0.53	0.33	0.61	0.25	0.14	0.29	0.76	0.68	0.85
MIDAS-SCAD	0.49	0.26	0.61	0.52	0.32	0.60	0.24	0.12	0.27	0.81	0.72	0.88
MIDAS-MC$+$	0.88	0.82	0.87	0.34	0.31	0.47	0.22	0.08	0.22	0.70	0.68	0.81
BMIDAS-comb	0.65	0.49	0.63	0.16	0.08	0.13	0.13	0.13	0.08	0.82	0.83	0.82
BMA-MIDAS	0.37	0.26	0.34	0.52	0.22	0.60	0.27	0.36	0.22	0.81	0.89	0.84
BMS-MIDAS	0.19	0.11	0.20	0.13	0.12	0.15	0.26	0.24	0.23	0.82	0.72	0.83
Factor-BMIDAS PCA	0.34	0.68	0.55	0.17	0.36	0.24	0.19	0.21	0.17	0.22	0.16	0.29
Factor-BMIDAS QML	0.16	0.40	0.19	0.12	0.22	0.09	0.14	0.19	0.09	0.09	0.07	0.12
$t$-Factor-BMIDAS PCA	0.90	0.93	0.92	0.61	0.55	0.53	0.87	0.66	0.76	0.76	0.68	0.80
$t$-Factor-BMIDAS QML	0.92	0.87	0.89	0.73	0.74	0.62	0.65	0.66	0.46	0.68	0.72	0.79
Notes: $p$-values of the DMW test for unconditional equal predictive accuracy of model $i$ compared to the BMIDAS-AGL-SS model, according to the one-sided $t$-statistic version of the test. Bold values denote $p$-values $\le 10\%$.

Equations

$$y_{t}=\alpha+\sum_{k=1}^{K}\sum_{c=0}^{C-1}B\left(c;\boldsymbol{\theta}_{k}\right)L^{\scalebox{0.55}{$\left.c\middle/m\right.$}}x_{k,t-h}^{(m)}+\epsilon_{t},$$

$$y_t=\alpha +{\mathbf{\theta }}^{\prime }{\mathbf{Z}}_{t-h}^{(m)}+{ϵ}_t,$$

$${\hat{y}}_{T|T-h}=\hat{\alpha }+{\hat{\mathbf{\theta }}}^{\prime }{\mathbf{Z}}_{T-h}^{(m)},$$

$${\mathcal{Q}}_{\text{AGL}}(\mathbf{\theta }):=T^{-1}{\mathcal{L}}_T(\mathbf{\theta })+\sum _{j=1}^G{\lambda }_j\Vert {\mathbf{\theta }}_j{\Vert }_2,$$

$$\Pi ({\mathbf{\theta }}_j|{\sigma }^2)$$

$$\propto {\left(\frac{{\lambda }_j^2}{{\sigma }^2}\right)}^{\frac{g_j}{2}}\exp \left(-\frac{{\lambda }_j}{\sigma }\Vert {\mathbf{\theta }}_j{\Vert }_2\right).$$

$$y|\mathbf{Z},\mathbf{\theta },{\sigma }^2$$

$${\mathbf{\theta }}_j|{\tau }_j^2,{\sigma }^2$$

$${\tau }_j^2$$

$${\sigma }^2$$

$$\Pi (\mathbf{\theta },\mathbf{\tau },{\sigma }^2|\mathbf{\lambda },y,\mathbf{Z})$$

$$\times \prod _{j=1}^G{\left(\frac{1}{2\pi {\sigma }^2{\tau }_j^2}\right)}^{\frac{g_j}{2}}\exp \left(-\frac{\Vert {\mathbf{\theta }}_j{\Vert }_2^2}{2{\sigma }^2{\tau }_j^2}\right)$$

$$\times \prod _{j=1}^G{\left({\lambda }_j^2\right)}^{\frac{g_j+1}{2}}{\left({\tau }_j^2\right)}^{\frac{g_j+1}{2}-1}\exp \left(-\frac{{\lambda }_j^2}{2}{\tau }_j^2\right).$$

$$y|\mathbf{Z},\mathbf{\theta },{\sigma }^2$$

$${\mathbf{\theta }}_j|{\tau }_j^2,{\sigma }^2,{\pi }_0$$

$${\tau }_j^2$$

$${\sigma }^2$$

$${\pi }_0$$

$$\Pi (\mathbf{\theta },\mathbf{\tau },{\sigma }^2,{\pi }_0|\mathbf{\lambda },y,\mathbf{Z})$$

$$\times \prod _{j=1}^G\left[(1-{\pi }_0){\left(\frac{1}{2\pi {\sigma }^2{\tau }_j^2}\right)}^{\frac{g_j}{2}}\exp \left(-\frac{\Vert {\mathbf{\theta }}_j{\Vert }_2^2}{2{\sigma }^2{\tau }_j^2}\right){\mathbf{I}}_{{\mathbf{\theta }}_j\ne 0}+{\pi }_0{\delta }_0({\mathbf{\theta }}_j)\right]$$

$$\times \prod _{j=1}^G{\left({\lambda }_j^2\right)}^{\frac{g_j+1}{2}}{\left({\tau }_j^2\right)}^{\frac{g_j+1}{2}-1}\exp \left(-\frac{{\lambda }_j^2}{2}{\tau }_j^2\right).$$

$$\Vert \mathbf{Z}{\Vert }_0\underline{c}/\sqrt{T}\le {\lambda }_{\min }\le {\lambda }_{\max }\le \overline{C}\min \left\{\sqrt{\log (T)},\sqrt{\log (G)}\Vert \mathbf{Z}{\Vert }_o\right\}$$

$$\underset{{\mathbf{\theta }}_0\in {\overline{\Theta }}_0,{\sigma }_0^2\in [{\underline{\sigma }}^2,{\overline{\sigma }}^2]}{\sup }{\mathbf{E}}_0\left[\Pi \left(\mathbf{\theta }\in {\mathbb{R}}^{\tilde{g}};\Vert \mathbf{Z}(\mathbf{\theta }-{\mathbf{\theta }}_0){\Vert }_2^2\ge MT{ϵ}^2|y,\mathbf{Z}\right)\right]\to 0.$$

$$\tilde{\varphi }(s):=\inf \left\{\frac{\Vert \mathbf{Z}\mathbf{\theta }{\Vert }_2^2}{\Vert \mathbf{Z}{\Vert }_o^2\Vert \mathbf{\theta }{\Vert }_2^2},0\le s_{\mathbf{\theta }}\le s\right\}.$$

$$\underset{{\mathbf{\theta }}_0\in {\overline{\Theta }}_0,{\sigma }_0^2\in [{\underline{\sigma }}^2,{\overline{\sigma }}^2]}{\sup }{\mathbf{E}}_0\left[\Pi \left(\mathbf{\theta }\in {\mathbb{R}}^{\tilde{g}};\Vert \mathbf{\theta }-{\mathbf{\theta }}_0{\Vert }_2^2\ge \frac{MT{ϵ}^2}{\tilde{\varphi }(s_0+M_2{\tilde{s}}_0)\Vert \mathbf{Z}{\Vert }_o^2}|y,\mathbf{Z}\right)\right]\to 0.$$

$$f_{{\mathbf{Z}}_{\tau -h},y,\mathbf{Z}}(y_{\tau }):=f(y_{\tau }|{\mathbf{Z}}_{\tau -h},y,\mathbf{Z}):=\int f_0(y_{\tau }|\mathbf{\theta },{\sigma }^2,{\mathbf{Z}}_{\tau -h})\Pi (\mathbf{\theta },{\sigma }^2|y,\mathbf{Z})d\mathbf{\theta }d{\sigma }^2,\forall \tau >T,$$

$${\eta }_0:=\underset{\{\mathbf{\theta };0\le s_{\mathbf{\theta }}\le M_2{\tilde{s}}_0+s_0\}}{\sup }\frac{|{\mathbf{Z}}_{\tau -h}^{\prime }\mathbf{\theta }{|}^2}{{\sum }_{j=1}^G\Vert {\mathbf{Z}}_j^{\prime }{\mathbf{\theta }}_j/\sqrt{T}{\Vert }_2^2}$$

$$\underset{{\mathbf{\theta }}_0\in {\overline{\Theta }}_0,{\sigma }_0^2\in [{\underline{\sigma }}^2,{\overline{\sigma }}^2]}{\sup }{\mathbf{E}}_0\left[\Pi \left(\mathbf{\theta }\in {\mathbb{R}}^{\tilde{g}};|{\mathbf{Z}}_{\tau -h}^{\prime }(\mathbf{\theta }-{\mathbf{\theta }}_0){|}^2\ge {\eta }_0\frac{M{ϵ}^2}{\tilde{\varphi }(s_0+M_2{\tilde{s}}_0)}|y,\mathbf{Z}\right)\right]\to 0.$$

$$\underset{{\mathbf{\theta }}_0\in {\overline{\Theta }}_0,{\sigma }_0^2\in [{\underline{\sigma }}^2,{\overline{\sigma }}^2]}{\sup }{\mathbf{E}}_{{\mathbf{Z}}_{\tau -h}}\left[{\mathbf{E}}_0\left[\Vert f_{{\mathbf{Z}}_{\tau -h},y,\mathbf{Z}}-f_{{\mathbf{Z}}_{\tau -h},{\mathbf{\theta }}_0,{\sigma }_0^2}{\Vert }_{TV}\right]\right]\to 0.$$

$$\begin{gathered} \int \Pi ({\mathbf{\theta }}_0,{\sigma }_0^2)d{\mathbf{\theta }}_{0}d{\sigma }_0^2{\mathbf{E}}_0\left[K\left(f_{{\mathbf{Z}}_{\tau -h},{\mathbf{\theta }}_0,{\sigma }_0^2}\Vert {\tilde{f}}_{{\mathbf{Z}}_{\tau -h},y,\mathbf{Z}}(y_{\tau })\right)-K\left(f_{{\mathbf{Z}}_{\tau -h},{\mathbf{\theta }}_0,{\sigma }_0^2}\Vert f_{{\mathbf{Z}}_{\tau -h},y,\mathbf{Z}}(y_{\tau })\right)\right] \\ =\int \Pi ({\mathbf{\theta }}_0,{\sigma }_0^2)d{\mathbf{\theta }}_{0}d{\sigma }_0^2\int f_{\mathbf{Z},{\mathbf{\theta }}_0,{\sigma }_0^2}(y)dy\int f_{{\mathbf{Z}}_{\tau -h},{\mathbf{\theta }}_0,{\sigma }_0^2}(y_{\tau })\log \frac{f_{{\mathbf{Z}}_{\tau -h},y,\mathbf{Z}}(y_{\tau })}{{\tilde{f}}_{{\mathbf{Z}}_{\tau -h},y,\mathbf{Z}}(y_{\tau })}d{y}_{\tau } \\ =\int f_{\mathbf{Z}}(y)dy\int f_{{\mathbf{Z}}_{\tau -h},y,\mathbf{Z}}(y_{\tau })\log \frac{f_{{\mathbf{Z}}_{\tau -h},y,\mathbf{Z}}(y_{\tau })}{{\tilde{f}}_{{\mathbf{Z}}_{\tau -h},y,\mathbf{Z}}(y_{\tau })}d{y}_{\tau }, \end{gathered}$$

$${\mathbf{\theta }}_j|{\mathbf{\theta }}_{\text{leavevmode}to2.45pt\text{vbox}to4.45pt\text{pgfpicture}\text{makeatletter}\text{lower}-0.22499ptto0.0pt\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{definecolor}pgfstrokecolorrgb\mathrm{0,0,0}\text{pgfsys@color@rgb@stroke}000\text{pgfsys@invoke}\text{pgfsys@color@rgb@fill}000\text{pgfsys@invoke}\text{pgfsys@setlinewidth}0.4pt\text{pgfsys@invoke}\text{nullfont}to0.0pt\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{pgfsys@setlinewidth}0.45pt\text{pgfsys@invoke}\text{pgfsys@roundcap}\text{pgfsys@invoke}\text{pgfsys@moveto}2.0pt0.0pt\text{pgfsys@lineto}0.0pt4.0pt\text{pgfsys@stroke}\text{pgfsys@invoke}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{hss}\text{pgfsys@discardpath}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{hss}\text{lxSVG@closescope}\text{endpgfpicture}j},{\sigma }^2,\mathbf{\tau },\mathbf{\lambda },y,\mathbf{Z}$$

$${\tau }_j^{-2}|\mathbf{\theta },{\sigma }^2,\mathbf{\lambda },y,\mathbf{Z}$$

$${\sigma }^2|\mathbf{\theta },\mathbf{\tau },\mathbf{\lambda },y,\mathbf{Z}$$

$${\lambda }_j^2|\mathbf{\theta },{\sigma }^2,\mathbf{\tau },y,\mathbf{Z}$$

$${\mathbf{\theta }}_j|{\mathbf{\theta }}_{\text{leavevmode}to2.45pt\text{vbox}to4.45pt\text{pgfpicture}\text{makeatletter}\text{lower}-0.22499ptto0.0pt\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{definecolor}pgfstrokecolorrgb\mathrm{0,0,0}\text{pgfsys@color@rgb@stroke}000\text{pgfsys@invoke}\text{pgfsys@color@rgb@fill}000\text{pgfsys@invoke}\text{pgfsys@setlinewidth}0.4pt\text{pgfsys@invoke}\text{nullfont}to0.0pt\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{pgfsys@setlinewidth}0.45pt\text{pgfsys@invoke}\text{pgfsys@roundcap}\text{pgfsys@invoke}\text{pgfsys@moveto}2.0pt0.0pt\text{pgfsys@lineto}0.0pt4.0pt\text{pgfsys@stroke}\text{pgfsys@invoke}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{hss}\text{pgfsys@discardpath}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{hss}\text{lxSVG@closescope}\text{endpgfpicture}j},{\sigma }^2,\mathbf{\tau },\mathbf{\lambda },{\mathbf{\pi }}_1,{\pi }_0,y,\mathbf{Z}$$

$${\tau }_j^{-2}|\mathbf{\theta },{\sigma }^2,\mathbf{\lambda },{\mathbf{\pi }}_1,{\pi }_0,y,\mathbf{Z}$$

$${\sigma }^2|\mathbf{\theta },\mathbf{\tau },\mathbf{\lambda },{\mathbf{\pi }}_1,{\pi }_0,y,\mathbf{Z}$$

$${\lambda }_j^2|\mathbf{\theta },{\sigma }^2,\mathbf{\tau },{\mathbf{\pi }}_1,{\pi }_0,y,\mathbf{Z}$$

$${\pi }_0|\mathbf{\theta },{\sigma }^2,\mathbf{\tau },\mathbf{\lambda },{\mathbf{\pi }}_1,y,\mathbf{Z}$$

$${\gamma }_j$$

Full text

1 Introduction

Mixed-data sampling (MIDAS) regressions provide a parsimonious and theoretically efficient treatment of the time-aggregation problem when data are sampled at different frequencies (Ghysels et al., 2005, 2007; Andreou et al., 2010). They have been therefore intensively and successfully used to forecast low frequency series, such as GDP, using monthly, weekly or daily predictors (Clements and Galvão, 2008, 2009; Kuzin et al., 2011; Andreou et al., 2013). However, one important issue with MIDAS regressions is the estimation and prediction in high dimension, as the presence of many high-frequency covariates may easily lead to overparameterized models, in-sample overfitting, and poor predictive performance. This happens because the MIDAS approach can efficiently address the dimensionality issue arising from the number of high-frequency lags in the model, but not that arising from the number of high-frequency variables. In this respect, a number of strategies have been proposed in the literature, such as unrestricted MIDAS regressions (U-MIDAS; Foroni et al., 2015) with automatic model selection of relevant predictors and high-frequency lags (Castle and Hendry, 2010; Bec and Mogliani, 2015), factor-MIDAS regressions (Marcellino and Schumacher, 2010), and targeted factor-MIDAS regressions (Bai and Ng, 2008; see also Bessec, 2013, Bulligan et al., 2015, and Girardi et al., 2017).

Recently, the literature has been increasingly focusing on Machine Learning and penalized regression techniques for macroeconomic applications in a high-dimensional environment (see Korobilis, 2013; Ng, 2013; Gefang, 2014; Ng, 2017; Koop et al., 2019; Korobilis and Pettenuzzo, 2019; to name only a few). Nevertheless, only a few contributions have paid attention to MIDAS regressions in high-dimension so far (Marsilli, 2014; Siliverstovs, 2017; Uematsu and Tanaka, 2019; Babii et al., 2020).111Marsilli (2014) proposes a functional MIDAS combined with a Lasso objective function, which is solved in 1-step through a non-linear optimization algorithm. Siliverstovs (2017) proposes a 2-step targeted factor-MIDAS approach, where the soft-thresholding rule is built around U-MIDAS regressions combined with an Elastic-Net objective function. Very recently, Babii et al. (2020) have proposed a sparse-group Lasso estimator for MIDAS regression based on Legendre polynomials, which is closely related to our approach. In particular, Uematsu and Tanaka (2019) propose a general theoretical framework for penalized regressions where the number of covariates diverges sub-exponentially from the number of observations. Their framework is then especially fitted for U-MIDAS regressions, as in this case the number of parameters to estimate grows with both the number of high-frequency regressors and the length of the unrestricted lag function. However, this approach might not be generally suited for macroeconomic applications with monthly or daily high-frequency predictors, as penalized U-MIDAS regressions can saturate in presence of more predictors than observations, while unrestricted lags of the high-frequency predictors, by construction highly correlated, may be subject to random selection, leaving most of the remaining lag coefficients shrunk to zero.

In the present paper, we address these issues by proposing a novel penalized Bayesian MIDAS approach, based on Almon lag polynomials (a classic approximating function of the distributed lag model) and a prior that induces a Group Lasso penalty. The group penalty seems particularly attractive for the MIDAS framework, as one can assign the lag polynomial of a given predictor into one single group. Shrinkage is then performed simultaneously over the entire lag polynomial, rather than on individual terms of the MIDAS weighting function, overcoming the problem of extremely high correlation within the distributed lags. We then introduce two models: the Bayesian MIDAS Adaptive Group Lasso (BMIDAS-AGL) and the Bayesian MIDAS Adaptive Group Lasso with spike-and-slab prior (BMIDAS-AGL-SS). The latter combines the penalized likelihood approach of the Group Lasso prior with a spike-and-slab prior at the group level, which is expected to improve the sparsity recovery ability of the model (Zhang et al., 2014; Zhao and Sarkar, 2015; Xu and Ghosh, 2015; Ročková and George, 2018). Finally, penalty hyper-parameters governing the group selection are here tuned through a computationally efficient and data-driven approach based on stochastic approximations (Atchadé, 2011; Atchadé et al., 2011) that does not resort to extremely time consuming pilot runs.

We theoretically validate our procedures by establishing good frequentist asymptotic properties of both the posterior predictive error (in- and out-of-sample) and the posterior predictive distribution. For this purpose, we focus on the BMIDAS-AGL-SS model and we adopt a frequentist point of view, i.e. we admit the existence of a true parameter that generates the data. The asymptotic theory is developed by allowing the in-sample size $T$, the number of groups $G$, and their maximum size $g_{\max}$, to increase to infinity. First, we establish consistency and derive the posterior contraction rate for both in- and out-of-sample prediction error. When the number of nonzero groups times $\log(G)$ is larger than $\log(T)$, our rate is the same (up to a logarithmic factor) as the minimax rate over a class of group sparse vectors derived in Lounici et al. (2011). Then, we obtain the posterior contraction rate for the parameter of interest and, under stronger assumptions, we establish consistent posterior dimension recovery of the true model. Our asymptotic results for in-sample predictive error and parameter recovery are similar to those reported in Ning et al. (2020) for a multivariate linear regression model with group sparsity, although their assumptions and prior structure differ from ours. In-sample posterior contraction and dimension recovery for a continuous variant of the spike-and-slab prior, which is a mixture of two Laplace densities, have been also recently considered in Ročková and George (2018) and Bai et al. (in press). Concerning the posterior predictive distribution, we provide two results. The first one holds asymptotically and shows consistency of the predictive distribution. The second one establishes optimality of the predictive distribution, in the sense that if one has a prior then our posterior predictive distribution dominates any other estimator of the distribution of a new observation.

Small-sample estimation, selection, and predictive performance is assessed numerically through Monte Carlo simulations. Results show that the proposed Bayesian MIDAS penalized models present very good in-sample properties. In particular, variable selection is achieved with high probability in a very sparse setting, quite irrespective of the size of the design matrix (up to 50 high-frequency predictors in the Monte Carlo experiments) and the choice of the shape of the weighting scheme in the DGP. However, both estimation and selection performance generally deteriorate with very high cross-correlation among the original high-frequency predictors, as the Group Lasso is not designed to handle strong collinearity between regressors. Simulations also point to good out-of-sample performance, especially in comparison with alternative penalized regressions, such as Lasso, Elastic-Net, SCAD, and MC$+$. Finally, we illustrate our approach in an empirical forecasting application to U.S. GDP growth with 134 real and financial predictors sampled at monthly, weekly, and daily frequencies. We show that our models can provide superior point and density forecasts at short-term horizons (nowcasting and 1-quarter-ahead) compared to simple as well as sophisticated competing models, such as folded-concave penalties, Bayesian Model Averaging, optimally combined univariate Bayesian MIDAS models, and Factor MIDAS models.

The paper is structured as follows. Sections 2 and 3 introduce the MIDAS penalized regressions and the Bayesian MIDAS framework. Section 4 presents the theoretical analysis. In Section 5 we describe the Gibbs sampling and we discuss the Empirical Bayes approach used to automatically tune the penalty hyper-parameters. Section 6 presents simulation experiments and Section 7 provides an empirical application to the U.S. GDP growth. Finally, Section 8 concludes. The appendix contains technical proofs on the theorems reported in Section 4, while additional results and proofs of technical lemmas are collected in the Supplementary Appendix available on-line.

2 MIDAS penalized regressions

2.1 Basic MIDAS setup

Consider the variable $y_{t}$, which is observed at discrete times (i.e. only once between $t-1$ and $t$), and suppose that we want to use information stemming from a set of $K$ predictors $\mathbf{x}_{t}^{(m)}:=(x_{1,t}^{(m)},\dots,x_{K,t}^{(m)})^{\prime}$, which are observed $m$ times between $t-1$ and $t$, for forecasting purposes. The variables $y_{t}$ and $x_{k,t}^{(m)}$, for $k=1,\dots,K$, are said to be sampled at different frequencies. For instance, quarterly and monthly frequencies, respectively, in which case $m=3$. Let us define the high-frequency lag operator L^{\scalebox{0.55}{\left.1\middle/m\right.}}, such that L^{\scalebox{0.55}{\left.1\middle/m\right.}}x_{k,t}^{(m)}=x_{k,t-\scalebox{0.55}{\left.1\middle/m\right.}}^{(m)}. Further, let $h=0,1/m,2/m,3/m,\dots$ be an (arbitrary) forecast horizon, where $h=0$ denotes a nowcast with high-frequency information fully matching the low-frequency sample.

The traditional, and most simple, way of dealing with mixed-frequency data is to aggregate the high-frequency predictors by averaging (i.e. assigning equal weights), and estimate the regression through least-squares. However, this approach (flat-weighting scheme) may result in an omitted variable bias if the true weighting scheme is not flat. For a set of processes governing the high-frequency predictors, Andreou et al. (2010) show analytically and numerically that the flat least squares estimator is asymptotically inefficient (in terms of asymptotic bias and/or asymptotic variance) compared to an alternative estimator (linear or non-linear), provided for instance within the MIDAS framework, that can account for a curvature in the weighting scheme. More specifically, the MIDAS approach plugs-in the high-frequency lagged structure of predictors $x_{k,t-h}^{(m)}$ in a regression model for the low-frequency response variable $y_{t}$ as follows:

[TABLE]

where $\epsilon_{t}$ is i.i.d. with mean zero and variance $\sigma^{2}<\infty$, and $B\left(c;\boldsymbol{\theta}_{k}\right)$ is a weighting function which depends on a vector of parameters $\boldsymbol{\theta}_{k}$ and a lag order $c=0,\dots,C-1$. In this study, we consider a simple polynomial approximation of the underlying true weighting structure provided by the Almon lag polynomial $B\left(c;\boldsymbol{\theta}_{k}\right)=\sum_{i=0}^{p_{k}}\theta_{k,i}c^{i}$, where $\boldsymbol{\theta}_{k}:=(\theta_{k,0},\theta_{k,1},\dots,\theta_{k,p_{k}})^{\prime}$. Under the so-called “direct method” (Cooper, 1972), Equation (1) with Almon lag polynomials can be reparameterized as:

[TABLE]

where $\boldsymbol{\theta}:=\left(\boldsymbol{\theta}_{1}^{\prime},\dots,\boldsymbol{\theta}_{K}^{\prime}\right)^{\prime}$ is a vector featuring $\sum_{k=1}^{K}(p_{k}+1)$ parameters and $\mathbf{Z}_{t}^{(m)}:=\left(\mathbf{z}_{1,t}^{(m)\prime},\dots,\mathbf{z}_{K,t}^{(m)\prime}\right)^{\prime}$ a $\left(\sum_{k=1}^{K}(p_{k}+1)\times 1\right)$ vector of linear combinations of the observed high-frequency regressors, where each sub-vector is defined as $\mathbf{z}_{k,t}^{(m)}:=\mathbf{Q}_{k}\mathbf{x}_{k,t}^{(m)}$, with \mathbf{x}_{k,t}^{(m)}:=\left(x_{k,t}^{(m)},x_{k,t-\scalebox{0.55}{\left.1\middle/m\right.}}^{(m)},\dots,x_{k,t-\scalebox{0.55}{\left.(C-1)\middle/m\right.}}^{(m)}\right)^{\prime} a $(C\times 1)$ vector of high-frequency lags, and $\mathbf{Q}_{k}$ is a $(p_{k}+1\times C)$ polynomial weighting matrix, with $(i+1)$-th row $[0^{i},1^{i},2^{i},\dots,(C-1)^{i}]$ for $i=0,\dots,p_{k}$. The $h$-step-ahead direct forecast $\widehat{y}_{T|T-h}$, conditional on sample information available up to $T-h$, can be hence obtained using (2):

[TABLE]

for some point estimators $\widehat{\alpha}$ and $\widehat{\boldsymbol{\theta}}$ of $\alpha$ and $\boldsymbol{\theta}$, respectively. The main advantage of the Almon lag polynomial is that (2) is linear in $\boldsymbol{\theta}$ and parsimonious, as it depends only on $\sum_{k=1}^{K}(p_{k}+1)$ parameters, and can be estimated consistently and efficiently via standard methods. However, two additional features make the Almon lag polynomial particularly attractive in the present framework. First, linear restrictions on the value and slope of the lag polynomial $B\left(c;\boldsymbol{\theta}_{k}\right)$ may be placed for any $c\in(0,C-1)$. Endpoint restrictions, such as $B\left(C-1;\boldsymbol{\theta}_{k}\right)=0$ and $\nabla_{c}B\left(c;\boldsymbol{\theta}_{k}\right)|_{c=C-1}=0$, may be desirable and economically meaningful, as they jointly constrain the weighting structure to tail off slowly to zero. This can be obtained by modifying the $\mathbf{Q}_{k}$ matrix consistently with the form and the number of restrictions considered (Smith and Giles, 1976). As a result, the number of parameters in (2) reduces from $\sum_{k=1}^{K}(p_{k}+1)$ to $\sum_{k=1}^{K}(p_{k}-r_{k}+1)$, where $r_{k}\leq p_{k}$ is the number of restrictions. Second, a slope coefficient that captures the overall impact of lagged values of $x_{k,t-h}^{(m)}$ on $y_{t}$ can be easily computed as $\widehat{\beta}_{k}=\widehat{\boldsymbol{\theta}}_{k}^{\prime}\mathbf{Q}_{k}\boldsymbol{\iota}_{C}$, where $\boldsymbol{\iota}_{C}$ is a $(C\times 1)$ vector of ones.

2.2 MIDAS penalized regressions

Although appealing, the MIDAS regression presented above may be easily affected by over-parameterization and multicollinarity in presence of a large number of potentially correlated predictors.222The direct method used in regression (2) may be also hampered by multicollinearity in the artificial variables $\mathbf{Z}_{t}^{(m)}$ (Cooper, 1972). However, if $p$ is small, the imprecision arising from multicollinearity may be compensated by the lower number of coefficients to be estimated. To achieve variable selection and parameter estimation simultaneously, penalized least squares procedures such as the Lasso (Tibshirani, 1996) have been recently investigated in a high-dimensional mixed-frequency framework by Uematsu and Tanaka (2019) and Lima et al. (2020). Cognizant of the limits of the Lasso (estimation bias and lack of selection consistency), the Adaptive Lasso (Zou, 2006; AL hereafter), which enjoys selection consistency under typically weaker assumptions than the so-called irrepresentable condition on the design matrix, represents a tempting alternative to the Lasso that has not yet been explored in this literature.333The irrepresentable condition states that the predictors not in the model are not representable by predictors in the true model (i.e. the irrelevant predictors are roughly orthogonal to the relevant ones). This represents a necessary and sufficient condition for the Lasso for exact recovery of the non-zero coefficients, but it can be easily violated in cases where the design matrix exhibits too strong (empirical) correlations (collinearity between predictors).

However, in this paper we argue that the AL might still not be suited in a mixed-frequency framework. The rationale is that since the distributed lags in the MIDAS regression (2) are by construction highly correlated, the AL estimator would tend to randomly select only one term of each lag polynomial in the active set and shrink the remaining (relevant) coefficients to zero. The theoretical rationale for a failure in the selection ability of the AL in our mixed-frequency setting is hence similar to the one pointed out by Efron et al. (2004) and Zou and Hastie (2005), and it is mostly related to the lack of strict convexity in the Lasso penalty. To address this issue, we propose a solution based on the Adaptive Group Lasso estimator (Yuan and Lin, 2006; Wang and Leng, 2008; AGL hereafter). This approach introduces a penalty to a group of regressors, rather than a single regressor, that may lead (if the group structure is carefully set by the researcher) to a finite sample improvement of the AL. In the present framework, it seems reasonable to define a group as each of the $k$ vectors of lag polynomials in the model. This grouping structure is motivated by the fact that if one high-frequency predictor is irrelevant, it is also expected that zero-coefficients occur in all the parameters of its lag polynomial. Hence, this strategy should overcome, at least in part, the limitations of both Lasso and AL in presence of strong correlation in the design matrix arising from the correlation among lags of the transformed high-frequency predictors $\mathbf{Z}_{t}^{(m)}$.

Let us now assume that $y_{t}$ is centered at [math] and regressors $\mathbf{Z}_{t}^{(m)}$ in (2) are standardized. Let us also consider, for ease of exposition, that the order of the Almon lag polynomial and the number of linear restrictions are the same across variables, such that $p_{k}=p$ and $r_{k}=r$, for $k=1,\dots,K$, and the total number of parameters is $K(p-r+1)$. Further, partition the parameter vector $\boldsymbol{\theta}$ into $G$ disjoint groups, $\boldsymbol{\theta}_{j}$, for $j=1,\dots,G$, each of size $g_{j}$ and including the lag polynomial of a given predictor. Despite the change in notation (necessary to avoid confusion), it is straightforward that $G=K$, $\boldsymbol{\theta}_{j}=\boldsymbol{\theta}_{k}$, $g_{j}=p-r+1$, and $\widetilde{g}:=\sum_{j=1}^{G}g_{j}=K(p-r+1)$. The objective function of the AGL takes the form:

[TABLE]

where $\mathcal{L}_{T}(\boldsymbol{\theta})$ denotes the negative log-likelihood function, $\|\boldsymbol{\theta}_{j}\|_{2}=(\boldsymbol{\theta}_{j}^{\prime}\boldsymbol{\theta}_{j})^{1/2}$ the $\ell_{2}$-norm, and $\lambda_{j}$ the group penalty parameter. Estimation and selection consistency of the AGL estimator are established by Wang and Leng (2008). However, as suggested by Callot and Kock (2014), the AGL possesses a variant of the oracle property only if the predictors are partitioned in the correct groups. This happens because selection consistency concerns the inactive groups (i.e. including only parameters whose true value is zero), while the asymptotic distribution of those parameters whose true value is zero but are located in active groups is equivalent to the one of least squares including all variables. Hence, the AGL only performs better than least squares including all variables if one is able to identify the true inactive groups. In the present framework, we expect that partitioning each lag polynomial into one group should attenuate this issue.

3 Bayesian MIDAS penalized regressions

Since the work of Yuan and Lin (2006), several algorithms based on convex optimization methods (e.g. block-wise and proximal gradient descent) have been proposed to compute the solution paths of the Adaptive Group Lasso (Wang and Leng, 2008). In this paper, we consider instead a Bayesian hierarchical approach based on MCMC methods (Kyung et al., 2010), which has several advantages compared to the frequentist approach. First, Bayesian methods exploit model inference via posterior distributions of parameters, which usually provide a valid measure of standard errors based on a geometrically ergodic Markov chain (Khare and Hobert, 2013).444It is nevertheless worth noting that the results in Khare and Hobert (2013) hold as long as the penalty hyper-parameters are assumed fixed, while convergence properties of the MCMC algorithm for the full Bayesian penalized regression models are still unknown (see also Roy and Chakraborty, 2017). Second, they provide a flexible way of estimating the penalty parameters, along with other parameters in the model. Lastly, they provide forecasts via predictive distributions. In what follows, we present in detail the hierarchical structure of the proposed Bayesian MIDAS penalized models.

3.1 Bayesian MIDAS Adaptive Group Lasso

As noted by Tibshirani (1996), the Lasso estimator can be interpreted as the Bayes posterior mode using normal likelihood and independent Laplace (double-exponential) prior for the regression coefficients. Park and Casella (2008) propose a Bayesian Lasso where the $\ell_{1}$ penalty corresponds to a conditional Laplace prior that can be represented as a scale mixture of Normals with an exponential mixing density. For the Bayesian Adaptive Group Lasso, we consider a multivariate generalization of the double exponential prior as in Kyung et al. (2010) and Leng et al. (2014), such that the conditional prior of $\boldsymbol{\theta}$ given $\sigma^{2}$ can be expressed as a scale mixture of Normals with Gamma hyper-priors. The conditional prior of $\boldsymbol{\theta}_{j}$, with $j=1,\dots,G$, is then (see Appendix A.1):

[TABLE]

This suggests the following hierarchical representation of the Bayesian MIDAS Adaptive Group Lasso model (BMIDAS-AGL):

[TABLE]

where $y:=(y_{1},\dots,y_{T})^{\prime}$ is centered at [math] and $\mathbf{Z}:=\left(\mathbf{Z}_{1}^{(m)},\dots,\mathbf{Z}_{T}^{(m)}\right)^{\prime}$ is the $T\times\widetilde{g}$ matrix of standardized regressors, and $\mathbf{I}_{g_{j}}$ is an identity matrix of order $g_{j}$, for $j=1,\dots,G$. By integrating out $\tau_{j}^{2}$ in the hierarchy above, we obtain that the marginal prior for $\boldsymbol{\theta}_{j}$, given $\sigma^{2}$, follows a $g_{j}$-dimensional Multi-Laplace distribution with density function $\text{M-Laplace}(\boldsymbol{\theta}_{j}|\mathbf{0},\sigma/\lambda_{j})\propto(\lambda_{j}/\sigma)^{g_{j}}\exp(-\lambda_{j}\|\boldsymbol{\theta}_{j}\|_{2}/\sigma)$ (see Appendix A.1). Let $\boldsymbol{\tau}:=(\tau_{1}^{2},\dots,\tau_{G}^{2})$ and $\boldsymbol{\lambda}:=(\lambda_{1}^{2},\dots,\lambda_{G}^{2})$. The full posterior distribution of the unknown parameters conditional on the data and for some penalty hyper-parameters $\boldsymbol{\lambda}$ is:

[TABLE]

3.2 Bayesian MIDAS Adaptive Group Lasso with Spike-and-Slab Prior

The Bayesian model outlined above would typically induce sparsity by shrinking the coefficients of the inactive set towards zero, but usually not exactly to zero. To achieve exact sparse solutions, recent literature has increasingly focused on introducing probabilistic sparse recovery by adding a point mass mixture prior to penalized regressions (Zhang et al., 2014; Zhao and Sarkar, 2015; Ročková and George, 2018). In the present study, we follow Xu and Ghosh (2015) and we consider a Bayesian Group Lasso with spike-and-slab prior for group variable selection, which provides two shrinkage effects: a point mass at $\mathbf{0}$ leading to exact zero coefficients (the spike part of the prior) and a Group Lasso prior on the slab part. The combination of these two components together is expected to facilitate variable selection at the group level and to shrink coefficients in the selected groups simultaneously. Similarly to the BMIDAS-AGL, the hierarchical Bayesian MIDAS Adaptive Group Lasso with spike-and-slab prior (BMIDAS-AGL-SS) is:

[TABLE]

where $\delta_{0}(\boldsymbol{\theta}_{j})$ denotes a point mass at $\mathbf{0}\in\mathbb{R}^{g_{j}}$, for $j=1,\dots,G$. Note that our spike-and-slab prior conditional on $\sigma^{2}$ is the same as in Ning et al. (2020), except that i) they specify an independent prior for $\boldsymbol{\theta}$ and the model variance parameter, and they have the same penalty hyper-parameter across the groups, ii) their prior for the dimension is not necessarily a Binomial distribution as in our case, and iii) we place a conjugate Beta prior on $\pi_{0}$ while they consider $\pi_{0}$ as fixed. Adding a Beta prior to the hierarchy allows for mixing over the sparsity level. We use typical non-informative priors for the error variance $\sigma^{2}$ and we set the hyper-parameters $c$ and $d$ according to the theoretical results presented is Section 4 (see also Section 5.1 for the exact values used in the present paper). By integrating out $\tau_{j}^{2}$, we now obtain that the marginal prior for $\boldsymbol{\theta}_{j}$ is a mixture of point mass at $\mathbf{0}\in\mathbb{R}^{g_{j}}$ and a $g_{j}$-dimensional Multi-Laplace distribution (see Appendix A.1). The full posterior distribution of the unknown parameters conditional on the data and for some penalty hyper-parameters $\boldsymbol{\lambda}$ is:

[TABLE]

4 Theoretical properties

This section provides the theoretical validation of our Bayesian MIDAS procedure. In Section 4.1 we investigate the asymptotic properties of our procedure by adopting a frequentist point of view. In Section 4.2 we address the question of optimality of the corresponding posterior predictive density. Further results on dimension recovery of the model and distributional approximation of the posterior of $\boldsymbol{\theta}$ are provided in the Supplementary Appendix B.2. As we assume sparsity in the data generating process, in the asymptotic analysis we focus on the BMIDAS-AGL-SS model, which can provide exact sparse recovery of the true model. We use the notation $\Pi(\cdot|y,\mathbf{Z})$ (resp. $\Pi(\cdot)$) to denote both the posterior (resp. prior) distribution and its Lebesgue density function.

4.1 Asymptotic analysis

In this section we establish the posterior contraction rate for both the in-sample and the out-of-sample prediction error, as well as the consistency of both the marginal posterior of $\boldsymbol{\theta}$ and the posterior predictive density of a new observation. We adopt a frequentist point of view, in the sense that we admit the existence of a true value of $(\boldsymbol{\theta},\sigma^{2})$, denoted by $(\boldsymbol{\theta}_{0},\sigma_{0}^{2})$, that generates the data. We denote by $\mathbf{E}_{0}[\cdot]$ the expectation taken with respect to the true data distribution $\mathcal{N}(\mathbf{Z}\boldsymbol{\theta}_{0},\sigma_{0}^{2}\mathbf{I}_{T})$ conditional on $(\mathbf{Z},\boldsymbol{\theta}_{0},\sigma_{0}^{2})$.

Recall that $g_{j}$ is the size of the $j$-th group and $\widetilde{g}:=\sum_{j=1}^{G}g_{j}$ is the total number of slope parameters in (2). Let $g_{\max}:=\max_{1\leq j\leq G}g_{j}$ be the maximum group size and let $S_{0}\subseteq\{1,\ldots,G\}$ denote the set containing the indices of the true nonzero groups with $s_{0}:=|S_{0}|$ its cardinality. For a generic vector $\mathbf{v}$, denote by $s_{\mathbf{v}}$ the number of nonzero groups of $\mathbf{v}$. Moreover, let $\|\boldsymbol{\theta}\|_{\infty}:=\max_{1\leq j\leq\widetilde{g}}|\theta_{j}|$, $\|\mathbf{Z}\|_{op}$ be the spectral norm and $\|\mathbf{Z}\|_{o}:=\max\{\|\mathbf{Z}_{j}\|_{op};1\leq j\leq G\}$ where $\mathbf{Z}_{j}$ is the $(T\times g_{j})$-submatrix of $\mathbf{Z}$ made of all the rows and the columns corresponding to the indices in the $j$-th block. We introduce the following assumptions, which are quite mild in that they put weak restrictions on the prior hyper-parameters and allow $g_{\max}$ and the largest element of $\boldsymbol{\theta}_{0}$ to increase with $T\rightarrow\infty$.

Assumption 4.1

For positive and bounded constants $\underline{\sigma}^{2}$, $\overline{\sigma}^{2}$, $c_{g}$ and $c_{\boldsymbol{\theta}}$, suppose that: (i) $0<\underline{\sigma}^{2}\leq\sigma_{0}^{2}\leq\overline{\sigma}^{2}<\infty$; (ii) $g_{\max}\leq c_{g}\log(G)/\log(T)$; (iii) $\|\boldsymbol{\theta}_{0}\|_{\infty}\leq c_{\boldsymbol{\theta}}\sqrt{\log(G)}$.

Define the region $\overline{\Theta}_{0}:=\{\boldsymbol{\theta}\in\mathbb{R}^{\widetilde{g}};\,\|\boldsymbol{\theta}\|_{\infty}\leq c_{\boldsymbol{\theta}}\sqrt{\log(G)}\}$, where $c_{\boldsymbol{\theta}}$ is as defined in Assumption 4.1. By Assumption 4.1 (iii), $\boldsymbol{\theta}_{0}\in\overline{\Theta}_{0}$.

Assumption 4.2

Let $\lambda_{\min}:=\min\{\lambda_{j};\,j\leq G\}$ and $\lambda_{\max}:=\max\{\lambda_{j};\,j\leq G\}$ and assume that $\|\mathbf{Z}\|_{0}/\sqrt{T}<\sqrt{\log(T)}$ with probability $1$. The scale parameters $\lambda_{j}$ are allowed to change with $T$ and belong to the range:

[TABLE]

for two positive constants $1<\underline{c}<\overline{C}<\infty$.

Assumption 4.3

The hyper-parameters $c,d$ of the Beta prior for $\pi_{0}$ satisfy: $(d+j-1)/(c+G-j)\leq j/[G^{u}(G-j+1)]$ for every $0<u<s_{0}$ and $\forall j\in\{1,\ldots,G\}\subseteq\mathbb{N}$, $d\leq s_{0}$ and $\log(G+c)\leq c_{\pi_{0}}\log(G)$ for some fixed constant $0<c_{\pi_{0}}<\infty$. Moreover, the hyperparameter $a_{1}$ of the prior for $\sigma^{2}$ satisfies $a_{1}>1$.

Assumption 4.1 (ii) restricts the rate at which the maximum group size increases to infinity. Assumption 4.1 (iii) restricts the growth rate of the maximum component of the true $\boldsymbol{\theta}_{0}$. Assumption 4.1 (i) excludes degenerate cases by restricting the model variance. In the following we establish results that hold uniformly over $(\boldsymbol{\theta}_{0},\sigma_{0}^{2})\in\overline{\Theta}_{0}\times[\underline{\sigma}^{2},\overline{\sigma}^{2}]$. Assumption 4.2 restricts the limiting behaviour of the parameters $\boldsymbol{\lambda}$ of the M-Laplace distribution. The upper bound prevents to shrink the non-zero groups too much towards zero while the lower bound prevents many false signals. Assumption 4.3 is a mild assumption on the prior on $\pi_{0}$, it allows $c$ of the form $c=\bar{\kappa}G^{v}$, for constants $\bar{\kappa}>1$ and $v>1$ as in Castillo et al. (2015, Example 5) but it excludes a $v$ that increases with $G$.

4.1.1 In-sample asymptotic properties

We now investigate the in-sample asymptotic properties of our model $y|\mathbf{Z},\mathbf{\boldsymbol{\theta}},\sigma^{2}\sim\mathcal{N}_{T}(\mathbf{Z}\mathbf{\boldsymbol{\theta}},\sigma^{2}\mathbf{I}_{T})$, where the parameter $(\boldsymbol{\theta}^{\prime},\sigma^{2})^{\prime}$ is endowed with the BMIDAS-AGL-SS prior (A2). The first theorem establishes posterior consistency of the in-sample prediction error.

Theorem 4.1

Suppose Assumptions 4.1 - 4.3 hold, and let $\epsilon:=\max\{\sqrt{s_{0}\log(G)/T},\sqrt{\log(T)/T}\}\allowbreak\rightarrow 0$. Then, for a sufficiently large $M>0$:

[TABLE]

The assumption $\epsilon\rightarrow 0$ in Theorem 4.1 implies that $\log(G)=o(T)$. The latter becomes $\log(dimension(\boldsymbol{\theta}))=o(T)$ when the number of groups is equal to the number of parameters, which has been shown in the literature to be a necessary condition for sparse recovery (see e.g. Lounici et al., 2011). The contraction rate $\epsilon$ is the same as the one in Ning et al. (2020) under our assumptions. The first term of the rate coincides (up to a logarithmic factor) with the minimax rate over a class of group sparse vectors derived in Lounici et al. (2011). As discussed in that paper, the Group Lasso estimator has some advantages over the Lasso estimator in some important cases. One of these cases is the one we consider where $g_{\max}<\log(G)$ under our Assumption 4.1 (ii). In this case the upper bound $\epsilon$ of the posterior contraction rate for the in-sample prediction error is faster than the lower bound on the prediction of the Lasso estimator (see Lounici et al., 2011, Section 7). Finally, when the number of groups is equal to the number of parameters, then the first term of the rate corresponds to the optimal one for this problem (see e.g. Bühlmann and Van De Geer, 2011, and Castillo et al., 2015).

The proof of Theorem 4.1 is provided in Appendix A.2 and follows the strategy of Ning et al. (2020) adapted to our slightly different prior and model. Instead of relying on the general posterior contraction theory based on the Hellinger distance, which is not appropriate in this setting, this strategy of proof relies on the average Rényi divergence, which can be easily used to obtain the rate in terms of Euclidean distance. To obtain the posterior contraction rate with respect to the Rényi divergence, we use the general posterior contraction procedure for independent observations, as in Ghosal and van der Vaart (2007). That is, we first show that our prior (A2) puts enough mass in a shrinking neighborhood of the true density function. Then, we show that the prior puts a small mass, decreasing to zero, on the complement of a finite dimensional space with increasing dimension that approximates well the infinite dimensional model. Finally, we prove the existence of exponential tests of the true density against the complement of balls around the truth. These results are proved in Lemmas B.1.2 and B.1.5 in the Supplementary Appendix.

Given the result in Theorem 4.1, we are now ready to establish a result on parameter recovery of our procedure, that is, consistency of the marginal posterior of $\boldsymbol{\theta}$. As discussed in Bickel et al. (2009), the parameter $\boldsymbol{\theta}$ is not estimable without a condition on $\mathbf{Z}$ because of its large dimension. In particular, it is well known from the literature (see e.g. Castillo et al., 2015, and Ning et al., 2020) that if $\boldsymbol{\theta}$ is sparse, then local invertibility of the Gram matrix $\mathbf{Z}^{\prime}\mathbf{Z}$ is sufficient for estimability of $\boldsymbol{\theta}$. We then introduce the following quantity.

Definition 1

For every $s>0$, the smallest scaled sparse singular value of dimension $s$ is defined as

[TABLE]

The restricted eigenvalue condition requires that $\widetilde{\phi}(s)>0$. This means that the smallest eigenvalue for the sub-matrix of $\mathbf{Z}$ made of columns corresponding to the non-zero groups is strictly larger than zero. We call it “scaled” eigenvalue because we divide by the maximum operator norm of submatrices of $\mathbf{Z}$. We point out that since $\widetilde{g}$ is large, possibly larger than $T$, then it would be unrealistic to require that the minimum eigenvalue of $\mathbf{Z}^{\prime}\mathbf{Z}$ is strictly positive. Our restricted sparse eigenvalue condition is instead much weaker than this condition, and it is hence more realistic.

Let $\widetilde{s}_{0}:=\max\{s_{0},\log(T)/\log(G)\}$. Lemma B.1.4 in the Supplementary Appendix shows that $\mathbf{E}_{0}[\Pi(\boldsymbol{\theta};s_{\boldsymbol{\theta}}\geq M_{2}\widetilde{s}_{0}|y,\mathbf{Z})]\rightarrow 0$ uniformly over $(\boldsymbol{\theta}_{0},\sigma_{0}^{2})\in\overline{\Theta}_{0}\times[\underline{\sigma}^{2},\overline{\sigma}^{2}]$. This means that with probability approaching $1$, the posterior puts mass one on the $\boldsymbol{\theta}$’s with dimension $\widetilde{s}_{0}$. This implies that when $\log(T)/\log(G)>s_{0}$, the posterior can go beyond the true dimension $s_{0}$. In spite of this feature, the next theorem shows that the posterior is still able to recover the true parameter $\boldsymbol{\theta}_{0}$, similarly to Ning et al. (2020, Corollary 3.3).

Theorem 4.2

Suppose Assumptions 4.1 - 4.3 hold, and let $\epsilon:=\max\{\sqrt{s_{0}\log(G)/T},\sqrt{\log(T)/T}\}\allowbreak\rightarrow 0$. Then, for a sufficiently large $M>0$ and for $\widetilde{\phi}(s)$ as defined in (9) we have:

[TABLE]

The proof of Theorem 4.2 is provided in Appendix A.2. The first term of the rate coincides (up to a logarithmic factor) with the minimax rate over a class of group sparse vectors derived in Lounici et al. (2011). The posterior contraction rate provided in Theorem 4.2 deteriorates when $\widetilde{\phi}(s_{0}+M_{2}\widetilde{s}_{0})$ is small.

4.1.2 Out-of-sample asymptotic properties

We now analyse the behaviour of the out-of-sample prediction error, as well as the behaviour of the predictive density associated with our BMIDAS-AGL-SS model. For this, consider the posterior predictive density

[TABLE]

where $\mathbf{Z}_{\tau-h}$ is a $(\widetilde{g}\times 1)$ vector, $f_{0}(\cdot|\mathbf{Z}_{\tau-h},\boldsymbol{\theta},\sigma^{2})$ is the Lebesgue density of the one-dimensional $\mathcal{N}(\mathbf{Z}_{\tau-h}^{\prime}\boldsymbol{\theta},\sigma^{2})$ distribution, and denote by $\|\cdot\|_{TV}$ the Total Variation distance. In addition, we denote by $f_{\mathbf{Z}_{\tau-h},\boldsymbol{\theta}_{0},\sigma_{0}^{2}}(\cdot):=f_{0}(\cdot|\mathbf{Z}_{\tau-h},\boldsymbol{\theta}_{0},\sigma_{0}^{2})$ the Lebesgue density of the one-dimensional $\mathcal{N}(\mathbf{Z}_{\tau-h}^{\prime}\boldsymbol{\theta}_{0},\sigma_{0}^{2})$ distribution. As before, the expectation $\mathbf{E}_{0}[\cdot]$ denotes the expectation taken with respect to the true data distribution of $y|\mathbf{Z},\boldsymbol{\theta}_{0},\sigma_{0}^{2}\sim\mathcal{N}(\mathbf{Z}\boldsymbol{\theta}_{0},\sigma_{0}^{2}\mathbf{I}_{T})$ conditional on $(\mathbf{Z},\boldsymbol{\theta}_{0},\sigma_{0}^{2})$. In addition, $\mathbf{E}_{\mathbf{Z}_{\tau-h}}$ has to be understood as the expectation taken with respect to the distribution of $\mathbf{Z}_{\tau-h}$. Finally, we denote by $\mathbf{Z}_{\tau-h,j}$ the $(g_{j}\times 1)$-subvector of $\mathbf{Z}_{\tau-h}$ made of all the elements corresponding to the indices in the $j$-th block of $\mathbf{Z}_{\tau-h}$.

As before, we assume that $(\boldsymbol{\theta}^{\prime},\sigma^{2})^{\prime}$ is endowed with the BMIDAS-AGL-SS prior (A2). We start by establishing the contraction rate of the posterior of the out-of-sample prediction error.

Theorem 4.3

Suppose the assumptions of Theorem 4.2 be satisfied. For every $\tau>T$ and every fixed $\mathbf{Z}_{\tau-h}$ define:

[TABLE]

for a constant $M_{2}>\frac{2+C_{1}}{u}+1$ that does not depend on $(\boldsymbol{\theta}_{0},\sigma_{0}^{2})$ and $C_{1}>0$ the constant in Lemma B.1.1. Moreover, let $\epsilon:=\max\{\sqrt{s_{0}\log(G)/T},\sqrt{\log(T)/T}\}\rightarrow 0$. Then, $\forall\tau>T$:

[TABLE]

The proof of Theorem 4.3 is provided in Appendix A.2. Result (11) in the previous theorem provides a contraction rate for the out-of-sample prediction error for a given $\mathbf{Z}_{\tau-h}$. This rate might be slower than the rate for the in-sample prediction error if $\eta_{0}/\widetilde{\phi}(s_{0}+M_{2}\widetilde{s}_{0})$ is not bounded in probability. To get this rate we have exploited the sparsity structure and the property of the posterior distribution to concentrate on vectors $\boldsymbol{\theta}$ with dimension $\widetilde{s}_{0}$. On the other hand, if we assumed that $\mathbf{E}_{\mathbf{Z}_{\tau-h}}[\mathbf{Z}_{\tau-h}^{4}]<\kappa$ for some constant $\kappa$, as in Carrasco and Rossi (2016), which in turn implies $\|\mathbf{Z}_{\tau-h}\|_{2}=O_{p}(1)$, then we would get the same rate as in Theorem 4.2. Here we do not impose this restriction. It is worth noting that the numerator of $\eta_{0}$ can be upper bounded by $\max_{1\leq j\leq G}|\mathbf{Z}_{\tau-h,j}^{\prime}\boldsymbol{\theta}_{j}|^{2}(M_{2}\widetilde{s}_{0}+s_{0})^{2}$, which diverges. However, $\eta_{0}$ does not necessarily diverge if its denominator diverges at least at the same rate.

We now state the consistency of the posterior predictive density function.

Theorem 4.4

Suppose the assumptions of Theorem 4.2 be satisfied. Let $\epsilon:=\max\{\sqrt{s_{0}\log(G)/T},\allowbreak\sqrt{\log(T)/T}\}\rightarrow 0$. Then, $\forall\tau>T$ such that $\mathbf{E}_{\mathbf{Z}_{\tau-h}}[\mathbf{Z}_{\tau-h,j}\mathbf{Z}_{\tau-h,k}^{\prime}]=0$ for every $j\neq k$ and $\epsilon^{2}\frac{\left\|\mathbf{E}_{\mathbf{Z}_{\tau-h}}[\mathbf{Z}_{\tau-h}\mathbf{Z}_{\tau-h}^{\prime}]\right\|_{o}}{\|\mathbf{Z}/\sqrt{T}\|_{o}^{2}\widetilde{\phi}(s_{0}+M_{2}\widetilde{s}_{0})}=o_{p}(1)$ the following holds:

[TABLE]

The theorem establishes consistency of the posterior predictive distribution under the condition that the out-of-sample predictors are uncorrelated across groups. This assumption seems reasonable given the group structure. However, if this condition is not verified then the consistency result remains true either under the condition that $\epsilon^{2}\frac{\left\|\mathbf{E}_{\mathbf{Z}_{\tau-h}}[\mathbf{Z}_{\tau-h}\mathbf{Z}_{\tau-h}^{\prime}]\right\|_{op}}{\|\mathbf{Z}/\sqrt{T}\|_{o}^{2}\widetilde{\phi}(s_{0}+M_{2}\widetilde{s}_{0})}=o_{p}(1)$, which is a little bit stronger, or under the condition $\epsilon^{2}\sup_{\{\boldsymbol{\theta};0\leq s_{\boldsymbol{\theta}}\leq M_{2}\widetilde{s}_{0}+s_{0}\}}\frac{\mathbf{E}_{\mathbf{Z}_{\tau-h}}|\mathbf{Z}_{\tau-h}^{\prime}\boldsymbol{\theta}|^{2}}{\sum_{j=1}^{G}\|\mathbf{Z}_{j}^{\prime}\boldsymbol{\theta}_{j}/\sqrt{T}\|_{2}^{2}\widetilde{\phi}(s_{0}+M_{2}\widetilde{s}_{0})}=o_{p}(1)$, which involves a quantity similar to the rate in Theorem 4.3.

4.2 Optimality of the predictive density

While Theorem 4.4 states good asymptotic properties of the posterior predictive density function, it is also possible to argue that the posterior predictive density function has optimal properties in finite sample. Indeed, if one has two contenders for estimating $f_{\mathbf{Z}_{\tau-h},\boldsymbol{\theta}_{0},\sigma_{0}^{2}}(y_{\tau})$, say our predictive density $f_{\mathbf{Z}_{\tau-h},y,\mathbf{Z}}(y_{\tau})$ and another density $\widetilde{f}_{\mathbf{Z}_{\tau-h},y,\mathbf{Z}}(y_{\tau})$, hence $f_{\mathbf{Z}_{\tau-h},y,\mathbf{Z}}$ is better than $\widetilde{f}_{\mathbf{Z}_{\tau-h},y,\mathbf{Z}}$ if the following expected difference of the Kullback-Leibler divergences with respect to the true $f_{\mathbf{Z}_{\tau-h},\boldsymbol{\theta}_{0},\sigma_{0}^{2}}$ is positive:

[TABLE]

where $K(f_{1}\|f_{2}):=\int\log(f_{1}/f_{2})f_{1}$ is the Kullback-Leibler divergence between two probability measures with Lebesgue densities $f_{1}$ and $f_{2}$, $f_{\mathbf{Z},\boldsymbol{\theta}_{0},\sigma_{0}^{2}}(y)$ is the Lebesgue density of the $\mathcal{N}(\mathbf{Z}\boldsymbol{\theta}_{0},\sigma_{0}^{2}\mathbf{I}_{T})$ distribution, and $f_{\mathbf{Z}}(y)=\int f_{\mathbf{Z},\boldsymbol{\theta}_{0},\sigma_{0}^{2}}(y)\Pi(\boldsymbol{\theta}_{0},\sigma_{0}^{2}|y,\mathbf{Z})d\boldsymbol{\theta}_{0}d\sigma_{0}^{2}$. Since the right hand side of (13) is equal to $\int f_{\mathbf{Z}}(y)dy\,K(f_{\mathbf{Z}_{\tau-h},y,\mathbf{Z}}\|\widetilde{f}_{\mathbf{Z}_{\tau-h},y,\mathbf{Z}})$, then it is strictly positive. This shows that if one has a prior, and exploits it, then our posterior predictive density dominates every other density estimator. In particular, it dominates any density estimator that replaces $(\boldsymbol{\theta}^{\prime},\sigma^{2})^{\prime}$ by a consistent estimator. Of course, this comparison does not hold when no prior information is available. However, for some distributions, it is known that the Bayesian predictive density outperforms other density estimators even without a prior (see e.g. Aitchison, 1975).

5 Gibbs sampler

5.1 Conditional posteriors

We now provide details on the Gibbs sampler we use to simulate from the posterior distributions (3.1) and (3.2). For this purpose, we consider an efficient block Gibbs sampler (Hobert and Geyer, 1998). Let us denote $\boldsymbol{\theta}_{\mathbin{\mathchoice{\hbox{ \leavevmode\hbox to3.6pt{\vbox to6.6pt{\pgfpicture\makeatletter\hbox{\hskip 0.3pt\lower-0.3pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{} {{}{}}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.6pt}\pgfsys@invoke{ }\pgfsys@roundcap\pgfsys@invoke{ }{}\pgfsys@moveto{3.0pt}{0.0pt}\pgfsys@lineto{0.0pt}{6.0pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}{\hbox{ \leavevmode\hbox to3.6pt{\vbox to6.6pt{\pgfpicture\makeatletter\hbox{\hskip 0.3pt\lower-0.3pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{} {{}{}}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.6pt}\pgfsys@invoke{ }\pgfsys@roundcap\pgfsys@invoke{ }{}\pgfsys@moveto{3.0pt}{0.0pt}\pgfsys@lineto{0.0pt}{6.0pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}{\hbox{ \leavevmode\hbox to2.45pt{\vbox to4.45pt{\pgfpicture\makeatletter\hbox{\hskip 0.22499pt\lower-0.22499pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{} {{}{}}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.45pt}\pgfsys@invoke{ }\pgfsys@roundcap\pgfsys@invoke{ }{}\pgfsys@moveto{2.0pt}{0.0pt}\pgfsys@lineto{0.0pt}{4.0pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}{\hbox{ \leavevmode\hbox to1.9pt{\vbox to3.4pt{\pgfpicture\makeatletter\hbox{\hskip 0.2pt\lower-0.2pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{} {{}{}}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\pgfsys@roundcap\pgfsys@invoke{ }{}\pgfsys@moveto{1.5pt}{0.0pt}\pgfsys@lineto{0.0pt}{3.0pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}}j}:=(\boldsymbol{\theta}_{1}^{\prime},\dots,\boldsymbol{\theta}_{j-1}^{\prime},\boldsymbol{\theta}_{j+1}^{\prime},\dots,\boldsymbol{\theta}_{G}^{\prime})^{\prime}$ the $\boldsymbol{\theta}$ vector without the $j$th high-frequency lag polynomial, and $\mathbf{Z}_{j}$ and denote $\mathbf{Z}_{\mathbin{\mathchoice{\hbox{ \leavevmode\hbox to3.6pt{\vbox to6.6pt{\pgfpicture\makeatletter\hbox{\hskip 0.3pt\lower-0.3pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{} {{}{}}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.6pt}\pgfsys@invoke{ }\pgfsys@roundcap\pgfsys@invoke{ }{}\pgfsys@moveto{3.0pt}{0.0pt}\pgfsys@lineto{0.0pt}{6.0pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}{\hbox{ \leavevmode\hbox to3.6pt{\vbox to6.6pt{\pgfpicture\makeatletter\hbox{\hskip 0.3pt\lower-0.3pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{} {{}{}}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.6pt}\pgfsys@invoke{ }\pgfsys@roundcap\pgfsys@invoke{ }{}\pgfsys@moveto{3.0pt}{0.0pt}\pgfsys@lineto{0.0pt}{6.0pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}{\hbox{ \leavevmode\hbox to2.45pt{\vbox to4.45pt{\pgfpicture\makeatletter\hbox{\hskip 0.22499pt\lower-0.22499pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{} {{}{}}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.45pt}\pgfsys@invoke{ }\pgfsys@roundcap\pgfsys@invoke{ }{}\pgfsys@moveto{2.0pt}{0.0pt}\pgfsys@lineto{0.0pt}{4.0pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}{\hbox{ \leavevmode\hbox to1.9pt{\vbox to3.4pt{\pgfpicture\makeatletter\hbox{\hskip 0.2pt\lower-0.2pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{} {{}{}}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\pgfsys@roundcap\pgfsys@invoke{ }{}\pgfsys@moveto{1.5pt}{0.0pt}\pgfsys@lineto{0.0pt}{3.0pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}}j}$ partitions of the design matrix corresponding to $\boldsymbol{\theta}_{j}$ and $\boldsymbol{\theta}_{\mathbin{\mathchoice{\hbox{ \leavevmode\hbox to3.6pt{\vbox to6.6pt{\pgfpicture\makeatletter\hbox{\hskip 0.3pt\lower-0.3pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{} {{}{}}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.6pt}\pgfsys@invoke{ }\pgfsys@roundcap\pgfsys@invoke{ }{}\pgfsys@moveto{3.0pt}{0.0pt}\pgfsys@lineto{0.0pt}{6.0pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}{\hbox{ \leavevmode\hbox to3.6pt{\vbox to6.6pt{\pgfpicture\makeatletter\hbox{\hskip 0.3pt\lower-0.3pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{} {{}{}}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.6pt}\pgfsys@invoke{ }\pgfsys@roundcap\pgfsys@invoke{ }{}\pgfsys@moveto{3.0pt}{0.0pt}\pgfsys@lineto{0.0pt}{6.0pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}{\hbox{ \leavevmode\hbox to2.45pt{\vbox to4.45pt{\pgfpicture\makeatletter\hbox{\hskip 0.22499pt\lower-0.22499pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{} {{}{}}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.45pt}\pgfsys@invoke{ }\pgfsys@roundcap\pgfsys@invoke{ }{}\pgfsys@moveto{2.0pt}{0.0pt}\pgfsys@lineto{0.0pt}{4.0pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}{\hbox{ \leavevmode\hbox to1.9pt{\vbox to3.4pt{\pgfpicture\makeatletter\hbox{\hskip 0.2pt\lower-0.2pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{} {{}{}}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\pgfsys@roundcap\pgfsys@invoke{ }{}\pgfsys@moveto{1.5pt}{0.0pt}\pgfsys@lineto{0.0pt}{3.0pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}}j}$, respectively. With a conjugate Gamma prior placed on the penalty hyper-parameters, $\lambda^{2}_{j}\sim\textrm{Gamma}\left(a_{2},b_{2}\right)$, the full conditional posteriors for the BMIDAS-AGL model are:

[TABLE]

where $\mathbf{A}_{j}:=\mathbf{Z}_{j}^{\prime}\mathbf{Z}_{j}+\tau_{j}^{-2}\mathbf{I}_{g_{j}}$ and $\mathbf{C}_{j}:=\mathbf{Z}_{j}^{\prime}\left(y-\mathbf{Z}_{\mathbin{\mathchoice{\hbox{ \leavevmode\hbox to3.6pt{\vbox to6.6pt{\pgfpicture\makeatletter\hbox{\hskip 0.3pt\lower-0.3pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{} {{}{}}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.6pt}\pgfsys@invoke{ }\pgfsys@roundcap\pgfsys@invoke{ }{}\pgfsys@moveto{3.0pt}{0.0pt}\pgfsys@lineto{0.0pt}{6.0pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}{\hbox{ \leavevmode\hbox to3.6pt{\vbox to6.6pt{\pgfpicture\makeatletter\hbox{\hskip 0.3pt\lower-0.3pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{} {{}{}}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.6pt}\pgfsys@invoke{ }\pgfsys@roundcap\pgfsys@invoke{ }{}\pgfsys@moveto{3.0pt}{0.0pt}\pgfsys@lineto{0.0pt}{6.0pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}{\hbox{ \leavevmode\hbox to2.45pt{\vbox to4.45pt{\pgfpicture\makeatletter\hbox{\hskip 0.22499pt\lower-0.22499pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{} {{}{}}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.45pt}\pgfsys@invoke{ }\pgfsys@roundcap\pgfsys@invoke{ }{}\pgfsys@moveto{2.0pt}{0.0pt}\pgfsys@lineto{0.0pt}{4.0pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}{\hbox{ \leavevmode\hbox to1.9pt{\vbox to3.4pt{\pgfpicture\makeatletter\hbox{\hskip 0.2pt\lower-0.2pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{} {{}{}}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\pgfsys@roundcap\pgfsys@invoke{ }{}\pgfsys@moveto{1.5pt}{0.0pt}\pgfsys@lineto{0.0pt}{3.0pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}}j}\boldsymbol{\theta}_{\mathbin{\mathchoice{\hbox{ \leavevmode\hbox to3.6pt{\vbox to6.6pt{\pgfpicture\makeatletter\hbox{\hskip 0.3pt\lower-0.3pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{} {{}{}}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.6pt}\pgfsys@invoke{ }\pgfsys@roundcap\pgfsys@invoke{ }{}\pgfsys@moveto{3.0pt}{0.0pt}\pgfsys@lineto{0.0pt}{6.0pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}{\hbox{ \leavevmode\hbox to3.6pt{\vbox to6.6pt{\pgfpicture\makeatletter\hbox{\hskip 0.3pt\lower-0.3pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{} {{}{}}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.6pt}\pgfsys@invoke{ }\pgfsys@roundcap\pgfsys@invoke{ }{}\pgfsys@moveto{3.0pt}{0.0pt}\pgfsys@lineto{0.0pt}{6.0pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}{\hbox{ \leavevmode\hbox to2.45pt{\vbox to4.45pt{\pgfpicture\makeatletter\hbox{\hskip 0.22499pt\lower-0.22499pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{} {{}{}}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.45pt}\pgfsys@invoke{ }\pgfsys@roundcap\pgfsys@invoke{ }{}\pgfsys@moveto{2.0pt}{0.0pt}\pgfsys@lineto{0.0pt}{4.0pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}{\hbox{ \leavevmode\hbox to1.9pt{\vbox to3.4pt{\pgfpicture\makeatletter\hbox{\hskip 0.2pt\lower-0.2pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{} {{}{}}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\pgfsys@roundcap\pgfsys@invoke{ }{}\pgfsys@moveto{1.5pt}{0.0pt}\pgfsys@lineto{0.0pt}{3.0pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}}j}\right)$, for $j=1,\dots,G$.

For the BMIDAS-AGL-SS model, we place again a conjugate gamma prior on the penalty hyper-parameters. Then, for $j=1,\dots,G$, we have the following full conditional posteriors :

[TABLE]

where $\mathbf{A}_{j}:=\mathbf{Z}_{j}^{\prime}\mathbf{Z}_{j}+\tau_{j}^{-2}\mathbf{I}_{g_{j}}$, $\mathbf{C}_{j}:=\mathbf{Z}_{j}^{\prime}\left(y-\mathbf{Z}_{\mathbin{\mathchoice{\hbox{ \leavevmode\hbox to3.6pt{\vbox to6.6pt{\pgfpicture\makeatletter\hbox{\hskip 0.3pt\lower-0.3pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{} {{}{}}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.6pt}\pgfsys@invoke{ }\pgfsys@roundcap\pgfsys@invoke{ }{}\pgfsys@moveto{3.0pt}{0.0pt}\pgfsys@lineto{0.0pt}{6.0pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}{\hbox{ \leavevmode\hbox to3.6pt{\vbox to6.6pt{\pgfpicture\makeatletter\hbox{\hskip 0.3pt\lower-0.3pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{} {{}{}}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.6pt}\pgfsys@invoke{ }\pgfsys@roundcap\pgfsys@invoke{ }{}\pgfsys@moveto{3.0pt}{0.0pt}\pgfsys@lineto{0.0pt}{6.0pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}{\hbox{ \leavevmode\hbox to2.45pt{\vbox to4.45pt{\pgfpicture\makeatletter\hbox{\hskip 0.22499pt\lower-0.22499pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{} {{}{}}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.45pt}\pgfsys@invoke{ }\pgfsys@roundcap\pgfsys@invoke{ }{}\pgfsys@moveto{2.0pt}{0.0pt}\pgfsys@lineto{0.0pt}{4.0pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}{\hbox{ \leavevmode\hbox to1.9pt{\vbox to3.4pt{\pgfpicture\makeatletter\hbox{\hskip 0.2pt\lower-0.2pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{} {{}{}}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\pgfsys@roundcap\pgfsys@invoke{ }{}\pgfsys@moveto{1.5pt}{0.0pt}\pgfsys@lineto{0.0pt}{3.0pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}}j}\boldsymbol{\theta}_{\mathbin{\mathchoice{\hbox{ \leavevmode\hbox to3.6pt{\vbox to6.6pt{\pgfpicture\makeatletter\hbox{\hskip 0.3pt\lower-0.3pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{} {{}{}}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.6pt}\pgfsys@invoke{ }\pgfsys@roundcap\pgfsys@invoke{ }{}\pgfsys@moveto{3.0pt}{0.0pt}\pgfsys@lineto{0.0pt}{6.0pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}{\hbox{ \leavevmode\hbox to3.6pt{\vbox to6.6pt{\pgfpicture\makeatletter\hbox{\hskip 0.3pt\lower-0.3pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{} {{}{}}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.6pt}\pgfsys@invoke{ }\pgfsys@roundcap\pgfsys@invoke{ }{}\pgfsys@moveto{3.0pt}{0.0pt}\pgfsys@lineto{0.0pt}{6.0pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}{\hbox{ \leavevmode\hbox to2.45pt{\vbox to4.45pt{\pgfpicture\makeatletter\hbox{\hskip 0.22499pt\lower-0.22499pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{} {{}{}}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.45pt}\pgfsys@invoke{ }\pgfsys@roundcap\pgfsys@invoke{ }{}\pgfsys@moveto{2.0pt}{0.0pt}\pgfsys@lineto{0.0pt}{4.0pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}{\hbox{ \leavevmode\hbox to1.9pt{\vbox to3.4pt{\pgfpicture\makeatletter\hbox{\hskip 0.2pt\lower-0.2pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{} {{}{}}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\pgfsys@roundcap\pgfsys@invoke{ }{}\pgfsys@moveto{1.5pt}{0.0pt}\pgfsys@lineto{0.0pt}{3.0pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}}j}\right)$, $\boldsymbol{\pi}_{1}:=(\pi_{1,1},\dots,\pi_{1,G})$, $\widetilde{G}:=\sum_{j=1}^{G}g_{j}\gamma_{j}$, $d=1$ and $c=\bar{\kappa}G^{v}$ with $\bar{\kappa}=v=(1+G^{-1})$ (see Section 4.1, Assumption 4.3), and

[TABLE]

5.2 Tuning the penalty hyper-parameters

The hierarchical models presented above treat the penalty parameters as hyper-parameters, i.e. as random variables with Gamma prior distributions $\Pi(\boldsymbol{\lambda})$ and Gamma conditional posterior distributions $\Pi(\boldsymbol{\lambda}|\boldsymbol{\phi},y,\mathbf{Z})$, given a vector of parameters $\boldsymbol{\phi}$. However, the main drawback of this approach is that the resulting posterior distributions can be sensitive to the choice of the prior. Park and Casella (2008) and Kyung et al. (2010) suggest to address this issue by implementing the Monte Carlo EM algorithm (MCEM) proposed by Casella (2001), which complements the Gibbs sampler and provides marginal maximum likelihood estimates of the hyper-parameters when the marginal distribution $f(y|\boldsymbol{\lambda},\mathbf{Z})=\int f(y|\boldsymbol{\phi},\boldsymbol{\lambda},\mathbf{Z})\Pi(\boldsymbol{\phi}|\boldsymbol{\lambda})d\boldsymbol{\phi}$ is not available in closed form. This can be easily obtained by repeatedly maximizing $N$ times the expectation function

[TABLE]

and updating $\boldsymbol{\lambda}$ at each $n$-th iteration. A run of the Gibbs sampler is nevertheless required to simulate from the intractable distribution $\Pi(\boldsymbol{\phi}|\boldsymbol{\lambda}^{(n)},y,\mathbf{Z})$.

Although very attractive, the MCEM algorithm may be extremely costly from a computational point of view, as a large number of Monte Carlo iterations, each requiring a fully converged Gibbs sampling, is usually needed to maximize the marginal log-likelihood $\log f(y|\boldsymbol{\lambda},\mathbf{Z})$. Hence, a serious trade-off between computational efficiency ($N$ Monte Carlo iterations) and accuracy of the results ($S$ Gibbs iterations) may arise. In the present framework, careful attention must be paid to this feature, because the computational burden implied by the Group Lasso increases dramatically as the number of predictors increases (Yuan and Lin, 2006). To deal with this issue, we adopt an alternative Empirical Bayes approach that relies on the so-called internal adaptive MCMC algorithms (see Atchadé et al., 2011). Within this family of algorithms, the specific class of controlled MCMC resorts to stochastic approximation algorithms to solve maximization problems when the likelihood function is intractable, by mimicking standard iterative methods such as the gradient algorithm. This approach is therefore computationally efficient, because it requires only a single Monte Carlo run $(N=1)$. Following Atchadé (2011), using a stochastic approximation to solve the maximization problem and updating $\boldsymbol{\lambda}$ in each Gibbs iteration $s=1,\dots,S$, the solution to the EM algorithm takes the form:

[TABLE]

where $H(\boldsymbol{\lambda},\boldsymbol{\phi}):=\nabla_{\boldsymbol{\lambda}}\log\left[f(y|\boldsymbol{\phi},\boldsymbol{\lambda},\mathbf{Z})\Pi(\boldsymbol{\phi}|\boldsymbol{\lambda})\right]=\nabla_{\boldsymbol{\lambda}}\log\Pi(\boldsymbol{\phi}|\boldsymbol{\lambda})$, as the likelihood does not usually depend on the hyper-parameters $\boldsymbol{\lambda}$, and $a^{(s)}$ is a step-size taking a Robbins-Monro form $a^{(s)}=1/s^{q}$, with $q\in(0.5,1)$ (Lange, 1995). If the integral $\int H(\boldsymbol{\lambda}^{(s)},\boldsymbol{\phi})\Pi(\boldsymbol{\phi}|\boldsymbol{\lambda}^{(s)},y,\mathbf{Z})d\boldsymbol{\phi}$ is approximated by $H(\boldsymbol{\lambda}^{(s)},\boldsymbol{\phi}^{(s+1)})$, we get an approximate EM algorithm, where both E- and M-steps are approximately implemented. Hence, marginal maximum likelihood estimates of the hyper-parameters, $\widehat{\boldsymbol{\lambda}}$, and draws from the posterior distribution of the parameters, $\Pi(\boldsymbol{\phi}|\widehat{\boldsymbol{\lambda}},y,\mathbf{Z})$, are both obtained using a single run of the Gibbs sampler. In the present framework, making the transformation $\boldsymbol{\omega}=\frac{1}{2}\log(\boldsymbol{\lambda})$, from Section 3 the function $H(\boldsymbol{\omega},\boldsymbol{\phi})=\nabla_{\boldsymbol{\omega}}\log\Pi(\boldsymbol{\phi}|\boldsymbol{\omega})$ takes the form:

[TABLE]

where $\mathbf{g}=(g_{1},\dots,g_{G})^{\prime}$ and $\odot$ is the element-wise product. Hence, the updating rule for $\boldsymbol{\omega}$ is:

[TABLE]

from which we get $\boldsymbol{\lambda}^{(s+1)}=\exp(2\boldsymbol{\omega}^{(s+1)})$. The algorithm can be completed by allowing for a stabilization procedure (e.g. truncation on random boundaries; Andrieu et al., 2005; Atchadé, 2011) ensuring the convergence of $\boldsymbol{\lambda}$ and the posterior distribution of $\boldsymbol{\phi}$ towards $\widehat{\boldsymbol{\lambda}}$ and $\Pi(\boldsymbol{\phi}|\widehat{\boldsymbol{\lambda}},y,\mathbf{Z})$, respectively. Details on the stabilization algorithm are reported in Appendix A.3.

5.3 Numerical illustration

We illustrate the main features and the computational advantage of the proposed methodology using simulated data generated from the MIDAS model (1) with $K=4$, $m=3$, and $C=12$. We calibrate $B\left(c;\boldsymbol{\theta}\right)$ such that the MIDAS weights decay monotonically to almost zero after four periods. The regressors and the error term are i.i.d. draws from a standard normal distribution of length $T=500$. Only the second regressor $(k=2)$ enters the active set, with slope coefficient $\beta_{2}=1$, while $\beta_{1}=\beta_{3}=\beta_{4}=0$. We estimate the models presented in Section 3 using $p=3$ and we tune the penalty hyper-parameters $\boldsymbol{\lambda}$ using the stochastic approximation approach. We update $\boldsymbol{\lambda}$ in a single run of the Gibbs sampler by drawing $S=400,000$ samples. The analysis is carried out using MATLAB R2017a on a workstation with a 2.50GHz Intel Core i7-6500U CPU.

The evolution of $\boldsymbol{\lambda}$ across iterations (starting with $\lambda_{j}^{(0)}=1$) is reported in the first panel of Figure 1. Each point in the plots represents the $2000$th update of $\boldsymbol{\lambda}$ provided by the stochastic approximation approach for the BMIDAS-AGL model (solid lines) and the BMIDAS-AGL-SS model (dotted lines). For both models, the hyper-parameters converge to fairly similar values. However, while the convergence is steady and extremely fast for the active variable, the BMIDAS-AGL model displays slower convergence for the penalty terms of the inactive set compared to the BMIDAS-AGL-SS. Hence, it turns out that the spike-and-slab prior may not only improve the sparse recovery ability of the model but also enhance the convergence of the penalty hyper-parameters, reducing the variance of the posterior distribution for the parameters of the inactive set around the zero-point mass when draws are assigned (even with some low probability) to the slab part of the model. This is indeed confirmed from the inspection of the posterior densities of $\boldsymbol{\beta}$ in the second panel of Figure 1. The two models feature correct variable selection and consistent estimates of the regression coefficients, with virtually identical densities for $\beta_{2}$ and largest mass at zero for $\beta_{1}$, $\beta_{3}$, and $\beta_{4}$. However, compared to the BMIDAS-AGL, the BMIDAS-AGL-SS model displays the lowest variation around the point mass at exactly zero for the parameters of the inactive set. Finally, these outcomes are compared to those obtained by tuning the penalty hyper-parameters of the BMIDAS-AGL model using the MCEM algorithm with a fairly reasonable amount of Monte Carlo runs ($N=200$) and Gibbs draws ($S=50,000$). Looking at the posterior densities, the results for the MCEM algorithm (dashed lines) appear almost indistinguishable from those obtained using the stochastic approximation approach, with the only exception of $\beta_{4}$. However, the computational burden differs substantially across algorithms: for this simple simulation experiment and the settings described above, the analysis is performed in less than 2 minutes with stochastic approximations, against 30 minutes required by the MCEM algorithm.

6 Monte Carlo experiments

6.1 Design of the experiments

We evaluate the small-sample performance of the proposed models through Monte Carlo experiments. For this purpose, we consider the following DGP, involving $K=\{30,50\}$ predictors sampled at frequency $m=3$ and $T=200$ in-sample observations:

[TABLE]

where $\widetilde{B}\left(c;\boldsymbol{\theta}\right)$ denotes the normalized weigths (i.e. summing up to 1). Following Andreou et al. (2010), we investigate three alternative weighting schemes that correspond to fast-decaying weights (DGP 1), slow-decaying weights (DGP 2), and near-flat weights (DGP 3). These three weighting schemes are represented in Figure 2. In all simulations we set the lag length $C=24$. Note that the same weighting structure applies to all the predictors entering the DGPs. Further, for ease of analysis we assume $h=0$, i.e. a nowcasting model with high-frequency information fully matching the low frequency. In this specification, $\epsilon_{t}$ and $\boldsymbol{\varepsilon}_{t}:=(\varepsilon_{1,t},\dots,\varepsilon_{K,t})^{\prime}$ are i.i.d. with distribution:

[TABLE]

where $\boldsymbol{\Sigma}_{\varepsilon}$ has elements $\sigma_{\varepsilon}^{|k-k^{\prime}|}$, such that the diagonal elements are equal to one and the off-diagonal elements control for the correlation between $x_{k,t}^{(m)}$ and $x_{k^{\prime},t}^{(m)}$, with $k\neq k^{\prime}$. We set $\sigma_{\varepsilon}=\{0.50,0.95\}$, i.e. from moderate to extremely high correlation structure in the design matrix $\mathbf{x}_{t}^{(m)}$. As for the parameters in the DGP, we choose $\alpha=0.5$, $\mu=0.1$, $\rho=0.9$, and $\boldsymbol{\beta}=(0,0.3,0.5,0,0.3,0.5,0,0,0.8,\mathbf{0})^{\prime}$. The latter implies that only five out of $K$ predictors are relevant. Conditional on these parameters, we set $\sigma$ such that the noise-to-signal ratio of the mixed-frequency regression is 0.20.

We estimate our BMIDAS-AGL and BMIDAS-AGL-SS models with Almon lag polynomial order $p=3$ and $r=2$ endpoint restrictions (both tail and derivative; see Section 2.1). The hyper-parameters $\boldsymbol{\lambda}$ are tuned using the stochastic approximation approach described in Section 5.2, with step-size $a^{(s)}=1/s^{0.8}$ (preliminary results suggest that this sequence is sufficient to achieve convergence). We set the number of Monte Carlo replications at $R=500$. The Gibbs sampler is run for $S=300,000$ iterations, with the first $100,000$ used as a burn-in period, and every 10th draw is saved.

6.2 Estimation and selection performance

We evaluate the estimation performance of our models by computing the average mean squared error (MSE), the average variance (VAR), and the average squared bias (BIAS2) over $R$ Monte Carlo replications and the full set of $K$ estimated overall slope parameters $\widehat{\boldsymbol{\beta}}$.555For $R$ Monte Carlo replications, $K$ variables, and $S$ Gibbs draws, we have that:

\text{MSE}=\text{VAR}+\text{BIAS{}^{2}}=\frac{1}{RKS}\sum_{i=1}^{R}\sum_{k=1}^{K}\sum_{s=1}^{S}\left[\widehat{\beta}_{k,i}^{(s)}-\mathbb{E}\left(\widehat{\beta}_{k,i}\right)\right]^{2}+\frac{1}{RK}\sum_{i=1}^{R}\sum_{k=1}^{K}\left[\mathbb{E}\left(\widehat{\beta}_{k,i}\right)-\beta_{k}\right]^{2},

where $\mathbb{E}\left(\widehat{\beta}_{k}\right)=\frac{1}{S}\sum_{s=1}^{S}\widehat{\beta}_{k}^{(s)}$ and $\widehat{\beta}_{k}=\widehat{\boldsymbol{\theta}}_{k}^{\prime}\mathbf{Q}\boldsymbol{\iota}_{C}$, for $k=1,\dots,K$. Note that for the BMIDAS-AGL-SS model we use the median estimator.

To evaluate the selection ability of the models, we consider the True Positive Rate (TPR), the False Positive Rate (FPR), and the Matthews correlation coefficient (MCC), the latter measuring the overall quality of the classification. We follow different approaches to implement variable selection. For the BMIDAS-AGL model, we rely on the credible interval criterion suggested by Kyung et al. (2010), which excludes a predictor $k$, for $k=1,\dots,K$, from the estimated active set if the credible interval, at say 95% level, of the posterior distribution of its slope coefficient $\widehat{\beta}_{k}$ includes zero. For the BMIDAS-AGL-SS model, we resort to the posterior median estimator (Barbieri and Berger, 2004), that is, under some conditions, a soft thresholding estimator presenting model selection consistency and optimal asymptotic estimation rate (Xu and Ghosh, 2015).

6.3 Forecast evaluation

Forecasts are obtained from the following posterior predictive density for $y_{T|T-h}$:

[TABLE]

where $\Pi(\boldsymbol{\phi},\boldsymbol{\lambda}|\mathcal{D})$ denotes the joint posterior distribution of the BMIDAS parameters conditional on past available information, $\mathcal{D}$. According to the framework described in Sections 3 and 5.2, draws $y_{T|T-h}^{(s)}$, $s=1,\dots,S$, from the predictive distribution (14) can be obtained directly from the Gibbs sampler.666It is worth noting that we do not condition on a fixed value $\widehat{\boldsymbol{\lambda}}$, such as the maximum likelihood estimate that can be obtained, for instance, by averaging over the Gibbs samples of $\boldsymbol{\lambda}$, because this would ignore the uncertainty around the estimate of the penalty parameters. This leads to a distribution of predictions that can be used for out-of-sample evaluation of the model. Point forecasts are computed by averaging over these draws, i.e. $\widehat{y}_{T|T-h}=(S-\bar{s}+1)^{-1}\sum_{s=\bar{s}+1}^{S}y_{T|T-h}^{(s)}$, where $\bar{s}$ is the last burn-in iteration, and evaluated through the average root mean squared forecast error (RMSFE) over the $R$ Monte Carlo replications. Further, since draws from the predictive density are available, an evaluation of the entire predictive distribution is performed through the average log-score (LS), i.e. the average of the log of the predictive likelihood evaluated at the out-turn of the forecast (Mitchell and Wallis, 2011), and the average continuously ranked probability score (CRPS), which measures the average distance between the empirical CDF of the out-of-sample observations and the empirical CDF associated with the predictive density of each model (Gneiting and Raftery, 2007).

6.4 Monte Carlo results

Simulation results are reported in Table 1 and point to a number of interesting features. First, the models perform overall quite similarly in terms of MSE, although the BMIDAS-AGL-SS seems to perform somewhat better across DGPs by mainly providing the smallest bias. This leads to highest TPR and lowest FPR for this model, entailing better classification of active and inactive sets across simulations. Second, the MSE increases substantially with the degree of correlation in the design matrix (governed by the value of $\sigma_{\varepsilon}$). However, the MSE tends to decrease with more irrelevant predictors in the DGP. To understand the latter result, it is useful to look at the bias/variance breakdown of the MSE stemming from both the active $(\mathcal{A})$ and inactive set $(\mathcal{A}^{c})$ reported in Figure 3, where we refer to our models as AGL and AGL-SS. Results suggest that while the contribution of the inactive set to the total MSE is broadly stable when $K$ increases, the contribution of the active set decreases for both bias and variance. These findings can nevertheless be attributed, at least in part, to the lower relative weight of the active set $(w_{\mathcal{A}})$ in the total bias and variance, as the number of relevant predictors $(K_{\mathcal{A}})$ is fixed while the number of irrelevant predictors $(K_{\mathcal{A}^{c}})$ is allowed to increase.777For $w_{\mathcal{A}}=K_{\mathcal{A}}/K$ and $w_{\mathcal{A}^{c}}=1-w_{\mathcal{A}}$, we have that:

\text{MSE}=w_{\mathcal{A}}\text{MSE}(\mathcal{A})+w_{\mathcal{A}^{c}}\text{MSE}(\mathcal{A}^{c})=w_{\mathcal{A}}\left[\text{VAR}(\mathcal{A})+\text{BIAS{}^{2}}(\mathcal{A})\right]+w_{\mathcal{A}^{c}}\left[\text{VAR}(\mathcal{A}^{c})+\text{BIAS{}^{2}}(\mathcal{A}^{c})\right].

Looking at the MSE by active/inactive set reported in Table 1 (columns 5 and 6), we note that the $\text{MSE}(\mathcal{A})$ appears broadly stable or it increases only moderately with more irrelevant variables in DGPs 1 and 2.888For DGP 3, the $\text{MSE}(\mathcal{A})$ tends to increase more substantially with more irrelevant predictors. One explanation is that the linear restrictions imposed to the lag polynomials are incorrect under this DGP. Hence, our interpretation is that the performance of our BMIDAS models in selecting and estimating the coefficients in the active set seems only marginally affected by the increase in the degree of sparsity.

This result is confirmed by the TPR, which is relatively high and hovers around 90% for moderate correlation, and it’s overall stable across different values of $K$, suggesting that the models can select the correct sparsity pattern with a high probability even in finite samples. It is worth noting that the TPR drops to 30-50% with very high correlation in the design matrix, while the FPR remains overall very low. Note that this outcome is nevertheless not unexpected, as the Group Lasso can address the issue of strong collinearity within the lag polynomials but is not designed to handle strong collinearity between the high-frequency regressors. For comparison purposes, we also report in Figure 3 the breakdown of the MSE for the Bayesian MIDAS Adaptive Lasso (BMIDAS-AL) and the Oracle BMIDAS, the former estimated following the same approach as in Sections 3 and 5 and the latter using the algorithm described in Pettenuzzo et al. (2016) on the set of relevant variables only.999We consider the same restricted Almon lag polynomial as for our models. Further, for the Oracle BMIDAS, we follow Pettenuzzo et al. (2016) and we use relatively diffuse priors on both the coefficient covariance matrix and the regression variance. As for the prior mean coefficients, we set all the coefficients but the intercept to zero. A visual inspection of Figure 3 seems to confirm the intuition discussed in Section 2.2, that the AL may not be suited in the present framework. Results suggest that higher MSE provided by the AL can be mainly attributed to higher bias and variance in the inactive set. It is worth noting that while these findings appear quite conclusive for DGPs 1 and 3, evidence is less clear-cut for DGP 2, especially when high correlation in the design matrix is considered. When compared to the Oracle, our models perform fairly well overall. Not surprisingly, in most cases the main difference lies in the bias of the active set, as the Group Lasso would typically trade off more bias for less variance.

Third, the in-sample results shown in Table 1 deteriorate when the DGP with near-flat weights is considered, and mostly when $\sigma_{\varepsilon}=0.95$. This happens because the linear restrictions imposed on the lag polynomials force the weighting structure to tail off to zero, while DGP 3 assumes that the weighting scheme under the null is almost uniform over the lag window. It follows that relaxing the restrictions on the Almon lag polynomial should lead to an improvement of the results under this DGP. However, it is not clear how much of the selection results under DGPs 1 and 2 can be attributed to the imposed linear restrictions. Table 2 provides an answer to these questions by reporting the difference in TPR, FPR, and MCC obtained with restricted $(r=2)$ and unrestricted $(r=0)$ Almon lag polynomials. For DGPs with fast- or slow-decaying weights, the results suggest that imposing correct linear restrictions that are valid under the null seems to improve the selection ability of the models. The gain in terms of TPR ranges 10-25 percentage points for moderate correlation and 5-15 percentage points for very high correlation in the design matrix, while the gain in terms of FPR ranges 1-4 percentage points. Interestingly enough, the results show that the BMIDAS-AGL-SS model is relatively less affected than the BMIDAS-AGL model when correct linear restrictions are relaxed. For the DGP with near-flat weights we observe an opposite outcome, as expected. However, the magnitude of these results must be considered with care, as the number of relevant and irrelevant predictors in the simulated true model is here strongly asymmetric.

Finally, we report in Table 3 the forecasting performance of our models, in relative terms with respect to the Oracle for ease of comparison.101010We compute RMSFE and CRPS ratios, such that values greater than one indicate that the Oracle performs best. For the LS, we compute the log-score differentials, such that negative values indicate that the Oracle performs best. The results are broadly in line with those stemming from the in-sample analysis, and suggest that our BMIDAS models perform quite similarly in terms of point and density forecasts, although the BMIDAS-AGL-SS model seems to perform best overall. The performance of our models deteriorates substantially with higher correlation in the design matrix, although a higher average variance of the error process (see $\overline{\sigma}$ in Table 1) might explain, at least in part, the large differences spotted across DGPs. Nevertheless, the performance appears only mildly sensitive to the number of irrelevant predictors ($K_{\mathcal{A}^{c}}$ increasing). We compare these results to those obtained from a set of alternative penalized regressions. In addition to the BMIDAS-AL, we consider the following penalized MIDAS models, in the spirit of Uematsu and Tanaka (2019): the Lasso (L), the Elastic-Net (EN; Zou and Hastie, 2005), and two folded-concave penalizations, such as the smoothly clipped absolute deviation (SCAD; Fan and Li, 2001) and the minimax concave penalty (MC$+$; Zhang, 2010). In particular, the SCAD and MC$+$ penalties do not require the irrepresentable condition to achieve selection consistency and can attenuate the estimation bias problem of convex penalty functions (Fan and Lv, 2011). Compared to these alternative penalized regressions, our models seem to provide a substantially higher predictive performance. Results look quite significant for DGPs 1 and 3, irrespective of the degree of sparsity and the correlation in the design matrix, but appear less clear-cut for DGP 2. Similarly to Uematsu and Tanaka (2019), we point out that the Lasso often performs best among the set of competing penalties. Taken jointly, these findings are supportive of the proposed Adaptive Group Lasso prior in the present penalized mixed-frequency framework, and suggest that the chosen modeling approach, although computationally demanding, would likely pay off in terms of predictive accuracy.

7 Empirical application: forecasting U.S. GDP

We apply the proposed Bayesian MIDAS penalized regression to U.S. GDP data. Following the literature, we consider the annualized quarterly growth rate of GDP, $y_{t}=4\log(Y_{t}/Y_{t-1})$. As for the predictors, we consider 125 macroeconomic series sampled at monthly frequency and extracted from the FRED-MD database (McCracken and Ng, 2016).111111We consider the entire FRED-MD database (at the time of writing), except for the series of non-borrowed reserves of depository institutions, because of the extreme changes observed since 2008 (see Uematsu and Tanaka, 2019), as well as the effective Federal Funds rate and the spread with the 10-year government bond rate, as we already use these series at a daily frequency. Further, we also consider a set of daily and weekly financial data, which have proven to improve short- to medium-term macro forecasts (Andreou et al., 2013; Pettenuzzo et al., 2016; Adrian et al., 2019): the effective Federal Funds rate; the interest rate spread between the 10-year government bond rate and the Federal Funds rate; returns on the portfolio of small minus big stocks considered by Fama and French (1993); returns on the portfolio of high minus low book-to-market ratio stocks studied by Fama and French (1993); returns on a winner minus loser momentum spread portfolio; the Chicago Fed National Financial Conditions Index (NFCI), and in particular its three sub-indexes (risk, credit, and leverage). Finally, we consider the Aruoba-Diebold-Scotti (ADS) daily business conditions index (Aruoba et al., 2009) to track the real business cycle at high frequency. Overall, the total number of predictors entering the models is $K=134$. The data sample starts in 1980Q1, and we set $\underline{T}=\text{2000Q1}$ and $\overline{T}=\text{2017Q4}$ the first and last out-of-sample observations, respectively. Estimates are carried-out recursively using an expanding window, and $h$-step-ahead posterior predictive densities are generated from (14) through a direct forecast approach. We hence dispose of $(\overline{T}-\underline{T}+1)=72$ out-of-sample observations. We consider several forecast horizons, leading to a sequence of 3 nowcasting ($h=0,1/3,2/3$) and 5 short- and medium-term forecasting ($h=1,4/3,5/3,2,4$) exercises. We keep the empirical application as much realistic as possible, but for ease of analysis we do not take into account real-time issues (ragged/jagged-edge data and revisions). The dataset is hence compiled using the latest vintages available at the time of writing.

Forecasts are compared to those from a benchmark model represented by a simple random-walk (RW). Point and density forecasts are evaluated again by the means of relative RMSFE ratios, average LS differentials, and average CRPS ratios. For both the RMSFE and the CRPS, values less than one shall indicate that our penalized mixed-frequency models outperform (in either point or density forecast sense) the RW. For the LS, positive values shall indicate that our models produce more accurate density forecasts than the RW. To account for sample uncertainty underlying the observed forecast differences, we report results for the Diebold and Mariano (1995) and West (1996) test (DMW hereafter), which posits the null hypothesis of unconditional equal predictive accuracy. The resulting test statistic is computed using HAC standard errors (for $h>1$) and a small-sample adjustment to the consistent estimate of the variance, and compared with critical values from the Student’s $t$ distribution with $(\overline{T}-\underline{T})$ degrees of freedom (Harvey et al., 1997). As a robustness check, we further consider forecasts from the following competing models:

AR(1) model.

• 2.

BMIDAS-AL model.

• 3.

Penalized MIDAS regressions with Lasso, Elastic-Net, SCAD, and MC$+$ penalties. Regularization parameters are tuned by optimizing the BIC over a fine grid of values.

• 4.

Combination of $K$ single-indicator Bayesian MIDAS models (BMIDAS-comb) as in Pettenuzzo et al. (2016). The combination weights are computed using a discounted version of the optimal prediction pool proposed by Geweke and Amisano (2011), with discount factor $\delta=0.9$ (Stock and Watson, 2004; Andreou et al., 2013).

• 5.

Bayesian Model Averaging and Model Selection with MIDAS (BMA-MIDAS and BMS-MIDAS, respectively). We consider a standard $g$-BRIC prior and a reversible-jump MC3 algorithm, the latter modified to account for groups of lag polynomials in the addition/deletion/swaption moves. This ensures that model proposals are based on selection of individual predictors rather than isolated terms of lag polynomials.

• 6.

Dynamic Factor and targeted-Factor Bayesian MIDAS regressions (Factor-BMIDAS and $t$-Factor-BMIDAS, respectively). We set the number of dynamic factors as in Bai and Ng (2002, 2007). Targeted predictors are selected using the Elastic-Net soft-thresholding rule (Bai and Ng, 2008). Factors are estimated using either principal components (PCA) or the quasi-maximum likelihood approach (QML) of Doz et al. (2012).

[FIGURE:]

[FIGURE:]

All the models considered in the application include one lag of the growth rate of GDP, which is hence excluded from the selection procedures. To match the sample frequencies, we consider again a restricted Almon lag polynomial, with $p=3$ and $r=2$ endpoint restrictions, and twelve months of past high-frequency observations. As for the MCMC, the Gibbs sampler is run for $S=600,000$ iterations, with the first $200,000$ used as a burn-in period, and every 10th draw is saved. For the BMA/BMS model, we increase the number of iterations to $4,000,000$, in order to let the algorithm sufficiently explore the model space, which is fairly vast in the current application.

Results are reported in Tables 4 and 5. In the first row of Table 4, we report RMSFE, LS and CRPS for the benchmark RW, while in the other rows we report the relative scores for all models compared to the RW. For the short horizons between $h=0$ and $h=1$ (i.e. nowcast and 1-step-ahead forecast), the results suggest that the BMIDAS-AGL-SS model systematically outperforms all the competitors, with point and density predictive gains often statistically significant (at 10% level; see Table 5). The best competing results are given by the BMA/BMS-MIDAS models for $h=0,1/3,2/3$ and the $t$-Factor-BMIDAS models for $h=1$. The predictive gains provided by the BMIDAS-AGL-SS are admittedly small compared to these best competing models, but often fairly large compared to the other models considered. For $h>1$, the BMIDAS-AGL-SS model cannot provide the best predictive outcomes, but it is often ranked among the best models and can often outperform a large number of competitors, although the number of rejections in the DMW test edges down to only a few. Conversely, results for the BMIDAS-AGL model are overall more disappointing, as the model seems to perform reasonably well only up to $h=1$.

A number of additional interesting features arise from the out-of-sample analysis. First, results for the competing penalized MIDAS regressions (Lasso, EN, SCAD, and MC$+$) are very similar, and sometimes virtually identical, to each other but also fairly volatile, as their predictive performance can deteriorate substantially from one forecast horizon to another. Second, for $h>1$ best outcomes are provided by targeted-Factor models, except for $h=4$ where the MIDAS-L regression performs slightly better. The AR is never ranked first but, for relatively long horizons, this model becomes hard to beat, which is broadly in line with the empirical findings reported in previous studies. Finally, Factor MIDAS regressions and MIDAS penalized regressions seem to perform quite similarly, overall. These findings are broadly in line with those reported by Uematsu and Tanaka (2019), but they additionally reveal that some prior sparse selection (i.e. targeting) seems necessary to let Factor models systematically outperform penalized regressions. Nevertheless, to our opinion these findings lack of generality on this point, leaving the debate over sparse $vs$ dense modelling in presence of possibly highly correlated predictors still open (Giannone et al., 2017).

Figure 4 reports the variables inclusion probabilities obtained from the BMIDAS-AGL-SS model. Given the large number of variables considered in the application, for ease of exposition we aggregate these probabilities according to the nature of the regressors and/or their frequency, as well as to the classification used by McCracken and Ng (2016). The patterns reported in the figure show a systematic inclusion with very high probability of the ADS index for all forecast horizons between $h=0$ and $h=1$. The model tends to select also a bunch of high-frequency predictors related to consumption, output, and inventories, but with much lower probability. This is not unsurprising, as the ADS index already contains signals stemming from the real economy (initial jobless claims, real GDP, payroll employment, industrial production, real personal income less transfers, and real manufacturing and trade sales). The housing market seems to play a very limited role in the model ($h=1$), while virtually no financial indicators are selected for nowcasting purposes. However, this feature tends to progressively attenuate for $h>1$, and financial variables (financial condition indexes and stock market) are selected with somewhat higher probability up to $h=2$. This result seems broadly in line with recent literature (Andreou et al., 2013) and suggests that financial variables may convey some, although here very limited, short-term leading information which goes beyond the predictive content of real indicators.

8 Concluding remarks

We propose a new approach to modeling and forecasting mixed-frequency regressions (MIDAS) that addresses the issue of simultaneously estimating and selecting relevant high-frequency predictors in a high-dimensional environment. Our approach is based on MIDAS regressions resorting to Almon lag polynomials and an adaptive penalized regression approach, namely the Group Lasso objective function. The proposed models rely on Bayesian techniques for estimation and inference. In particular, the penalty hyper-parameters driving the model shrinkage are automatically tuned via an Empirical Bayes algorithm based on stochastic approximations. We establish the posterior contraction rate for the in-sample and the out-of-sample prediction error, as well as the consistency of the marginal posterior of parameters. Simulations show that the proposed models present very good in-sample and out-of-sample performance. When applied to a forecasting model of the U.S. GDP growth with high-frequency real and financial predictors, the results suggest that our models produce significant short-term predictive gains compared to several alternative models. Our findings point to a very limited short-term predictive content of high-frequency financial variables, which is broadly in line with the existing literature.

The models presented in the present paper could be extended in several ways, such as time-varying lag polynomials with stochastic volatility error dynamics (Carriero et al., 2015; Pettenuzzo et al., 2016), as well as quantile mixed-frequency regressions (Lima et al., 2020). Further, recent research has been focusing on factor-adjusted sparse regressions to deal with high-dimensional regressions and highly correlated predictors (Kneip and Sarda, 2011; Fan et al., 2020). We believe that these extensions to our models represent interesting paths for future research.

Appendix

A.1 Marginal prior for the Bayesian Adaptive Group Lasso

In this section we detail the computations to get the marginal prior for $\boldsymbol{\theta}$ conditional on $\sigma^{2}$ for the prior given in Section 3.1 or conditional on $(\sigma^{2},\pi_{0})$ for the prior given in Section 3.2. Given the specification of these priors, we have to marginalise the $\mathcal{N}(\mathbf{0},\sigma^{2}\tau_{j}^{2}\mathbf{I}_{g_{j}})$ distribution with respect to the prior of $\tau_{j}^{2}$:

[TABLE]

where $GIG(\tau_{j}^{2};a,b,p)$ denotes the pdf of a Generalized Inverse Gaussian distribution with parameters $a$, $b$ and $p$, and $K_{1/2}(\lambda_{j}\|\boldsymbol{\theta}_{j}\|_{2}/\sigma)=\sqrt{\frac{\pi\sigma}{2\lambda_{j}\|\boldsymbol{\theta}_{j}\|_{2}}}\exp\{-\lambda_{j}\|\boldsymbol{\theta}_{j}\|_{2}/\sigma\}$ is a modified Bessel function of the second kind. Remark that the constant $C_{j}$ can also be written as $C_{j}=\pi^{-g_{j}/2}\Gamma(g_{j}/2)/(2\Gamma(g_{j}))$ by using the property of the Gamma function $\Gamma(g_{j}/2+1/2)=2^{1-g_{j}}\sqrt{\pi}\Gamma(g_{j})/\Gamma(g_{j}/2)$. Hence, by integrating out $\tau_{j}^{2}$ we get:

[TABLE]

for the prior of the BMIDAS-AGL model in Section 3.1, or

[TABLE]

for the prior of the BMIDAS-AGL-SS model in Section 3.2. In both cases, $\textrm{M-Laplace}(\boldsymbol{\theta}_{j};\mathbf{0},\sigma/\lambda_{j})$ denotes a $g_{j}$-dimensional Multi-Laplace distribution.

A.2 Proofs

In this section we prove the theorems in Section 4. Technical lemmas and further results on dimension recovery and distributional approximation of the posterior of $\boldsymbol{\theta}$ are reported in the Supplementary Appendix. Recall the notation $y:=(y_{1},\ldots,y_{T})^{\prime}$ for the $T$-vector of observations on $y_{t}$ and $\mathbf{Z}$ for the $T\times\widetilde{g}$ matrix of observations on $\mathbf{Z}_{t}^{(m)}$. Moreover, hereafter we use the following notation. Let $f$ (resp. $f_{0}$) denote the joint conditional probability density with respect to the Lebesgue measure of $y$ given $(\mathbf{Z},\boldsymbol{\theta},\sigma^{2})$ for a generic value of the parameters $(\boldsymbol{\theta}^{\prime},\sigma^{2})^{\prime}$ (resp. of $y$ given $(\mathbf{Z},\boldsymbol{\theta}_{0},\sigma_{0}^{2})$ for the true value of the parameters). We denote by $\mathbf{E}_{f}[\cdot]$ the expectation taken with respect to the probability distribution with Lebesgue density $f$ and by $\mathbf{E}_{0}[\cdot]$ (resp. $\mathbf{P}_{0}$) the expectation (resp. the probability) taken with respect to the true conditional distribution of $y$ given $(\mathbf{Z},\boldsymbol{\theta}_{0},\sigma_{0}^{2})$. For two probability densities $f_{1}$ and $f_{2}$, the Kullback-Leibler divergence is denoted by $K(f_{1},f_{2}):=\int f_{1}\log(f_{1}/f_{2})$ while $V(f_{1},f_{2}):=\int f_{1}(\log(f_{1}/f_{2})-K(f_{1},f_{2}))^{2}$ denotes the Kullback-Leibler variation. The average Rényi divergence of order $1/2$ between two $T$-dimensional densities $f_{1}$ and $f_{2}$ is denoted by $d(f_{1},f_{2}):=-\frac{1}{T}\log\int\sqrt{f_{1}f_{2}}$.

Let $\widetilde{s}_{0}:=\max\{s_{0},\log(T)/\log(G)\}$ and, for a vector $\mathbf{v}$, we denote by $s_{\mathbf{v}}$ the number of groups in $\mathbf{v}$ with nonzero components and by $S_{\mathbf{v}}$ the indices of these groups. For two constants $M_{2}>0$ and $C_{2}>0$ define the set $\Theta_{\widetilde{s}_{0}}:=\{\boldsymbol{\theta}\in\Theta;\,s_{\boldsymbol{\theta}}<M_{2}\widetilde{s}_{0}\}$ and the sieve

[TABLE]

where $\xi:=T\widetilde{s}_{0}\log(G)/\|\mathbf{Z}\|_{o}$. Let $f=\mathcal{N}_{T}(\mathbf{Z}\boldsymbol{\theta},\sigma^{2}\mathbf{I}_{T})$ and $f_{0}=\mathcal{N}_{T}(\mathbf{Z}\boldsymbol{\theta}_{0},\sigma_{0}^{2}\mathbf{I}_{T})$. The notation $f\in\mathcal{F}_{T}(C_{2})$ must be understood as: $f=\mathcal{N}_{T}(\mathbf{Z}\boldsymbol{\theta},\sigma^{2}\mathbf{I}_{T})$ with $(\boldsymbol{\theta},\sigma^{2})\in\mathcal{F}_{T}(C_{2})$. We denote by $N(\eta,\mathcal{F},\rho)$ the $\eta$-covering number of a set $\mathcal{F}$ with respect to a metric $\rho$, which is the minimal number of $\eta$-balls in the $\rho$-metric needed to cover the set $\mathcal{F}$.

A.2.1 Proof of Theorem 4.1

The proof of this theorem follows Ning et al. (2020). As in the latter, it is made of two parts. In the first part we obtain the posterior contraction rate with respect to the divergence $d(f_{0},f)$. In the second part we use this result to show the result of the theorem. However, while the general idea of the proof is the same as in Ning et al. (2020), many technical intermediate steps are different due to the fact that we have a different prior distribution and different assumptions. We provide all the details for completeness.

Let us consider the set $\Theta_{\widetilde{s}_{0}}:=\{\boldsymbol{\theta}\in\Theta;s_{\boldsymbol{\theta}}<M_{2}\widetilde{s}_{0}\}$ with $M_{2}>\frac{2+C_{1}}{u}+1$, where $C_{1}>0$ is a constant that depends only on $a_{1},d,c_{\pi_{0}},\overline{C},\underline{\sigma}^{2},c_{\boldsymbol{\theta}},c_{g}$ and that satisfies (B3) in Lemma B.1.2. For every $\varepsilon>0$ let $\mathcal{A}_{\varepsilon}:=\{(\boldsymbol{\theta},\sigma^{2})\in\Theta_{\widetilde{s}_{0}}\times\mathbb{R}_{+};\,d(f_{0},f)>\varepsilon\}$. Then, for every $\varepsilon>0$:

[TABLE]

The second term on the right hand side goes to zero by Lemma B.1.4 in the Supplementary Appendix. We then focus on the first term. By Lemma B.1.5 in the Supplementary Appendix there exists a test $\phi_{T}$ such that for some constant $M_{1}>2M_{2}(1+8c_{g})+4+2C_{2}$ that does not depend on $(\boldsymbol{\theta}_{0},\sigma_{0}^{2})$,

[TABLE]

Hence, by using the first result in (A5) we get the upper bound:

[TABLE]

and we just need to upper bound the second term in (A6). To do this and for the constant $C_{1}$ introduced above, define the event $\mathcal{A}_{1}:=\{\int(f/f_{0})\overline{\Pi}(\boldsymbol{\theta},\sigma^{2})d\boldsymbol{\theta}d\sigma^{2}\geq\exp\{-(1+C_{1})T\epsilon^{2}\}\}$, where $\overline{\Pi}(\boldsymbol{\theta},\sigma^{2})$ denotes the prior in (A2) marginalized with respect to the prior of $\pi_{0}$ and restricted to have support on $\mathcal{B}_{0}(\epsilon):=\{(\boldsymbol{\theta},\sigma^{2});K(f_{0},f)\leq T\epsilon^{2},V(f_{0},f)\leq T\epsilon^{2}\}$. By using the results of Lemmas B.1.1 and B.1.2 in the Supplementary Appendix we obtain that for large $T$:

[TABLE]

Moreover, let $\mathcal{A}_{\varepsilon,T}:=\{(\boldsymbol{\theta},\sigma^{2})\in\mathcal{F}_{T}(C_{2});\,d(f_{0},f)>\varepsilon\}$ where $\mathcal{F}_{T}(C_{2})$ is as defined in (A3) with $C_{2}>2C_{1}+1$ and remark that $\mathcal{A}_{\varepsilon}\subseteq\mathcal{A}_{\varepsilon,T}\cup\{(\boldsymbol{\theta},\sigma^{2})\in(\Theta_{\widetilde{s}_{0}}\times\mathbb{R}_{+})\setminus\mathcal{F}_{T}(C_{2})\}$. Then, by letting $\varepsilon=M_{1}\epsilon^{2}$, we get:

[TABLE]

by using (A5) and result (B17) of Lemma B.1.5. By putting together (A4) and (A6)-(A8) and by Lemma B.1.4 we obtain:

[TABLE]

which converges to zero for the values of $M_{1},M_{2}$ and $C_{2}$ specified above. Finally, since all the constants in the right hand side of (A10) do not depend on $(\boldsymbol{\theta}_{0},\sigma_{0}^{2})$ we conclude that

[TABLE]

We now develop the second part of the proof. Denote $\mathcal{A}_{2}:=\{(\boldsymbol{\theta},\sigma^{2})\in\mathbb{R}^{\widetilde{g}}\times\mathbb{R}_{+};\frac{1}{4}\log\left(\frac{\sigma^{2}+\sigma_{0}^{2}}{2\sigma_{0}^{2}}\right)\leq M_{1}\epsilon^{2}\}$. Remark that under Assumption 4.1 (i) this event is feasible. First, since

[TABLE]

where $f_{t}=\mathcal{N}(Z_{t}\boldsymbol{\theta},\sigma^{2})$ and $f_{0,t}=\mathcal{N}(Z_{t}^{\prime}\boldsymbol{\theta}_{0},\sigma_{0}^{2})$, then (A10) implies that uniformly on $(\boldsymbol{\theta}_{0},\sigma_{0}^{2})\in\overline{\Theta}_{0}\times[\underline{\sigma}^{2},\overline{\sigma}^{2}]$:

[TABLE]

Remark that on $\mathcal{A}_{2}$: $\frac{\sigma^{2}+\sigma_{0}^{2}}{2\sigma_{0}^{2}}\leq e^{4M_{1}\epsilon^{2}}$. Then on $\mathcal{A}_{2}$:

[TABLE]

It follows from this and the law of total probability that

[TABLE]

for every $M\geq 2M_{1}4\overline{\sigma}^{2}e^{4M_{1}\epsilon^{2}}$. Finally, (A14) converges to one by (A12) and (A13) uniformly over $(\boldsymbol{\theta}_{0},\sigma_{0}^{2})\in\overline{\Theta}_{0}\times[\underline{\sigma}^{2},\overline{\sigma}^{2}]$.

A.2.2 Proof of Theorem 4.2

By definition of $\Theta_{\widetilde{s}_{0}}$, for every $\boldsymbol{\theta}\in\Theta_{\widetilde{s}_{0}}$ it holds that $s_{(\boldsymbol{\theta}-\boldsymbol{\theta}_{0})}\leq s_{\boldsymbol{\theta}}+s_{\boldsymbol{\theta}_{0}}\leq M_{2}\widetilde{s}_{0}+s_{0}$. Therefore, $\forall\boldsymbol{\theta}\in\Theta_{\widetilde{s}_{0}}$

[TABLE]

By using this inequality to get the inequality in the fourth line below, we obtain:

[TABLE]

The two terms on the right hand side converge to zero uniformly over $(\boldsymbol{\theta}_{0},\sigma_{0}^{2})\in\overline{\Theta}_{0}\times[\underline{\sigma}^{2},\overline{\sigma}^{2}]$ due to Theorem 4.1 and Lemma B.1.4 in the Supplementary Appendix.

A.2.3 Proof of Theorem 4.3

First, remark that $\|\mathbf{Z}\|_{0}^{2}\|\boldsymbol{\theta}-\boldsymbol{\theta}\|_{2}^{2}=\max_{1\leq j\leq G}\|\mathbf{Z}_{j}\|_{op}^{2}\sum_{j=1}^{G}\|(\boldsymbol{\theta}-\boldsymbol{\theta}_{0})_{j}\|_{2}^{2}\geq\sum_{j=1}^{G}\|\mathbf{Z}_{j}\|_{op}^{2}\|(\boldsymbol{\theta}-\boldsymbol{\theta}_{0})_{j}\|_{2}^{2}\geq\sum_{j=1}^{G}\|\mathbf{Z}_{j}^{\prime}(\boldsymbol{\theta}-\boldsymbol{\theta}_{0})_{j}\|_{2}^{2}$. On the set $\Theta_{\widetilde{s}_{0}}$ the following upper bound holds:

[TABLE]

since on $\Theta_{\widetilde{s}_{0}}$, $s_{(\boldsymbol{\theta}-\boldsymbol{\theta}_{0})}\leq s_{\boldsymbol{\theta}}+s_{\boldsymbol{\theta}_{0}}\leq M_{2}\widetilde{s}_{0}+s_{0}$. We can then use this result to get the second inequality below:

[TABLE]

The last two terms converge to zero uniformly over $(\boldsymbol{\theta}_{0},\sigma_{0}^{2})\in\overline{\Theta}_{0}\times[\underline{\sigma}^{2},\overline{\sigma}^{2}]$ due to Theorem 4.1 and Lemma B.1.4 in the Supplementary Appendix.

A.2.4 Proof of Theorem 4.4

By using the definition of $f_{\mathbf{Z}_{\tau-h},y,\mathbf{Z}}(y_{\tau})$, $f(\cdot|\mathbf{Z}_{\tau-h},y,\mathbf{Z})$, $f_{\mathbf{Z}_{\tau-h},\boldsymbol{\theta}_{0},\sigma_{0}^{2}}(\cdot)$ and the Total Variation distance, we have:

[TABLE]

where we have used the Fubini’s Theorem to get the equality in the last line. Let $\mathcal{A}$ be a subset of the parameter space that will be defined below and recall the notation $\Theta_{\widetilde{s}_{0}}:=\{\boldsymbol{\theta}\in\Theta;\,s_{\boldsymbol{\theta}}<M_{2}\widetilde{s}_{0}\}$, where $\widetilde{s}_{0}:=\max\{s_{0},\log(T)/\log(G)\}$. Hence, by using the upper bound of the Total variation distance: $\|f_{\mathbf{Z}_{\tau-h},\boldsymbol{\theta},\sigma^{2}}-f_{\mathbf{Z}_{\tau-h},\boldsymbol{\theta}_{0},\sigma_{0}^{2}}\|_{TV}=2-2\int\min\{f_{\mathbf{Z}_{\tau-h},y,\mathbf{Z}}(y_{\tau}),f_{\mathbf{Z}_{\tau-h},\boldsymbol{\theta}_{0},\sigma_{0}^{2}}(y_{\tau})\}dy_{\tau}\leq 2$ we obtain:

[TABLE]

By using Lemma B.1(i) and Lemma B.5 (i) and (vi) in Ghosal and van der Vaart (2017) we get the inequalities

[TABLE]

where $d_{H}^{2}$ denotes the squared Hellinger distance and $d$ is the Renyi divergence of order $1/2$. Moreover, on $\Theta_{\widetilde{s}_{0}}$ and by assuming that $\mathbf{E}_{\mathbf{Z}_{\tau-h}}[\mathbf{Z}_{\tau-h,j}\mathbf{Z}_{\tau-h,k}^{\prime}]=0$ for every $j\neq k$, the following holds

[TABLE]

since $s_{\boldsymbol{\theta}-\boldsymbol{\theta}_{0}}\leq s_{\boldsymbol{\theta}}+s_{\boldsymbol{\theta}_{0}}\leq M_{2}\widetilde{s}_{0}+s_{0}$, $\forall\boldsymbol{\theta}\in\Theta_{\widetilde{s}_{0}}$ and where we have used the notation $\eta:=\frac{\left\|\mathbf{E}_{\mathbf{Z}_{\tau-h}}[\mathbf{Z}_{\tau-h}\mathbf{Z}_{\tau-h}^{\prime}]\right\|_{o}}{\|\mathbf{Z}/\sqrt{T}\|_{o}^{2}\widetilde{\phi}(s_{0}+M_{2}\widetilde{s}_{0})}$. By using this result, we have that on $\Theta_{\widetilde{s}_{0}}$:

[TABLE]

where $\widetilde{\eta}:=\max\{\eta,1\}$ and $d(f_{0},f)$ is given in (A11).

Define $\mathcal{A}:=\{(\boldsymbol{\theta},\sigma^{2})\in\mathbb{R}^{\widetilde{g}}\times\mathbb{R}_{+};d(f_{0},f)\leq M_{1}\epsilon^{2}\}$. Hence, from (A16), (A17) and the Fubini’s theorem, we obtain:

[TABLE]

By putting together (A15) and (A18) we obtain:

[TABLE]

which converges to zero uniformly over $(\boldsymbol{\theta}_{0},\sigma_{0}^{2})\in\overline{\Theta}_{0}\times[\underline{\sigma}^{2},\overline{\sigma}^{2}]$ under the condition of the theorem and because $\mathbf{E}_{0}\Pi(\mathcal{A}^{c}|y,\mathbf{Z})\rightarrow 0$ by (A10) and $\mathbf{E}_{0}\Pi(\Theta_{\widetilde{s}_{0}}^{c}|y,\mathbf{Z})\rightarrow 0$ by Lemma B.1.4 (uniformly over $(\boldsymbol{\theta}_{0},\sigma_{0}^{2})\in\overline{\Theta}_{0}\times[\underline{\sigma}^{2},\overline{\sigma}^{2}]$).

A.3 Stabilization algorithm

The stabilization algorithm used in the paper is a slightly modified version of the algorithm proposed in Andrieu et al. (2005) and discussed in Atchadé (2011), which uses re-projections on randomly varying compact sets. Recall that the updating (approximate EM) algorithm described in Section 5.2 is:

[TABLE]

where we use the transformation $\boldsymbol{\omega}=0.5\log(\boldsymbol{\lambda})$. Let $\{a^{(s)},s\geq 0\}$ and $\{e^{(s)},s\geq 0\}$ be two monotone non-increasing sequences of positive numbers. Here we choose $a^{(s)}=1/s^{q}$, with $q=0.8$, and $e^{(s)}=\overline{e}+(1-\overline{e})(1-\varsigma_{s}^{-\alpha_{e}})$, with $\overline{e}=3$ and $\alpha_{e}=0.1$. Let $\{\mathsf{K}^{(s)},s\geq 0\}$ be a monotone increasing sequence of compact subsets of $\boldsymbol{\Omega}$ such that $\bigcup_{s\geq 0}\mathsf{K}^{(s)}=\boldsymbol{\Omega}$. Here we set compact subsets of the form $\mathsf{K}^{(s)}=[\max(-\kappa_{s}-1,-c),\kappa_{s}+1]$, where $c>0$. To avoid unstable outcomes due to extremely small numbers in $\boldsymbol{\lambda}$, we set $c=5$. Let $\widetilde{\boldsymbol{\Omega}}\times\widetilde{\boldsymbol{\Phi}}\subset\mathsf{K}^{(s)}\times\boldsymbol{\Phi}$ and $\Pi:\boldsymbol{\Omega}\times\boldsymbol{\Phi}\rightarrow\widetilde{\boldsymbol{\Omega}}\times\widetilde{\boldsymbol{\Phi}}$ be a re-projection function, such as $\widetilde{\boldsymbol{\Omega}}\times\widetilde{\boldsymbol{\Phi}}=(\widetilde{\boldsymbol{\phi}},\widetilde{\boldsymbol{\omega}})$ for an arbitrary point $(\widetilde{\boldsymbol{\omega}},\widetilde{\boldsymbol{\phi}})\in\mathsf{K}^{(s)}\times\widetilde{\boldsymbol{\Phi}}$. Let $\varphi$ be a function such that $\varphi(w)=1-w$, for all $w\geq 0$.

With this algorithm, $\kappa_{s}$ is the index of the active truncation set (also equal to the number of restarts before $s$), $\nu_{s}$ is the number of iterations since the last restart, and $\varsigma_{s}$ is the current index in the step-size sequence. We set $\varphi(w)=1$ for all $w\in\mathbb{N}$, such that $\varsigma_{s}=s$. Hence, if $\overline{\boldsymbol{\omega}}\notin\mathsf{K}^{(\kappa_{s-1})}$ or $|\overline{\boldsymbol{\omega}}-\boldsymbol{\omega}^{(s-1)}|>e^{(\varsigma_{s-1})}$, we re-initialize the algorithm starting from $(\widetilde{\boldsymbol{\omega}},\widetilde{\boldsymbol{\phi}})$, which are obtained by drawing from:

[TABLE]

where $\mathsf{K}^{(s-1)}_{u}=\kappa_{s-1}+1$ and $\mathsf{K}^{(s-1)}_{l}=\max(-\kappa_{s-1}-1,-c)$, and parameters $\boldsymbol{\phi}|\widetilde{\boldsymbol{\omega}}$ are drawn from the prior distributions described in Sections 3.1 and 3.2. We then iterate until the acceptance conditions stated in the algorithm are met. Finally, we set the new compact subsets to $\mathsf{K}^{(\kappa_{s-1}+1)}$ and the new sequence of step-size.