Bayesian Bi-level Sparse Group Regressions for Macroeconomic Density Forecasting

singlespace Bayesian Bi-level Sparse Group Regressions for Macroeconomic Density Forecasting singlespace

\setstcolor{red}

\newgeometry{lmargin=1.20in,rmargin=1.20in,tmargin=1.20in,bmargin=1.20in}

\setcounter{page}{2} \setcounter{footnote}{0}

Introduction

Density forecasts of macroeconomic and financial time-series are of crucial importance for policymakers, as they provide an assessment of future (upside or downside) tail risks GarrattEtAl2003,AdrianEtAl2019,Adams2021,Carriero2024. Nowadays, researchers dispose of rich datasets, released by official or alternative sources, that potentially contain valuable information for predicting tail events. On the one hand, handling these datasets might be challenging, mainly because of their large dimension (i.e. the number of variables is large compared to the number of available observations in the time dimension). On the other hand, the variables in these datasets are (or can be) often organised in groups, a feature that may help the researcher in reducing the dimensionality. The research question of interest is then: how can we optimally exploit group-structures to construct accurate density forecast? To the best of our knowledge, there is no clear answer to this question in the literature for general macroeconomic forecasting models. The contribution of this paper is to fill this gap by constructing optimal density forecasts that take advantage of the group-structure for simultaneously i) reducing the potentially large dimension in the dataset, and ii) evaluating the tail risks probabilities. By detecting, in a data-driven way, the relevant driving factors and groups, our procedure leads to economically interpretable forecasting models and results, which is a crucial feature for decision-taking by policymakers.

Groups might be formed by predictors that present strong covariation, as well as common characteristics and patterns. Some datasets may be organised in groups by the researcher according to economic criteria. For example, real economic activity can be analysed by using a large number of official series arranged in homogeneous blocks of indicators, such as production, employment, consumption, and housing (see e.g. McCracken2016). Alternatively, datasets may be organised in a group-structure directly by the data provider. This is the case, for instance, of the Google Search data used in FerraraSimoniJBES. In this paper, we widen the relevance of group-structures beyond datasets organised in groups, and we point out that specific forecasting models -- like models with unknown nonlinearities or with covariates observed at mixed frequencies -- might present a group-structure.

To allow for different sources of group-structure, we consider a general model flexible enough to include, among others, i) linear regression models with many predictors organised in groups, ii) mixed-frequency regression models with nonparametric weighting functions, and iii) nonlinear forecasting models as special cases. The nonparametric models ii) and iii) are suitable for capturing complex relationships between the target variable and the predictors, and the group-structure arises there from generalized Fourier expansions of these relationships.

While the true model is allowed to have an infinite number of elements in each group, we approximate the model by restricting each group to have no more than $g\geq 1$ components. This introduces a bias in finite samples, but this bias vanishes asymptotically as $g$ is allowed to increase with the sample size $T$. The approximate model is assumed to be bi-level sparse: only $s_0^{gr}$ among the $N$ groups are active, and within each of the $s_0^{gr}$ active groups, only a few elements have a non-zero impact on the target variable. This means that some groups and some predictors within an active group can be irrelevant for modelling and forecasting the target variable, conditional on the remaining predictors.

Our Bayesian procedure is based on a hierarchical prior that generates a group structure with bi-level sparsity and charges only the approximate model. To induce exact bi-level sparsity, we use $(g+1)N$ spike-and-slab priors specified as mixtures between a continuous distribution and a Dirac distribution at zero, the latter inducing exact zeros with non-negligible probability. Density forecasts are obtained from the posterior predictive distribution, which is well-known to dominate plug-in predictive distributions when there is prior information. We establish frequentist asymptotic optimality of our Bayesian procedure in a setting where the true data generation process (DGP) has a bi-level sparse group-structure. Importantly, optimality does not require orthogonality between the groups of predictors, such that cross-correlation between covariates in different groups is allowed by our theory. Asymptotic properties are established for the sample size $T$ increasing to infinity, as well as for the number of groups $N$, active groups $s_0^{gr}$, components per groups $g$, and active predictors in the model $s_0$ increasing to infinity with $T$. As a by-product, we provide parameter recovery and point forecasts, as well as the contraction rate of the posterior distribution of the in-sample prediction error. Our posterior contraction rate attains the minimax rate established in Cai2022IEEE and LiZhangYin2024 for bi-level sparse settings.

The bi-level sparse group-structure that we consider in this paper has remarkable advantages over the sparse group-structure considered e.g. in MoglianiSimoni2021, which only accounts for sparsity among groups but not within groups. First, if the true model is bi-level sparse, then the posterior of MoglianiSimoni2021 obtained under group sparsity would not put zero mass asymptotically on models with dense groups. Therefore, the true sparsity pattern is not recovered with positive probability, implying a worsening of both selection of relevant driving factors and out-of-sample predictive accuracy. This issue is illustrated in our first Monte Carlo simulation. Second, the contraction rate with bi-level sparsity is faster than the rate with group sparsity if $g\gg \log(N)/\log(T)$ and if $s_0^{gr}g \gg \log(s_0^{gr} g) s_0$. It follows that, in this case, the required sample size for consistency would be large. This is not appealing for macroeconomic forecasting, where the time-series dimension $T$ is usually not large.

Finally, our framework can be easily extended to account for stochastic volatility and ARMA errors, which are important features in several macroeconomic applications and have been proven to improve predictive results Chan2014,Zhang2020,CCMM2024. We use this extended setting in an empirical application on density nowcasting the US quarter-on-quarter GDP growth rate.

The paper is structured as follows. In Sections (ref) and (ref) we present respectively the model and the prior. Conditional posterior and predictive distributions are provided in Section (ref). The asymptotic properties of our procedure are analysed in Section (ref), while Monte Carlo experiments are discussed in Section (ref). Section (ref) provides an empirical application. Section (ref) concludes. Additional theoretical results and all the proofs are reported in the Online Appendix.

The model

The sampling model

Let $y_t$ be the series of interest to be predicted $h$-steps ahead. At each time $t$, a large number of predictors, organized into $N$ groups and denoted by $\mbox{\bf x}_{j,t} := \{x_{j,t,i}\}_{i\geq 1}$ for every $j\in\{1,\ldots,N\}$, are available to the forecaster. A group $j$ might be formed, for instance, by indicators belonging to a given sector or category, by lagged values of one predictors or the dependent variable, or by functional transformations of one predictor. Density forecasts of $y_T$, conditional on information available up to $T-h$, are based on the model:

for every $h \geq 0$ and $t=1,\ldots,T$. Examples of model (ref) will be presented in Section (ref) and an extension to the case where one group contains lagged values of $y_t$ is considered in Appendix A.7. For every $j\in\{1,\ldots, N\}$, $\varphi_j(\cdot)$ denotes a $j$-specific unknown function of $\mbox{\bf x}_{j,t-h}$, belonging to a separable Hilbert space $\mathcal{H}_j$ and taking values in $\mathbb{R}$. Depending on the dimension of $\mbox{\bf x}_{j,t-h}$, each function $\varphi_j$ might depend on a potentially infinite number of arguments. We assume that $\varepsilon_t$ are independent and identically distributed according to a $\mathcal{N}(0,\sigma^2)$ distribution. Therefore, by introducing the vector $\mbox{\bf x}_t := (\mbox{\bf x}_{1,t}', \ldots,\mbox{\bf x}_{N,t}')'$ of potentially infinite dimension, the matrix $\mbox{\bf X} := (\mbox{\bf x}_{1-h}, \ldots, \mbox{\bf x}_{T-h})'$ with $T$ rows, and the $N$-vectors $\varphi(\mbox{\bf x}_t):= (\varphi_1(\mbox{\bf x}_{1,t}), \ldots, \varphi_N(\mbox{\bf x}_{N,t}))'$ and $\varphi:=(\varphi_1,\ldots,\varphi_N)'$, we write the sampling model as:

