Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
121,771 characters · 23 sections · 114 citation commands
Bayesian MIDAS Penalized Regressions: Estimation, Selection, and Prediction
\setcounter{page}{1} \setcounter{footnote}{0}
Mixed-data sampling (MIDAS) regressions provide a parsimonious and theoretically efficient treatment of the time-aggregation problem when data are sampled at different frequencies Ghysels2005,Ghysels2007,Andreou2010. They have been therefore intensively and successfully used to forecast low frequency series, such as GDP, using monthly, weekly or daily predictors Clements2008,Clements2009,Kuzin2011,Andreou2013. 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; Foroni2015) with automatic model selection of relevant predictors and high-frequency lags Castle2010,Bec2015, factor-MIDAS regressions Marcellino2010, and targeted factor-MIDAS regressions (Bai2008; see also Bessec2013, Bulligan2015, and Girardi2017).
Recently, the literature has been increasingly focusing on Machine Learning and penalized regression techniques for macroeconomic applications in a high-dimensional environment (see Korobilis2013; NG2013752; Gefang2014; ng_2017; Koop2019; Korobilis2019; to name only a few). Nevertheless, only a few contributions have paid attention to MIDAS regressions in high-dimension so far Marsilli2014,Siliverstovs2017,Uematsu2019,Babii2020.\footnote{Marsilli2014 proposes a functional MIDAS combined with a Lasso objective function, which is solved in 1-step through a non-linear optimization algorithm. Siliverstovs2017 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, Babii2020 have proposed a sparse-group Lasso estimator for MIDAS regression based on Legendre polynomials, which is closely related to our approach.} In particular, Uematsu2019 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 Zhang2014,Zhao2015,Xu2015,Rockova2018. Finally, penalty hyper-parameters governing the group selection are here tuned through a computationally efficient and data-driven approach based on stochastic approximations Atchade2011,Atchade2011a 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 lounici2011. 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 Ning2020 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 Rockova2018 and Baiinpress. 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 (ref) and (ref) introduce the MIDAS penalized regressions and the Bayesian MIDAS framework. Section (ref) presents the theoretical analysis. In Section (ref) we describe the Gibbs sampling and we discuss the Empirical Bayes approach used to automatically tune the penalty hyper-parameters. Section (ref) presents simulation experiments and Section (ref) provides an empirical application to the U.S. GDP growth. Finally, Section (ref) concludes. The appendix contains technical proofs on the theorems reported in Section (ref), while additional results and proofs of technical lemmas are collected in the Supplementary Appendix available on-line.
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, Andreou2010 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:
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” Cooper1972, Equation (ref) with Almon lag polynomials can be reparameterized as:
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\vert T-h}^{}$, conditional on sample information available up to $T-h$, can be hence obtained using (ref):
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 (ref) 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)\vert_{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 Smith1976. As a result, the number of parameters in (ref) reduces from $\sum_{k=1}^{K}(p_{k}+1)$ to $\sum_{k=1}^{K}(p_{k}-r_{k}+1)$, where $r_{k}\le 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.
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.\footnote{The direct method used in regression (ref) may be also hampered by multicollinearity in the artificial variables $\mathbf{Z}_{t}^{(m)}$ Cooper1972. 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 Tibshirani1996 have been recently investigated in a high-dimensional mixed-frequency framework by Uematsu2019 and Lima2020. Cognizant of the limits of the Lasso (estimation bias and lack of selection consistency), the Adaptive Lasso (Zou2006; 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.\footnote{The 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 (ref) 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 Efron2004 and Zou2005, 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 (Yuan2006, Wang2008; 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 $0$ and regressors $\mathbf{Z}_{t}^{(m)}$ in (ref) 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:
where $\mathcal{L}^{}_{T}(\boldsymbol\theta)$ denotes the negative log-likelihood function, $\Vert\boldsymbol\theta_{j}\Vert_{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 Wang2008. However, as suggested by Callot2014, 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.
Since the work of Yuan2006, 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 Wang2008. In this paper, we consider instead a Bayesian hierarchical approach based on MCMC methods Kyung2010, 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 Khare2013.\footnote{It is nevertheless worth noting that the results in Khare2013 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 Roy2017).} 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.
As noted by Tibshirani1996, the Lasso estimator can be interpreted as the Bayes posterior mode using normal likelihood and independent Laplace (double-exponential) prior for the regression coefficients. Park2008 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 Kyung2010 and Leng2014, 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 (ref)):
This suggests the following hierarchical representation of the Bayesian MIDAS Adaptive Group Lasso model (BMIDAS-AGL):
where $y:=(y_{1},\dots,y_{T})^{\prime}$ is centered at $0$ 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} \vert \mathbf{0},\sigma/\lambda_{j}) \propto (\lambda_{j}/\sigma)^{g_{j}}\exp(-\lambda_{j}\Vert\boldsymbol\theta_{j}\Vert_{2}/\sigma)$ (see Appendix (ref)). 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:
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 Zhang2014,Zhao2015,Rockova2018. In the present study, we follow Xu2015 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:
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 Ning2020, 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 (ref) (see also Section (ref) 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 (ref)). The full posterior distribution of the unknown parameters conditional on the data and for some penalty hyper-parameters $\boldsymbol\lambda$ is:
This section provides the theoretical validation of our Bayesian MIDAS procedure. In Section (ref) we investigate the asymptotic properties of our procedure by adopting a frequentist point of view. In Section (ref) 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.
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 (ref). 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 := \vert S_0\vert$ 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$.
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 (ref). By Assumption (ref) (iii), $\boldsymbol{\theta}_0\in\overline\Theta_0$.
Assumption (ref) (ii) restricts the rate at which the maximum group size increases to infinity. Assumption (ref) (iii) restricts the growth rate of the maximum component of the true $\boldsymbol{\theta}_0$. Assumption (ref) (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 (ref) 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 (ref) 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 Castillo2015 but it excludes a $v$ that increases with $G$.
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}',\sigma^2)'$ is endowed with the BMIDAS-AGL-SS prior (ref). The first theorem establishes posterior consistency of the in-sample prediction error.
The assumption $\epsilon \rightarrow 0$ in Theorem (ref) 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. lounici2011). The contraction rate $\epsilon$ is the same as the one in Ning2020 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 lounici2011. 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 (ref) (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 lounici2011). 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. buhlmann2011, and Castillo2015).\\ The proof of Theorem (ref) is provided in Appendix (ref) and follows the strategy of Ning2020 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\'{e}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\'{e}nyi divergence, we use the general posterior contraction procedure for independent observations, as in ghosal2007. That is, we first show that our prior (ref) 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 (ref), 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 bickel2009, 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. Castillo2015, and Ning2020) that if $\boldsymbol{\theta}$ is sparse, then local invertibility of the Gram matrix $\mathbf{Z}'\mathbf{Z}$ is sufficient for estimability of $\boldsymbol{\theta}$. We then introduce the following quantity.
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}'\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 Ning2020.
The proof of Theorem (ref) is provided in Appendix (ref). 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 lounici2011. The posterior contraction rate provided in Theorem (ref) deteriorates when $\widetilde\phi(s_0 + M_2\widetilde s_0)$ is small.
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 $$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}|\boldsymbol{\theta},\sigma^2,\mathbf{Z}_{\tau-h})\Pi(\boldsymbol{\theta},\sigma^2|y,\mathbf{Z})d\boldsymbol{\theta} d\sigma^2,\qquad \forall \tau > T,$$ 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}'\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}'\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}',\sigma^2)'$ is endowed with the BMIDAS-AGL-SS prior (ref). We start by establishing the contraction rate of the posterior of the out-of-sample prediction error.
The proof of Theorem (ref) is provided in Appendix (ref). Result (ref) 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 CarrascoRossi2016, which in turn implies $\|\mathbf{Z}_{\tau - h}\|_2 = O_p(1)$, then we would get the same rate as in Theorem (ref). 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}'\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.
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}']\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}'\boldsymbol{\theta}|^2}{\sum_{j=1}^G \|\mathbf{Z}_j' \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 (ref).\\
While Theorem (ref) 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:
where $K(f_1\Vert 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 (ref) is equal to $\int f_{\mathbf{Z}}(y) dy \,K(f_{\mathbf{Z}_{\tau-h},y,\mathbf{Z}} \Vert \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}',\sigma^2)'$ 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. Aitchison1975).
We now provide details on the Gibbs sampler we use to simulate from the posterior distributions (ref) and (ref). For this purpose, we consider an efficient block Gibbs sampler Hobert1998. Let us denote $\boldsymbol\theta_{\mathbin{\mathchoice{\hbox{\tikz{\draw[line width=0.6pt,line cap=round] (3pt,0) -- (0,6pt);}}}{\hbox{\tikz{\draw[line width=0.6pt,line cap=round] (3pt,0) -- (0,6pt);}}}{\hbox{\tikz{\draw[line width=0.45pt,line cap=round] (2pt,0) -- (0,4pt);}}}{\hbox{\tikz{\draw[line width=0.4pt,line cap=round] (1.5pt,0) -- (0,3pt);}}}} 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{\tikz{\draw[line width=0.6pt,line cap=round] (3pt,0) -- (0,6pt);}}}{\hbox{\tikz{\draw[line width=0.6pt,line cap=round] (3pt,0) -- (0,6pt);}}}{\hbox{\tikz{\draw[line width=0.45pt,line cap=round] (2pt,0) -- (0,4pt);}}}{\hbox{\tikz{\draw[line width=0.4pt,line cap=round] (1.5pt,0) -- (0,3pt);}}}} j}$ partitions of the design matrix corresponding to $\boldsymbol\theta_{j}$ and $\boldsymbol\theta_{\mathbin{\mathchoice{\hbox{\tikz{\draw[line width=0.6pt,line cap=round] (3pt,0) -- (0,6pt);}}}{\hbox{\tikz{\draw[line width=0.6pt,line cap=round] (3pt,0) -- (0,6pt);}}}{\hbox{\tikz{\draw[line width=0.45pt,line cap=round] (2pt,0) -- (0,4pt);}}}{\hbox{\tikz{\draw[line width=0.4pt,line cap=round] (1.5pt,0) -- (0,3pt);}}}} 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:
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{\tikz{\draw[line width=0.6pt,line cap=round] (3pt,0) -- (0,6pt);}}}{\hbox{\tikz{\draw[line width=0.6pt,line cap=round] (3pt,0) -- (0,6pt);}}}{\hbox{\tikz{\draw[line width=0.45pt,line cap=round] (2pt,0) -- (0,4pt);}}}{\hbox{\tikz{\draw[line width=0.4pt,line cap=round] (1.5pt,0) -- (0,3pt);}}}} j}\boldsymbol\theta_{\mathbin{\mathchoice{\hbox{\tikz{\draw[line width=0.6pt,line cap=round] (3pt,0) -- (0,6pt);}}}{\hbox{\tikz{\draw[line width=0.6pt,line cap=round] (3pt,0) -- (0,6pt);}}}{\hbox{\tikz{\draw[line width=0.45pt,line cap=round] (2pt,0) -- (0,4pt);}}}{\hbox{\tikz{\draw[line width=0.4pt,line cap=round] (1.5pt,0) -- (0,3pt);}}}} 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 :
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{\tikz{\draw[line width=0.6pt,line cap=round] (3pt,0) -- (0,6pt);}}}{\hbox{\tikz{\draw[line width=0.6pt,line cap=round] (3pt,0) -- (0,6pt);}}}{\hbox{\tikz{\draw[line width=0.45pt,line cap=round] (2pt,0) -- (0,4pt);}}}{\hbox{\tikz{\draw[line width=0.4pt,line cap=round] (1.5pt,0) -- (0,3pt);}}}} j}\boldsymbol\theta_{\mathbin{\mathchoice{\hbox{\tikz{\draw[line width=0.6pt,line cap=round] (3pt,0) -- (0,6pt);}}}{\hbox{\tikz{\draw[line width=0.6pt,line cap=round] (3pt,0) -- (0,6pt);}}}{\hbox{\tikz{\draw[line width=0.45pt,line cap=round] (2pt,0) -- (0,4pt);}}}{\hbox{\tikz{\draw[line width=0.4pt,line cap=round] (1.5pt,0) -- (0,3pt);}}}} 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 (ref), Assumption (ref)), and
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\vert \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. Park2008 and Kyung2010 suggest to address this issue by implementing the Monte Carlo EM algorithm (MCEM) proposed by Casella2001, which complements the Gibbs sampler and provides marginal maximum likelihood estimates of the hyper-parameters when the marginal distribution $f(y\vert\boldsymbol\lambda,\mathbf{Z})=\int f(y\vert\boldsymbol\phi,\boldsymbol\lambda,\mathbf{Z})\Pi(\boldsymbol\phi\vert\boldsymbol\lambda)d\boldsymbol\phi$ is not available in closed form. This can be easily obtained by repeatedly maximizing $N$ times the expectation function
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\vert \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\vert\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 Yuan2006. To deal with this issue, we adopt an alternative Empirical Bayes approach that relies on the so-called internal adaptive MCMC algorithms (see Atchade2011a). 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 Atchade2011, 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:
where $H(\boldsymbol\lambda,\boldsymbol\phi):=\nabla_{\boldsymbol\lambda}\log\left[f(y\vert\boldsymbol\phi,\boldsymbol\lambda,\mathbf{Z})\Pi(\boldsymbol\phi\vert\boldsymbol\lambda)\right]=\nabla_{\boldsymbol\lambda}\log\Pi(\boldsymbol\phi\vert\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)$ Lange1995. If the integral $\int H(\boldsymbol\lambda^{(s)},\boldsymbol\phi)\Pi(\boldsymbol\phi\vert \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\vert \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 (ref) the function $H(\boldsymbol\omega,\boldsymbol\phi)=\nabla_{\boldsymbol\omega}\log\Pi(\boldsymbol\phi\vert\boldsymbol\omega)$ takes the form:
where $\mathbf{g}=(g_{1}^{},\dots,g_{G}^{})^{\prime}$ and $\odot$ is the element-wise product. Hence, the updating rule for $\boldsymbol\omega$ is:
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; Andrieu2005; Atchade2011) ensuring the convergence of $\boldsymbol\lambda$ and the posterior distribution of $\boldsymbol\phi$ towards $\widehat{\boldsymbol\lambda}$ and $\Pi(\boldsymbol\phi\vert \widehat{\boldsymbol\lambda},y,\mathbf{Z})$, respectively. Details on the stabilization algorithm are reported in Appendix (ref).
We illustrate the main features and the computational advantage of the proposed methodology using simulated data generated from the MIDAS model (ref) 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 (ref) 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 (ref). 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 (ref). 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.
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:
where $\widetilde{B}\left(c;\boldsymbol\theta\right)$ denotes the normalized weigths (i.e. summing up to 1). Following Andreou2010, 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 (ref). 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:
where $\boldsymbol\Sigma_{\varepsilon}$ has elements $\sigma_{\varepsilon}^{\vert k-k^{\prime} \vert}$, 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 (ref)). The hyper-parameters $\boldsymbol\lambda$ are tuned using the stochastic approximation approach described in Section (ref), 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.
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 (BIAS$^2$) over $R$ Monte Carlo replications and the full set of $K$ estimated overall slope parameters $\widehat{\boldsymbol\beta}$.\footnote{For $R$ Monte Carlo replications, $K$ variables, and $S$ Gibbs draws, we have that:
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 Kyung2010, 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 Barbieri2004, that is, under some conditions, a soft thresholding estimator presenting model selection consistency and optimal asymptotic estimation rate Xu2015.
Forecasts are obtained from the following posterior predictive density for $y_{T\vert T-h}^{}$:
where $\Pi(\boldsymbol\phi,\boldsymbol\lambda\vert\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 (ref) and (ref), draws $y_{T\vert T-h}^{(s)}$, $s=1,\dots,S$, from the predictive distribution (ref) can be obtained directly from the Gibbs sampler.\footnote{It 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\vert T-h}^{}=(S-\bar{s}+1)^{-1}\sum_{s=\bar{s}+1}^{S}y_{T\vert 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 Mitchell2011, 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 Gneiting2007.
Simulation results are reported in Table (ref) 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 (ref), 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.\footnote{For $w_{\mathcal{A}}=K_{\mathcal{A}}/K$ and $w_{\mathcal{A}^{c}}=1-w_{\mathcal{A}}$, we have that:
} Looking at the MSE by active/inactive set reported in Table (ref) (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.\footnote{For 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 (ref) 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 (ref) and (ref) and the latter using the algorithm described in Pettenuzzo2016 on the set of relevant variables only.\footnote{We consider the same restricted Almon lag polynomial as for our models. Further, for the Oracle BMIDAS, we follow Pettenuzzo2016 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 (ref) seems to confirm the intuition discussed in Section (ref), 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 (ref) 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 (ref) 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 (ref) the forecasting performance of our models, in relative terms with respect to the Oracle for ease of comparison.\footnote{We 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 (ref)) 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 Uematsu2019: the Lasso (L), the Elastic-Net (EN; Zou2005), and two folded-concave penalizations, such as the smoothly clipped absolute deviation (SCAD; Fan2001) and the minimax concave penalty (MC$+$; Zhang2010). 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 Fan2011. 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 Uematsu2019, 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.
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 McCracken2016.\footnote{We 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 Uematsu2019), 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 Andreou2013,Pettenuzzo2016,Adrian2019: 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 Fama1993; returns on the portfolio of high minus low book-to-market ratio stocks studied by Fama1993; 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 Aruoba2009 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 (ref) 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 Diebold1995 and West1996 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 Harvey1997. As a robustness check, we further consider forecasts from the following competing models:
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 (ref) and (ref). In the first row of Table (ref), 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 (ref)). 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 Uematsu2019, 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 Giannone2017.
Figure (ref) 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 McCracken2016. 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 Andreou2013 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.
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 Carriero2015, Pettenuzzo2016, as well as quantile mixed-frequency regressions Lima2020. Further, recent research has been focusing on factor-adjusted sparse regressions to deal with high-dimensional regressions and highly correlated predictors Kneip2011,Fan2020. We believe that these extensions to our models represent interesting paths for future research.
\setcounter{equation}{0} \setcounter{table}{0} \setcounter{figure}{0} \setcounter{subsection}{0}