and the joint sampling distribution of $y := (y_{1},\ldots, y_{T})'$ conditional on $\mbox{\bf X}$ is $\prod_{t=1}^T \\\mathcal{N}\left(\sum_{j=1}^N \varphi_{j}(\mbox{\bf x}_{j,t-h}),\sigma^2\right)$.\\ For $j=1,\ldots,N$, let $\{z_{j,t-h,i}\}_{i\geq 1}$ be transformations of the elements of $\mbox{\bf x}_{j,t-h}$ such that the function $\varphi_j(\mbox{\bf x}_{j,t-h})$ writes as $\varphi_j(\mbox{\bf x}_{j,t-h}) = \sum_{i=1}^{\infty} \theta_{j,i}z_{j,t-h,i}$. In the remainder of this paper, we shall use this expression for the function $\varphi_j(\cdot)$. To reduce the dimension of the model, each function $\varphi_j(\mbox{\bf x}_{j,t-h})$ is then approximated by $\sum_{i=1}^g \theta_{j,i} z_{j,t-h,i}$, where $g \geq 1$ is a truncation parameter. We introduce the following notation associated with this approximation: for every $j=1,\ldots,N$, define $\boldsymbol{\theta}_j := (\theta_{j,1},\ldots,\theta_{j,g})' \in \mathbb{R}^g$, $\boldsymbol{\theta} := (\boldsymbol{\theta}_1', \ldots, \boldsymbol{\theta}_N')'\in\Theta \subset \mathbb{R}^{Ng}$, $\mbox{\bf z}_{j,t-h}:= (z_{j,t-h,1},\ldots,z_{j,t-h,g})'$, $\mbox{\bf z}_{t} := (\mbox{\bf z}_{1,t}^{'},\ldots,\mbox{\bf z}_{N,t}^{'})'$ a $(Ng\times 1)$ vector, and $\mbox{\bf Z} := (\mbox{\bf z}_{1-h},\ldots, \mbox{\bf z}_{T-h})'$ a $(T \times Ng)$ matrix. By using this notation, the approximation bias in the mean is:

and $B(g):= (B_1(g), \ldots, B_T(g))'$ is a $T$-vector. Therefore, $B_{t,j}(g) := \sum_{i> g} \theta_{j,i}z_{j,t-h,i}$. In this paper we adopt a Bayesian approach and specify a convenient prior for $\boldsymbol{\theta}$ such that the induced prior for $B_t(g)$ is degenerate at zero (see Section (ref)).

The bi-level sparsity structure

Let $(\varphi_{0},\sigma_0^2)$ be the true value of $(\varphi,\sigma^2)$ that generates the data. Under the gaussianity assumption of the error term, the true conditional distribution of $y$ given $\mbox{\bf X}$ is $\prod_{t=1}^T\mathcal{N}\left(\sum_{j=1}^N \varphi_{0,j},\sigma_0^2\right)$, with Lebesgue density denoted by $f_0$. The approximation bias of the true sampling mean is denoted by $B_0(g)= (B_{0,1}(g),\ldots, B_{0,T}(g))$ and the true value of $\boldsymbol{\theta}$ in the approximation of $\varphi_0(\mbox{\bf x}_t)$ by $\boldsymbol{\theta}_{0}:=(\boldsymbol{\theta}_{0,1}',\ldots, \boldsymbol{\theta}_{0,N}')'$. The latter is assumed to be bi-level group sparse in the sense explained below.\\ Exact bi-level group sparsity is the feature of the model that guarantees the existence of an approximation $\mbox{\bf z}_{t-h}'\boldsymbol{\theta}_0 \equiv \sum_{j=1}^N \mbox{\bf z}_{j,t-h}'\boldsymbol{\theta}_{0,j}$ to $\sum_{j=1}^N \varphi_{0,j}(\mbox{\bf x}_{j,t-h})$ in (ref), with a small number of active groups and of non-zero coefficients for each active group such that the approximation bias $B_{0,t}(g)$ is small relative to the estimation error. To the best of our knowledge, the assumption of exact bi-level sparsity dates back to the recent contributions of Huang2012, Simon2013, and Breheny2015, which propose algorithms for estimation and sparse recovery in cross-section models. Our approach is instead designed for time-series models, with an explicit predictive density motivation and a theoretical foundation.\\ To clarify the concept of exact bi-level sparsity of $\boldsymbol{\theta}_0$, let $S_{0}^{gr} := \{1 \leq j \leq N; \|\boldsymbol{\theta}_{0,j}\|_2 > 0\} \subset \{1,2,\ldots,N\}$ be the set of indices of the groups (subvectors) in $\boldsymbol{\theta}_0$ with at least one nonzero component (active groups), and let $s_{0}^{gr}:=|S_{0}^{gr}|$ denote the cardinality of $S_{0}^{gr}$. If $S_{0}^{gr} \neq \varnothing$, then for every $j \in S_{0}^{gr}$ let $S_{0,j}$ be the set of the indices of the nonzero elements in $\boldsymbol{\theta}_{0,j}$, such that $S_0 = \bigcup_{j\in S_0^{gr}}|S_{0,j}|$ and $s_{0} := |S_0| = \sum_{j\in S_{0}^{gr}}|S_{0,j}|$. Remark that $s_{0} \geq 1$ since there is at least one active group under the assumption $S_{0}^{gr} \neq \varnothing$. If $S_{0}^{gr} = \varnothing$, then $s_{0} = 0$. Moreover, $s_0 \equiv s_0(g)$ is a non-decreasing function of $g$. Hence, we say that $\boldsymbol{\theta}_0$ is $(s_0,s_0^{gr})$-sparse.\\ The next assumption supposes exact bi-level sparsity of $\boldsymbol{\theta}_0$ and controls the approximation bias $B_{0,t}(g)$ relative to the estimation error. The latter is similar to ChernozhukovBelloniHansen2014.

In Section (ref), we will let $N$, $g$, $s_0$ and $s_0^{gr}$ to increase with $T$. This, together with Assumption (ref), will allow the size of the approximation model to grow with the sample size $T$. The constant value in the denominator of the upper bound of the approximation bias can be replaced by any constant larger or equal than $16$.

Examples

Linear regression models with grouped predictors

Many datasets used for macroeconomic forecasting display a large panel of real and financial indicators. Theses series can be often organized in homogeneous groups by following either economic information or statistical procedures. For example, a nominal group for different price indicators, an output group for supply and production indicators, a financial group for interest rates and stock prices, etc. In this paper, we propose a procedure for density forecasting that takes into account the group structure in large datasets.

We suppose that each group of covariates contains at most $g$ elements. Let $y_t$ be the target variable to be predicted $h$-steps ahead, $j\in\{1,\ldots,N\}$ be the group index, and $\mbox{\bf x}_{j,t}$ be the $g$-vector of variables in each group $j$. Then, by assuming a linear model: $\forall t=1,\ldots,T$ and $h>0$,

which can be cast in model (ref) with $\mbox{\bf z}_{j,t} = \mbox{\bf x}_{j,t}$, $\varphi_j(\mbox{\bf x}_{j,t-h}) = \mbox{\bf x}_{j,t-h}'\boldsymbol{\theta}_j$ and zero approximation bias.

Nonparametric mixed-frequency regression models

Consider a low-frequency variable $y_{t}^{L}$, where $t=1,\ldots,T$ indexes the low frequency time unit, and consider $(N-1)$ high-frequency variables $x_{j,t}^H$ for $j=1,\ldots, N-1$. By denoting with $m$ the number of times the higher sampling frequency appears in the low-frequency time unit $t$, then $t - k/m$ denotes the $k$-th past high frequency period for $k = 0,1,2,\ldots$. To simplify the notation, we set the same $m$ for all the groups, but our framework can accommodate the case with a $m_j$ specific to each group. 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_{j,t}^{H}=x_{j,t-\scalebox{0.55}{$\left.1\middle/m\right.$}}^{H}$. Further, let $h=0,1/m,2/m,3/m,\dots$ be an (arbitrary) forecast horizon. For given orders $p_y, p_x\geq 0$, the general Mixed Data Sampling (MIDAS) regression model can be written as follows: $\forall t=1,\ldots,T$,

where $\Psi_j(L^{1/m})$ is the high-frequency lag polynomial

This model can be cast in model (ref) with $N$ groups, where $\varphi_j(\mbox{\bf x}_{j,t-h}) = \Psi_j(L^{1/m})x_{j,t-h}^{H}$ for every $j=1,\ldots,N$ with $\mbox{\bf x}_{j,t-h} = (x_{j,t-h}^H, \ldots, x_{j,t-h - p_x/m}^H)^{\prime}$.\footnote{Model (ref) can be generalized to accommodate groups of indicators sampled at the same frequency as the target variable $y_{t}$.} Previous literature has considered different parameterizations of the weighting function $\psi_j(u)$ in (ref). Parametric specifications have been considered, for instance, by Foroni2015 and Ghysels2007. MoglianiSimoni2021 use (non-orthogonalized) algebraic power polynomials which correspond to restricted Almon lag polynomials, while Babii2022 use shifted orthogonal polynomials such as Jacobi and Legendre. With the method developed in the present paper, we can allow $\psi_j(\cdot)$ to belong to a separable Hilbert space $\mathcal{H}_j$ with a countable orthonormal basis $\{\phi_1(\cdot), \phi_2(\cdot),\ldots\}$ and construct density forecasts for MIDAS models. Therefore, for any $\psi_j(\cdot)\in\mathcal{H}_j$, we can write

where $\theta_{ji} := \langle \psi_j,\phi_i\rangle$ is the $i$-th Fourier coefficient. Hence,

where $\Phi_i := (\phi_i(0),\phi_i(1), \ldots, \phi_i(p_x))^{\prime}$. By cutting the sum in $i$ at some $g > 0$, we get $ \varphi_j(\mbox{\bf x}_{j, t-h}) = \Psi_j(L^{1/m})x_{j,t-h}^H = \sum_{i=1}^{g}\theta_{ji}\Phi_i^{\prime} \mbox{\bf x}_{j, t-h} + B_{t,j}(g)$, which yields a mixed-frequency regression model with approximated high-frequency polynomial $\sum_{i=1}^{g}\theta_{ji}\Phi_i^{\prime} \mbox{\bf x}_{j, t-h}$, $z_{j,t-h,i} = \Phi_i^{\prime} \mbox{\bf x}_{j, t-h}$, $\mbox{\bf z}_{j,t-h} = (\mbox{\bf x}_{j, t-h}^{\prime}\Phi_1, \ldots, \mbox{\bf x}_{j, t-h}^{\prime}\Phi_g)'$, and an approximation bias at time $t$ given by $B_{t,j}(g) = \sum_{i > g}\theta_{ji}\Phi_i^{\prime} \mbox{\bf x}_{j, t-h}$.

Nonlinear predictive regression models with interaction effects

In many settings, it may be desirable to account for possible nonlinearities in the conditional mean function of the target variable $y_t$. Model (ref) can accommodate different types of nonlinearities. The simplest case is when, given $p$ covariates, the effect of each covariate on $y_t$ can be separated in $p$ nonlinear functions and interaction effects are taken into account by using an additive partially linear model: $\forall t=1,\ldots,T$ and $h>0$,

where $z_{t-h} := \{x_{j,t-h}x_{k,t-h}\}_{k>j}$ is a $p(p-1)/2$-vector that includes the interactions between covariates. Here, we have $N=p+1$ groups, with the last group having a very high number of components (of the order $p^2$), and $\varphi_j$ is a function of only one covariate taking values in a separable Hilbert space $\mathcal{H}_j$ with a countable orthonormal basis $\{\phi_{j,1}(\cdot), \phi_{j,2}(\cdot),\ldots\}$: $\varphi_j\in\mathcal{H}_j$. Then, for every $j\in\{1,\ldots,p\}$, $\varphi_j(x_{j,t-h}) = \sum_{i=1}^{\infty}\theta_{ji}\phi_{j,i}(x_{j,t-h})$, where $\theta_{ji} := \langle \varphi_j,\phi_{j,i}\rangle$ is the $i$-th Fourier coefficient and $\langle \cdot,\cdot\rangle$ denotes the inner product in $\mathcal{H}_j$. For some $g>0$ we can approximate $\varphi_j(x_{j,t-h})$ as

which yields an approximation bias at time $t$ given by $B_{t,j}(g) = \sum_{i > g}\theta_{ji}\phi_{ji}(x_{j,t-h})$. Alternatively, we can replace $z_{t-h}'\boldsymbol{\theta}_{p+1}$ with $\sum_{j=1}^p z_{j,t-h}'\boldsymbol{\theta}_{p+j}$ in (ref), where $z_{j,t-h} := \{x_{j,t-h}x_{k,t-h}\}_{k\neq j}$ is a $(p-1)$-vector. The total number of groups would hence be $N=2p$, with the last $p$ groups of dimension $p-1$ each, which is much lower than $p^2$. It follows that the choice of the group structure may imply a trade-off between the number of groups and the number of components of each group.

The BSGS-SS prior

With the approximation model in Section (ref) and the Assumption (ref) in mind, we elicit a prior that puts all its mass on the approximation $\mbox{\bf z}_{t-h}'\boldsymbol{\theta}$, conditional on $\mbox{\bf z}_{t-h}$, and that induces sparsity at the group level and within groups. Let us first define, for every group $j=1,\ldots,N$, the following quantities: $\boldsymbol{\theta}_{j} = \mathbf{V}_j^{1/2}\mbox{\bf b}_{j}$, $\mbox{\bf b}_{j} := (b_{j1},\ldots,b_{jg})'$, $\mathbf{V}_j^{1/2} := \mathrm{diag}(v_{j1},\ldots, v_{j g})$ and $v_{ji} \geq 0$ for $i=1,\ldots, g$.\\ We treat the truncation parameter $g$ as deterministic and, under Assumption (ref), it may depend on $s_0$. Our proposed Bayesian Sparse Group Selection with Spike-and-Slab prior (BSGS-SS henceforth) is specified as: $\forall j=1,\ldots, N $,

where $\mathcal{N}^+(0,\tau_{j}^{2})$ denotes a $\mathcal{N}(0,\tau_{j}^{2})$ distribution truncated below at zero, $\delta_0(\cdot)$ denotes a Dirac distribution at zero, $\mathcal{G}$ (resp. $\mathcal{G}^{-1}$) denotes the Gamma (resp. Inverse-Gamma) distribution, $\lambda_{1,j}$ in the prior for $\tau_j$ is the scale parameter, and $\mathcal{B}$ denotes a Beta distribution. The conditional priors for $\mbox{\bf b}_{j}|g,\pi_0$ and for $v_{ji}|\pi_1,\tau_j^2$ are specified as hard spike-and-slab priors. A similar prior has been considered in Xu2015. However, a subtle difference here stems from the non-conjugate Gamma prior for the standard deviation hyperparameter $\tau_j$. We prove in Section (ref) that this specification is crucial to ensure good asymptotic properties for the corresponding posterior. Compared with the priors in ChenChuYuanWu2016, the proposed prior is computationally convenient as we can simulate from the posterior distribution by using a relatively simple Gibbs sampler with one Metropolis-Hastings step.\\ By integrating out $\{v_{ji}\}_{i=1}^g$, the induced prior on $\theta_{ji}|\pi_0,\pi_1,\tau_j$ is: $\forall j=1,\ldots,N$,

where the Lebesgue density $f_{\theta_{ji}}$ is upper bounded by the Lebesgue density of a $\mathcal{G}(1/2,\tau_j)$ distribution. Its explicit expression is given in the Online Appendix A.1, together with a proof of this result. In the appendix we also discuss the conditional prior on the random function $\varphi_j$, given $\{v_{ji}\}_{i=1}^g$ and $g$, induced by the prior (ref)-(ref).\\

The BSGS-SS predictive density and the MCMC algorithm

Let $\Pi(\cdot|y,\mathbf{X})$ denote the posterior of the model's parameters $(\varphi,\sigma^2)$ obtained from the sampling model (ref) and the BSGS-SS prior in Section (ref). We construct density forecasts based on the posterior predictive density of $y_{\tau}|\mbox{\bf x}_{\tau-h}$, denoted by $\widehat{f}(y_{\tau}|\mbox{\bf x}_{\tau - h},y,\mbox{\bf X})$ and conditional on a new value of the covariates $\mbox{\bf x}_{\tau - h}$, for $\tau>T$:

where $f_0(y_{\tau}|\mbox{\bf x}_{\tau - h},\varphi,\sigma^2) $ denotes the Lebesgue density associated with the sampling distribution in (ref). Draws from the predictive distribution (ref) can be obtained directly from the MCMC sampler detailed in Algorithm (ref), where we denote by “rest” the conditioning parameters and data in each conditional distribution.

The hyperparameters of the priors on $\pi_{0}$ and $\pi_{1}$ must be chosen very carefully as they control the overall amount of between-group and within-group prior sparsity. In this paper, we propose to select $c_0$ and $c_1$ by using a data-driven approach based on the Deviance Information Criterion (DIC; Spiegelhalter2002), defined as:

where $(\boldsymbol{\theta}^{\prime}_{\mathbf{c}},\sigma^{2}_{\mathbf{c}})^{\prime}$ denotes the model parameters for a given set of hyperparameters $\mathbf{c}=(c_{0},c_{1})$, and $(\widehat{\boldsymbol{\theta}}_{\mathbf{c}},\widehat\sigma_{\mathbf{c}}^2)$ is a point estimate of $(\boldsymbol{\theta}^{\prime}_{\mathbf{c}},\sigma^{2}_{\mathbf{c}})^{\prime}$.\footnote{The first part of (ref) is the posterior mean deviance, which is here estimated by averaging the log-likelihood function, $\log f(y\vert \boldsymbol\boldsymbol{\theta},\sigma^2,\mbox{\bf X})$, over the posterior draws of $(\boldsymbol{\theta}^{\prime}_{\mathbf{c}},\sigma^{2}_{\mathbf{c}})^{\prime}$. $(\widehat{\boldsymbol{\theta}}_{\mathbf{c}},\widehat\sigma_{\mathbf{c}}^2)$ is computed using the posterior median for $\boldsymbol{\theta}$ and the posterior mean for $\sigma^2$.} This metrics of fit presents the advantage of being easily available from the MCMC output, with no additional costs in terms of computational burden.\footnote{We select the $c_{0}$ and $c_{1}$ that minimize the DIC by performing a random search on a fine 2-dimensional grid, with lower-bounds set as described in (A.1.3). In empirical applications, $s_{0}^{gr}$ and $s_{0}$ are unknown, but we recommend to tune them at some arbitrary values such that the inequalities in (A.1.3) are satisfied for low values of $c_{0}$ and $c_{1}$. Therefore, the price to pay for not knowing $s_{0}^{gr}$ and $s_{0}$ is simply the burden of searching on a larger grid.}

Extending the model to a general error process

One of the main advantage of the proposed prior is its flexibility on the treatment of the error process. Unlike other shrinkage priors (Group Lasso, Sparse Group Lasso), it is quite straightforward to extend the prior in Section (ref) to account, for instance, for stochastic volatility and ARMA errors. The relevance of these features in macroeconomic applications has been extensively stressed by the literature (e.g., Clark2011,Chan2014,Carriero2019). Further, recent works CCMM2024,LenzaPrimiceri2022 have suggested that time-varying volatility could attenuate inferential issues (e.g., explosive patterns and widening uncertainty for forecasts and Impulse Response Functions) when the estimation sample includes extreme observations such as those recorded during the COVID-19 pandemic.\\ We then modify the homoskedastic model in (ref) to account for a general error process as follows:

where we assume a stationary and invertible ARMA($p$,$q$) specification for the errors $\varepsilon_{t}$ and heavy-tails innovations $u_t$, the latter featuring time-varying log-volatility $\zeta_{t}$ and outliers component $\omega_{t}$ in the variance. Since the Student-$t$ distribution $t_{\nu}$ in (ref) can be expressed as a scale mixture of Normal distributions with scale term $\tau_{t} \sim \mathcal{G}^{-1}(0.5\nu,0.5\nu)$, then, conditional on $(\tau_{t}, \omega_{t}, \zeta_{t})$, the innovations have distribution $u_{t} \sim \mathcal{N}\left(0,\tau_{t}\omega_{t}\exp(\zeta_{t})\right)$. The log-volatility $\zeta_{t}$ is assumed to evolve according to a stationary AR(1) process: $\zeta_{t} \sim \mathcal{N}\left(\mu_{\zeta} + \phi_{\zeta}(\zeta_{t-1}-\mu_{\zeta}),\sigma^{2}_{\zeta}\right),$ with initial conditions $\zeta_{0} \sim \mathcal{N}\left(\mu_{\zeta},\sigma^{2}_{\zeta}/(1-\phi_{\zeta}^{2})\right)$ and $\phi_{\zeta}$ restricted to lie in the stationary region. While the scale term $\tau_{t}$ is intended to capture frequent, although limited, jumps in the volatility, the additional scale term $\omega_{t}$ is instead intended to capture infrequent, but extreme, jumps. Following Stock2016 and CCMM2024, the prior for $\omega_{t}$ has a two-part distribution, depending on whether an observation at period $t$ is considered either regular or outlier: $\omega_t =1$ with probability $(1-p_{\omega})$ and $\omega_t \sim \mathcal{U}(2,\bar{\omega})$ with probability $p_{\omega}$ for a constant $\bar{\omega} > 2$, where $\mathcal{U}(2,\bar{\omega})$ denotes the Uniform distribution with support $(2,\bar{\omega})$.

The implementation of these features implies only minor modifications to the MCMC Algorithm (ref) (see the Online Appendix A.3 and Section (ref) for an application). Conditional posteriors for stochastic volatility and ARMA parameters are the same as in Kastner2014, Chan2013, and Zhang2020, and we refer the reader to these contributions for more details.

Theoretical properties

This section provides the theoretical validation of our procedure. In Section (ref) we establish consistency of the predictive density for density forecast. As a by-product, in Section (ref) we demonstrate consistency of the posterior distribution of the model's parameters $(\varphi,\sigma^2)$, denoted by $\Pi(\cdot|y,\mathbf{X})$. We adopt a frequentist point of view, in the sense that we admit the existence of a true value $(\varphi_{0},\sigma_0^2)$ that generates the data. We denote by $\mathbb{E}_0[\cdot]$ the expectation taken with respect to the true data distribution $P_0 = \prod_{t=1}^T \mathcal{N} \left(\sum_{j=1}^N \varphi_{0,j}(\mbox{\bf x}_{j,t-h}),\sigma_0^2\right)$ of $y$ conditional on $(\varphi_{0},\sigma_0^2)$ and $\mbox{\bf X}$ with Lebesgue density denoted by $f_0$. All the analysis is conditional on $\mbox{\bf Z}$. In the following, we denote by $\|\mbox{\bf Z}\|_{op}$ the operator norm and $\|\mbox{\bf Z}\|_o:=\max\{\|\mbox{\bf Z}_j\|_{op}; 1\leq j\leq N\}$, where $\mbox{\bf Z}_j$ is the $(T \times g)$-submatrix of $\mbox{\bf Z}$ made of all the rows and the columns corresponding to the indices in the $j$-th group. Moreover, $\|\cdot\|_2$ denotes the Euclidean norm. Denote $$\epsilon := \max\{\sqrt{s_0^{gr} \log(N)/T},\sqrt{s_0\log(T)/T},\sqrt{s_0\log(s_0^{gr}g)/T}\}$$ the rate of contraction of the posterior distribution of $\varphi$. The first term corresponds to the complexity of identifying $s_0^{gr}$ non-zero groups, while the third term corresponds to the complexity of estimating $s_0$ non-zero elements of a parameter distributed in $s_0^{gr}$ known groups. In absence of group-structure, the maximum between these two rates corresponds to the rate for recovery of sparse vectors over $\ell_0$-balls (see e.g. RaskuttiWainwright2011). This case can be obtained by setting either only one group (i.e. $N=1$ and $g = p$) or a number of groups equal to the number of $p$ parameters (i.e. $N=p$ and $g=1$).

For two sequences $a_T$ and $b_T$ and two constants $c_1,c_2>0$ we write $a_T \asymp b_T$ if $c_1 b_T\leq a_T \leq c_2 b_T$. By using this notation, if $\log(s_0^{gr} g) \asymp \log(s_0^{gr} g/s_0)$, $\log(N) \asymp \log(N/s_0^{gr})$, and $s_0\log(T) \leq \min\{s_0^{gr}\log(N),s_0\log(s_0^{gr} g)\}$, then $\epsilon$ corresponds to the minimax rate for recovering $\varphi$ given in Cai2022IEEE and in LiZhangYin2024. The conditions $\log(s_0^{gr} g) \asymp \log(s_0^{gr} g/s_0)$ and $\log(N) \asymp \log(N/s_0^{gr})$ guarantee a sparse setting, while the condition $s_0\log(T) \leq \min\{s_0^{gr}\log(N),s_0\log(s_0^{gr} g)\}$ guarantees a high-dimensional setting. If $s_0 = s_0^{gr} g$, which corresponds to the case with no sparsity within the groups, then the last condition corresponds to Assumption 4.1(ii) in MoglianiSimoni2021.\\ The following assumption restricts some of the parameters of the model.

Assumption (ref) (i) excludes degenerate cases by restricting the model variance. Assumption (ref) (ii) restricts the rate at which the number of groups can increase compared to $T$, $s_0^{gr}$ and $g$. It is violated for instance if $N>\max\{T,\exp\{s_0^{gr} g\}\}$. Assumption (ref) (iii) restricts the growth rate of the largest component of $\boldsymbol{\theta}_0$.

The next two assumptions concern the hyperparameters of the prior. Assumption (ref) allows $\lambda_{1,j}$ (the hyperparameter of the prior for $\tau_j$) to increase or decrease with $T$, $g$, $s_0^{gr}$ and $|S_{0,j}|$. It rules out a $\lambda_{1,j}$ decreasing too fast to zero or increasing too fast to infinity. In practice, any positive constant can be taken as a value for $\lambda_{1,j}$ and its choice depends on the desired tightness of the prior for $\tau_j$.

For the following assumption we introduce the function $(u)_+ := \max\{u,0\}$.

Assumption (ref) (i) and (ii) require that $c_0$ and $c_1$ increase together with $N$, $s_0^{gr}$ and $g$ and control their rate. To satisfy the assumption, if $d_0$ and $d_1$ are chosen equal to one and if $\kappa_0$ and $\kappa_1$ are chosen such that $\kappa_0 < N^{u_0}$ and $\kappa_1 < (s_0^{gr}g)^{u_1}$, then $c_0$ and $c_1$ have to be at least of the order of $1 - N +N^{u_0+1}/\kappa_0$ and $(s_0^{gr} g)^{u_1} Ng/\kappa_1 - Ng + 1$, respectively, for $u_0$, $u_1$ in the range of values given in the assumption and up to a constant. On the other hand, if $d_0$ and $d_1$ are chosen to increase with $T$, then $c_0$ and $c_1$ have to increase faster than $d_0N^{u_0 + 1}$ and $d_1(s_0^{gr} g)^{u_1}$. In practice, one can set the constants $\kappa_0$ and $\kappa_1$ at very small values, as long as they are fixed and do not increase with $T$. Assumptions (ref) (i) correspond to those in Castillo2015 for the special case of a Beta prior. We provide in Appendix A.2 a deeper analysis of this assumption.

Consistency of the predictive distribution

We now provide a validation of our procedure for density forecasts. The next theorem states that, for every $\tau>T$, $\widehat{f}(y_{\tau}|\mbox{\bf x}_{\tau - h},y,\mbox{\bf X})$ in (ref) converges to the Lebesgue density $f_0(y_{\tau}|\mbox{\bf x}_{\tau - h},\varphi_0,\sigma_0^2)$ of the true distribution $\mathcal{N} \left(\sum_{j=1}^N \varphi_{0,j}(\mbox{\bf x}_{j,\tau-h}),\sigma_0^2\right)$ with respect to the Hellinger distance denoted by $d_H(\cdot,\cdot)$.

We know from (A.6.18) and Theorem A.5.1 in the Online Appendix that the probabilities in (ref) go to zero. The assumption of the theorem is then just an assumption about the rate at which these posterior probabilities go to zero. We refer to HoffmannRousseauSH2015 for primitive conditions for this assumption.

Consistency of the posterior distribution

We establish here the posterior consistency uniformly over the infinite dimensional set $\mathcal{F}(s_0,s_0^{gr};\mbox{\bf Z})$ defined as:

for $s_0^{gr}, s_0 \in \mathbb{N}_+$ satisfying $s_0^{gr} \leq N$ and $s_0^{gr} \leq s_0 \leq g\, s_0^{gr}$. The next theorem establishes the in-sample consistency of our BSGS-SS procedure. It shows that the posterior contracts in a neighborhood of the true value defined by the in-sample prediction error. We use the notation $\varphi_{j}^{(T)}(\mbox{\bf X}) := (\varphi_j(\mbox{\bf x}_{j,-h+1}), \ldots, \varphi_j(\mbox{\bf x}_{j,T-h}))'$ for $j\in\{1,\ldots,N\}$ and $\mathcal{H} = \mathcal{H}_1\times \cdots \times \mathcal{H}_N$.

A similar result is obtained for the case with lagged values of the dependent variable among the covariates (see the Online Appendix A.7). The required sample size to achieve convergence to zero of the in-sample prediction error in Theorem (ref) is $T> C\max\{s_0^{gr}\log(N), s_0 \log(s_0^{gr} g), s_0\log(T)\}$ for some positive constant $C$. In the grouped predictors example of Section (ref), the norm in Theorem (ref) is the inner product weighted by the Gram matrix $\mbox{\bf X}'\mbox{\bf X}$, that is $\|\mbox{\bf X}(\boldsymbol{\theta} - \boldsymbol{\theta}_0)\|_2^2$. In the mixed-frequency example of Section (ref), the norm is: $$\left\|\sum_{j=1}^N \left(\varphi_j^{(T)}(\mbox{\bf X}) - \varphi_{0,j}^{(T)}(\mbox{\bf X})\right)\right\|_2^2 = \sum_{t=1}^T\left(\sum_{j=1}^N \sum_{i=1}^{\infty}\left(\theta_{j,i} - \theta_{0,ji}\right)\Phi_i^{\prime} \mbox{\bf x}_{j, t-h}\right)^2.$$

In addition to the entire function $\varphi$, it is worth focusing on the coefficients $\boldsymbol{\theta}$. The consistency of the marginal posterior of $\boldsymbol{\theta}$ is established in Theorem A.5.1 in the Online Appendix. The contraction rate for $\boldsymbol{\theta}$ coincides with the rate obtained in LiZhangYin2024 for a linear regression model.\\ We now look at the posterior mean of $\varphi(\mbox{\bf x}_t)$ as a possible point estimator for $\varphi_0(\mbox{\bf x}_t)$. Convergence of a Bayesian point estimator towards the true $\varphi_0$ is in general not implied by the result of Theorem (ref) if the loss function is not bounded. The next theorem analyzes the asymptotic behaviour of the posterior mean in terms of the Euclidean loss function $\ell_h(\widetilde\varphi,\varphi) := \frac{1}{T}\sum_{t=1}^T\left(\sum_{j=1}^N\left[\widetilde\varphi_j(\mbox{\bf x}_{t - h}) - \varphi_j(\mbox{\bf x}_{t - h})\right]\right)^2$, which is not uniformly bounded if the parameter space is not compact.

The result of the theorem holds under the assumption that the posterior of the complement of a ball centered on $\varphi_0$, with radius proportional to the contraction rate $\epsilon^2$, becomes sufficiently small for large $T$. Results of this type have been established for instance in HoffmannRousseauSH2015.

Monte Carlo experiments

We evaluate the finite sample performance of our procedure through two Monte Carlo experiments. First, we consider a DGP featuring grouped predictors, as in the example described in Section (ref). Second, we consider a DGP featuring mixed-frequency data, as in the example described in Section (ref). Simulations are based on 200 Monte Carlo iterations, each featuring 60'000 MCMC sweeps (with a burn-in sample of 10'000 sweeps and a chain-thinning parameter set at 5).

Experiment 1: DGP with grouped predictors

In this experiment, the DGP involves $Ng=\{100,300\}$ predictors (sampled at the same frequency as the target variable) structured in $N=\{5,10,20\}$ groups. The group active set has cardinality $s_{0}^{gr}=\{1,5,10\}$ and we fix $|S_{0,j}|=1$ for all simulations, such that $s_0 = s_0^{gr}$. The elements of $S_{0}^{gr}$ and $S_{0,j}$ are fixed realizations of random draws without replacement from $\{1,\dots,N\}$ and $\{1,\dots,g\}$, respectively. The number of in-sample observations is fixed at $T=200$, and the out-of-sample at $T_{oos}=50$. Simulated data are obtained from the following process:

where $\boldsymbol\epsilon_{t}:=(\epsilon_{1,t,1},\ldots,\epsilon_{N,t,g})'$, $\boldsymbol\Sigma_{\epsilon}=\mathcal{S}_{\epsilon}\mathcal{R}_{\epsilon}\mathcal{S}_{\epsilon}$ is a block-diagonal matrix, $\mathcal{S}_{\epsilon}$ a $(Ng \times Ng)$ diagonal matrix with elements $\sigma_{\epsilon}$, and $\mathcal{R}_{\epsilon}$ a block-diagonal Toeplitz correlation matrix, with $N$ blocks each of size $(g \times g)$ and featuring diagonal elements equal to one and off-diagonal elements $\rho_{\epsilon}^{|i-i^{\prime}|}$ for all $i\neq i^{\prime}$. The group structure is therefore given by the structure of this covariance matrix. We set $\sigma=0.50$ and $\rho_{\epsilon}=0.50$, and we calibrate $\sigma_{\epsilon}$ such that the noise-to-signal ratio of (ref) is $\text{NSR}=0.20$ for all the simulations. The coefficients of the active variables is set to $|\theta_{ji}|=0.5$, for each $j\in S_{0}^{gr}$ and $i\in S_{0,j}$. The sign of the coefficients is a fixed realization of random draws with replacement from $\{-1,1\}$. Finally, we set $\alpha=0.2$, $\beta_{y}=0.3$, and $\rho_{z}=0.9$.

Simulation results are reported in Table (ref), where we focus on both selection and density forecast performance. The former is evaluated by computing the Matthews correlation coefficient at the group level (MCC$_{N}$), which measures the overall quality of the groups classification, while the latter is evaluated by the means of the Continuously Ranked Probability Score (CRPS), averaged over the out-of-sample observations and expressed in relative terms with respect to an AR(1) benchmark. The results for our prior (BSGS-SS) are then compared to those obtained from the Bayesian Sparse Group Lasso (BSGL) prior of Xu2015 and the Bayesian Adaptive Group Lasso with Spike-and-Slab prior (BAGL-SS) prior of MoglianiSimoni2021. The values reported in the table denote average outcomes across Monte Carlo iterations and their bootstrap standard errors (in parentheses).

The results for our BSGS-SS prior suggest that true group active set can be correctly selected with extremely high probability (MCC$_{N}$ close to one) and point to accurate density forecasts (CRPS less than one). This holds irrespective of the number of groups $N$ and for DGPs involving a very sparse environment (i.e. when $s_{0}^{gr}$ and $s_{0}$ are small compared to $N$ and $g$). However, consistently with the contraction rates derived in Section (ref), selection and predictive accuracy tend to deteriorate in less sparse DGPs. Compared to the alternative priors considered, the results also point to a substantial outperformance of our BSGS-SS prior, in particular for highly sparse DGPs. Focusing on density forecasts, three main findings can be remarked. First, the predictive accuracy of our model is almost unaffected by the increase in the total number of predictors, while both the BSGL and the BAGL-SS priors display a fairly strong deterioration of their forecasting performance. Second, as mentioned above, our model fails to significantly outperform the AR(1) benchmark in a less sparse environment, but it still significantly outperforms the competing models. Third, the predictive performance of single-level group sparsity priors, such as the BAGL-SS, may quickly deteriorate in presence of complex sparsity structures (such as bi-level sparsity).

Additional simulation results and robustness checks are reported in the Online Appendix A.4. In particular, we show that in a very sparse environment, the mean squared error (MSE) is extremely low (with almost zero bias) and the true active set (both at group and variables level) is correctly selected with extremely high probability (according to both the True Positive Rate and the Matthews correlation coefficient). Further, we show that the performance of our model is robust to different correlation structures in $\boldsymbol\Sigma_{\epsilon}$ and alternative error distributions (e.g. Skew-Normal distribution). Overall, the results are very close to those reported in Table (ref), irrespective of the values for $Ng$, $N$, and $s_0^{gr}$, although they tend to deteriorate with a substantial increase in the noise-to-signal ratio.\footnote{We perform the following exercises: i) we increase the within-groups correlation parameter by setting $\rho_{\epsilon}=0.75$; ii) we allow for between-groups correlation by letting the Toeplitz correlation matrix $\mathcal{R}_{\epsilon}$ to be a full matrix, i.e. $\mathcal{R}_{\epsilon}=\rho_{\epsilon}^{|\iota-\iota^{\prime}|}$ for all $\iota\neq \iota^{\prime}$, with $\iota=1,\dots,Ng$; iii) we allow for asymmetric shocks in the DGP by letting $\varepsilon_{t}$ to follow a $\text{Skew-}\mathcal{N}(\xi,\omega^{2},\alpha)$ distribution Azzalini2014, with shape parameter $\alpha=-5$, and location and scale parameters $(\xi,\omega)$ calibrated such that $\varepsilon_{t}$ has mean zero and variance $\sigma^{2}$; iv) we increase the noise-to-signal ratio to 0.5.}

As a final exercise, we compare the variable selection accuracy of the BSGS-SS prior and the BAGL-SS prior of MoglianiSimoni2021. The latter relies on a single layer of sparsity (that is, sparsity at the group level only), and for this reason it is expected to be an appropriate prior only when dealing with specifications that are dense within each active group. Instead, when the true model is bi-level sparse, assuming a single-level group sparsity as in the BAGL-SS-prior, instead of bi-level group sparsity, does not ensure the posterior recovery of the correct sparsity pattern. This entails inaccurate selection and out-of-sample prediction. We investigate this issue numerically and the results are reported in Figure (ref), which shows the MCC computed at the variable level (MCC$_g$) and across various DGPs. The simulation results confirm the theoretical results and strongly suggest a substantially lower selection accuracy for the BAGL-SS (MCC$_g<1$) compared to our BSGS-SS prior (MCC$_g$ close to one) when the true model is in fact bi-level sparse.

Experiment 2: DGP with mixed-frequency data

In this example, we consider a DGP with mixed-frequency data, where $N=\{50,100\}$ predictors are sampled at a higher frequency compared to the target variable. Simulated data are obtained from the following process:

where $m=3$ (which is akin to a quarterly-monthly process) and $p_{x}=11$. We set the true weighting scheme $\psi_{j}\left(u\right)$ by relying on a three-parameters Beta function:

We investigate four alternative weighting schemes that correspond to bell-shaped weights (DGP 1) with $(a,b,c)=(5,15,0)$, fast-decaying weights (DGP 2) with $(a,b,c)=(1,10,0)$, slow-decaying weights (DGP 3) with $(a,b,c)=(1,4,0)$, and flat weights (DGP 4) with $(a,b,c)=(1,1,0)$.\footnote{The weights are set to exactly zero for values of the Beta function $<1\text{e-04}$, and normalized to sum up to 1.} Note that the same weighting scheme applies to all the predictors entering the DGPs.

We define $\varepsilon_{t}$ and $\boldsymbol\epsilon_{t} := (\epsilon_{1,t},\dots, \epsilon_{N,t})^{\prime}$ as i.i.d. draws from the normal multivariate distribution in (ref), where $\boldsymbol\Sigma_{\epsilon}=\mathcal{S}_{\epsilon}\mathcal{R}_{\epsilon}\mathcal{S}_{\epsilon}$, with $\mathcal{S}_{\epsilon}$ a diagonal matrix with elements $\sigma_{\epsilon}$ and $\mathcal{R}_{\epsilon}$ a Toeplitz correlation matrix with diagonal elements equal to one and off-diagonal elements $\rho_{\epsilon}^{|j-j^{\prime}|}$ for all $j\neq j^{\prime}$. We set $\rho_{\epsilon}=0.50$, $\alpha=0.50$, $\sigma=0.50$, $\rho_{x}=0.90$. The active set has cardinality $s_{0}^{gr}=\{5,10\}$ and the indices of active variables in $S_{0}^{gr}$ are fixed realizations of random draws without replacement from $\{1,\dots,N\}$. Conditional on these parameters, we calibrate $\sigma^2_{\epsilon}$ such that the noise-to-signal ratio of the mixed-frequency regression is $\text{NSR}=0.20$ for all the simulations. We assume $h=0$ (akin to a nowcasting model). The number of in- and out-of-sample observations is fixed at $T=200$ and $T_{oos}=50$, respectively.

In the simulations, we consider and evaluate a set of alternative approximating functions for the true weighting function $\psi_j\left(u\right)$: the U-MIDAS Foroni2015, Almon lag polynomials (restricted and unrestricted; MoglianiSimoni2021), and orthogonal lag polynomials (Legendre, Chebychev (first kind), and Bernstein orthogonal polynomials).

Selection and predictive accuracy performance of our BSGS-SS procedure are reported in Table (ref), where we report the True Positive Rate at the group level (TPR$_N$) and the CRPS. The results point to a number of interesting features. First, among the lag polynomials considered, the best results are provided by the unrestricted, the restricted Almon and the Bernstein polynomials. Second, consistently with the theory, the results for the best-performing polynomials are overall unaltered by the increase in the total number of high-frequency predictors. However, the performance is substantially affected by the increase in the number of true active variables. In this case, the restricted Almon performs surprisingly better than the other polynomials, which show, for some DGPs, very low selection and predictive accuracy. Third, the ranking of the best-performing polynomials may depend, at least in part, on the shape of the underlying true weighting function. For instance, the restricted Almon seems well suited for bell-shaped and slow-decaying weights, but somewhat less for fast-decaying weights, while unrestricted and Bernstein polynomials can perform fairly well irrespective of the underlying weighting structure.\footnote{However, we show that the U-MIDAS provides systematically less accurate parameter estimates compared to restricted Almon and Bernstein polynomials. These results are illustrated by the MSE in Figure A.2 in the Online Appendix.} The performance of Legendre and Chebychev polynomials improves substantially under flat weights. This weighting structure is nevertheless less likely to describe the actual temporal aggregation for most economic data.

Finally, we compare again the predictive accuracy of our BSGS-SS prior with respect to the BSGL and BAGL-SS benchmarks. Unlike Experiment 1, here we do not control directly for the degree of sparsity (in the coefficients of the approximating functions) within the groups, but this arises indirectly from the specification of DGPs 1 to 4 considered. The results reported in Figure (ref) suggest that the BSGS-SS prior provides a better predictive performance also in the mixed-frequency framework: the relative CRPS (with respect to the benchmarks) is on average less than one across the considered DGPs and polynomials (the best performing ones according to the results in Table (ref)), pointing to an average gain hovering around 5% to 10% overall.

Empirical application: nowcasting US GDP in a mixed-frequency framework

We use the proposed prior for a nowcasting exercise of US GDP in the mixed-frequency framework (ref) of Experiment 2, where $y_{t}=400\log(Y_{t}/Y_{t-1})$ is the annualized quarterly growth rate of GDP, and $\mathbf{x}_{t}$ is a vector of $122$ macroeconomic series sampled at monthly frequency and extracted from the FRED-MD database McCracken2016. The data sample starts in 1980Q1, while the pseudo out-of-sample analysis spans 2013Q1 to 2022Q4. Estimates are carried-out recursively using a rolling window of $T=132$ quarterly observations, and $h$-step-ahead posterior predictive densities are generated from (ref) through a direct forecast approach. We consider 3 nowcasting horizons ($h=0,1/3,2/3$) and the two lag polynomials that provide overall better results in the simulation analysis reported in Section (ref) -- namely the restricted Almon and the orthonormal Bernstein polynomials. Further, as the considered data sample spans periods characterised by strong and uneven macroeconomic fluctuations (permanent and transitory), such as the Great Moderation, the Great Recession, and the Covid crisis, we also allow for time-varying volatility with heavy tails and occasional outliers in the regression errors by using the model in Section (ref). Finally, we do not take into account real-time issues (ragged/jagged-edge data and revisions) and the data used for the analysis reflect the latest vintages available at the time of writing.\footnote{The GDP series used in the analysis was downloaded from FRED-MD database on July 2023 and we consider the 2023-06 vintage. We exclude from the analysis 5 series presenting sampling issues over the time span considered for the analysis (namely, new orders for consumer goods, non-borrowed reserves of depository institutions, 3-month AA financial commercial paper rate, 3-month commercial paper minus fed funds, consumer sentiment index).}\\ We consider several modelling strategies to exploit our bi-level sparsity prior approach. First, we estimate the forecasting model (ref) on the whole set of 122 indicators. Then, since the simulation results in Section (ref) point to a substantial deterioration of the selection and prediction performance in extreme sparse settings where $Ng\gg T$ then, we also estimate a separate model for each of the 8 groups of indicators partitioned as in McCracken2016.\footnote{The groups encompass i) output and income, ii) labour market, iii) housing, iv) consumption, orders, and inventories, v) money and credit, vi) interest and exchange rates, \textit{vii)} prices, and \textit{viii)} stock market. The number of indicators included in each group ranges between 5 and 31.} Summing up, we estimate a large set of alternative specifications, according to the 2 lag polynomials (Almon and Bernstein), the 5 volatility process (homoskedastic, SV, SV with Student-$t$ shocks, SV with outliers, SV with Student-$t$ shocks and outliers), and the 2 partition strategies (whole dataset \textit{vs} 8 groups). For each of these specifications we compute density forecasts and we combine them along several dimensions (across the 8 groups, the volatility process, and the lag polynomials, see Online Appendix A.4 for more details). Density forecasts combination is carried out using the optimal prediction pool proposed by Geweke2011, which relies on log-scores.

The combined density forecasts are then evaluated by the means of the average CRPS. The relative scores with respect to the AR(1) benchmark are reported in Table (ref). Three main findings can be outlined. First, Bernstein polynomials provide accurate density forecasts for $h=0$ and $h=1/3$, but the restricted Almon performs best at $h=2/3$. This outcome can be due to the way the different polynomials best describe the underlying weighting structure at each horizon. However, pooling these lag polynomials does not seem to provide here a systematic optimal strategy. Second, there is no clear-cut evidence in favour of a particular partition strategy. For $h=0$ and $h=2/3$, models including the whole set of indicators seem to perform better, or very closely, to those based on groups. Conversely, pooled forecasts from models based on group partitions tend to outperform for $h=1/3$. Overall, and consistently with the simulation results reported in Section (ref), this outcome suggests that the proposed procedure is robust to the presence of a large number of predictors. Finally, when compared to the BSGL prior (right hand side of Table (ref)), the results point to a considerable outperformance of our prior and substantial predictive gains for some horizons and pooling strategies.

Concluding remarks

This paper constructs optimal density forecasts for macroeconomic indicators in high-dimensional settings where the covariates present a known group structure. Such a structure arises, for instance, because variables in datasets are organised in groups or because of the flexible modelling of the relationship between the target variable and the set of covariates. The group-structure is important in order to reduce the dimensionality and we propose an optimal way to exploit it via the specification of a convenient prior distribution. By working under the assumption of bi-level sparsity, i.e. among and within groups, our procedure is able to detect the driving factors if only some regressors are relevant for forecast.

Density forecasts are constructed from the posterior predictive density, which we prove to present optimal asymptotic properties and to outperform many competitors in finite sample. We also provide parameter recovery and point forecasts, which we prove to be minimax-optimal for our procedure under some mild conditions. Unlike alternative approaches, our procedure does not require orthogonality among covariates belonging to different groups. Finally, we show that a more general error structure, such as stochastic volatility and ARMA, can be accommodated by only slightly modifying the proposed prior.

\paragraph{Acknowledgements.} The second author gratefully acknowledges financial support from Hi!Paris, and ANR-21-CE26-0003. The usual disclaimer applies. The views expressed in this paper are those of the authors and do not necessarily reflect those of the Banque de France or the Eurosystem.

\paragraph{Online Appendix.} It contains: additional elements about the prior, the Monte Carlo experiments and the application, guidelines for setting hyperparameters, the MCMC algorithm for the extended model with ARMA errors and stochastic volatility, additional theoretical results and all the proofs.

\onehalfspacing