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.
74,654 characters · 12 sections · 70 citation commands
Horseshoe Prior Bayesian Quantile Regression
\def\spacingset#1{ {#1}} \spacingset{1}
\def\spacingset#1{ {#1}} \spacingset{1}
\if00 \fi
\if10 {
} \fi
{\it Keywords:} Global-Local Priors, Monte Carlo, Shrinkage, Machine Learning, Quantile Regression. \\ {\it JEL:} C110, C530, C550, C63
\spacingset{1.45}
Quantile regression has been an important tool in the econometricians’ toolkit when estimating heterogeneous effects across the conditional response distribution since the seminal work of koenker1978regression. In contrast to least squares methods, it estimates quantiles of the dependent variables’ conditional distribution directly, which allows for richer inference than solely focusing on the conditional mean. While highly influential in the risk-management and finance literature in calculating risk measures such as VaR (i.e. the loss a portfolios value incurs at a specified probability level), it has experienced a recent surge in use especially in the macroeconomic literature to quantify downside risks of the aggreagate economy to financial conditions. This is of special interest to policy institutions due to a stronger macroprudential focus in the wake of the financial crisis. Quantile regression has recently been used in particular to quantify risks and vulnerabilities of output growth to indicators about financial health adrian2019vulnerable,figueres2020vulnerable,adams2020forecasting,hasenzagl2020financial, appraise likelihoods of scenario analyses prasad2019growth, monitor downside risks to output growth in real time hasenzagl2020financial,carriero2020nowcasting and forecast density construction korobilis2017quantile,mazzi2019nowcasting.
A challenge for these purposes is that sources of risk can be numerous such that simple quantile regression is often rendered imprecise or infeasible in high dimensions. While a variety of regularization and dimension reduction techniques have been proposed for macroeconomic forecasting, stock2002forecasting,stock2012generalized,kim2014forecasting,bai2008forecasting,de2008forecasting, extensions of high dimensional methods, in particular Bayesian methods, applied to quantile regression, remain relatively scant.
The Bayesian quantile regression approach, as popularized by yu2001bayesian, is based on the asymmetric Laplace likelihood (ALL), which has a special connection to the frequentist quantile regression solution, in that its maximum likelihood estimates are equivalent to traditional quantile regression with a check-loss function koenker2005. A hurdle in the Bayesian literature has been that ALL based methods result in improper posteriors with any but non-informative or exponential Laplace priors, where the latter results in the popular Bayesian Lasso quantile regression li2010bayesian,alhamzawi2013conjugate,alhamzawi2012bayesian,chen2013bayesian. The broader Bayesian shrinkage literature has shown, however, that global-local shrinkage priors such as the horseshoe carvalho2010horseshoe and Dirichlet-Laplace prior bhattacharya2016fast offer asymptotic as well as computational advantages over the former methods bhadra2019lasso. These methods have not yet been considered for the Bayesian quantile regression. The aim of this paper is to bridge this gap and extend the global local prior to quantile regression.
This paper's primary contribution is twofold. First, we derive the horseshoe prior of carvalho2010horseshoe for the Bayesian quantile regression framework (BQR) of yu2001bayesian. Second, we develop an efficient posterior sampler for the quantile specific regression coefficients based on data augmentation akin to bhattacharya2016fast which speeds up computation significantly for high dimensional quantile problems.
To showcase the performance of the horseshoe BQR (HS-BQR) we provide a large scale monte carlo study as well as a high dimensional VaR application to U.S. GDP (often called GaR in the literature). In the Monte Carlo study we show that the proposed estimator provides more stable and at worst, similar performance compared to a variety of Bayesian lasso quantile regression methods in terms of coefficient bias and forecast accuracy. We find that, particularly, tails of the distributions are consistently better estimated by the HS-BQR which echos findings from the Bayesian VaR literature chen2012forecasting. In the GaR application we show that the HS-BQR produces better calibrated forecast densities compared to the Bayesian alternatives and importantly provides the best performance for lower and upper tails which makes it a powerful tool for recession probability monitoring. The framework provided in this paper has the additional advantage that the derived algorithms can be directly applied to other global-local priors\footnote{For an overview of global-local priors see polson2010shrink.} that can be expressed as scale mixture of normals.
In what follows, we will first review the generic quantile regression framework. Then, we will present the main results in the Bayesian quantile as well as shrinkage literature which have motivated the form of the model. Following this, we will develop posteriors as well as the sampling algorithm. Lastly, we will provide evidence from Monte Carlo simulations and an empirical application of the favorable performance of the HS-BQR compared to alternative methods. We conclude with further generalizations of the algorithms provided and a discussion of our results.
Taking the linear model $Y=X\beta + \epsilon$ as our starting point, the conditional quantile function of $Y$ can be defined as
where $p \in (0,1)$, X is a $T \times K$ matrix of covariates, $\beta(p)$ is a $K \times 1$ vector of quantile specific regression coefficients, and $\epsilon$ is a $T \times 1$ vector of residuals which follow some unspecified distribution. Unlike in classical regression analysis, quantile regression does not make any parametric assumption about $\epsilon$ koenker2005 which allows for rich, non-symmetric inference about the conditional distribution of Y.
While it is possible to estimate an infinite amount of quantiles, in practice one only estimates a finite number of quantiles which are of interest, as the number of distinctly estimable quantiles increase roughly linearly with sample size buchinsky1998recent,davino2013quantile.
The objective function of the quantile regression can be represented by the following minimization problem:
where $\rho_{p}(.)$ is a loss function with the following form:
where I(.) is an indicator function taking on a value of 0 or 1 depending on whether the condition is satisfied. Equation ((ref)) determines the weight each observation receives in the minimization problem. It is often referred to as the check-loss function due to the weight profile it assigns depending on the quantile being estimated koenker2005. Note how $(y_t-x_t'\beta)$ is the residual of a regression model. The interpretation of the coefficients is thus similar to the classical regression case: $\beta_j(p)$ is the rate of change of the $p^{th}$ quantile of the dependent variable's distribution to a unit change in the $j^{th}$ regressor.
We assume the quantile regression model ((ref)) and a fixed design X. As shown by yu2001bayesian, $\beta(p)$ can be obtained as the maximum likelihood estimator for $\beta$ under the fully parametric model $y_t=x_t'\beta + \epsilon$ where $\{ \epsilon_t\}_{t=1}^T$ are assumed i.i.d. with common density given by
where $\mathbb{R}_{+}:=(0,\infty)$ and $\mathbb{R}_{-}:=(-\infty,0]$. The errors follow an asymmetric Laplace density with the $p^{th}$ quantile equal to zero. Assuming the linear model as above with error density (4), the joint likelihood $f(Y|\beta,\sigma)$ becomes:
It is apparent that using any non-trivial prior for $(\beta,\sigma)$, will result in an intractable posterior which will lead to inefficient accept and reject sampling algorithms yu2001bayesian. However, kozumi2011gibbs showed using the mixture representation of the asymmetric Laplace distribution provided by kotz2012laplace, that the likelihood in ((ref)) can be obtained by formulating the error process as:
where $z \sim exp(1)$, $u \sim N(0,1)$, while $\theta=\frac{1-2p}{p(1-p)}$ and $\tau^2=\frac{2}{p(1-p)}$ are deterministic quantile specific parameters. Let $\theta=\theta(p)$ and $\tau^2=\tau^2(p)$ be defined as above and let $\{(y_t,z_t)\}_{t=1}^T$ be independent random pairs. Now, to simplify the Gibbs sampler, we instead assume $z_t \sim exp(\sigma)$ such that given $z_t$, $y_t$ is in normally distributed as $y_t|z_t \sim N(x_t^{'}\beta+\theta z_t,z_t\sigma\tau^2)$. The joint density of $Y|Z$ is given by:
Although the scale of the likelihood is non-standard, kozumi2011gibbs show that independent normal-inverse-gamma (N-IG) priors result in conditionally conjugate posteriors which we exploit for the horseshoe prior adaptation.
In order to identify the posterior of a large dimensional coefficient vector in small samples, informative priors are needed. Ideally, these priors are able to separate noise variables from signals such that the noise is shrunk towards zero and signals attain their unrestricted parameter values. In the frequentist setting, this is usually achieved through penalized regression which forces variables to threshold to 0 or not. In the Bayesian approach, it is important to note that the assumption about sparsity is fundamentally different in that proper prior distributions have non-zero probability over sparse and non-sparse regions in the posterior batencourtsparse. In order, therefore, to achieve the desired separation between shrunk and unshrunk variables, the amount of shrinkage on a 0-1 scale should approach a bi-modal distribution where most of the mass is on 0 and 1 respectively. The horseshoe prior of carvalho2010horseshoe achieves such a shrinkage profile, while double-exponential based lasso priors do not.
The idea of the global-local family of shrinkage priors as defined by polson2010shrink is to apply a scale mixture of normal prior to the regression coefficients with global scale-prior $\nu^2$, controlling the overall shrinkage applied to the regression and a local scale $\lambda^2_j$ which allows for the local possibility of regressors to escape shrinkage when they have large effects on the response. Expressed differently, global-local shrinkage priors define a distribution for a shrinkage factor, $\kappa_j, \; \text{for} \; j=(1,\cdots,K)$, which is bounded to be between 0 and 1, and whose distribution is implicitly defined by the shape of the scales, $\nu$ and $\lambda$. Under certain conditions (see piironen2017sparsity), this shrinkage factor for generalized linear models for any global-local prior can be shown to have the following form:
where $s_j$ refers to the column wise variance of X. One can use this shrinkage factor to approximate the mean of the posterior coefficient vector as:
where $\hat{\beta}$ refers to the maximum likelihood estimate. The distribution of $\kappa_j$ implied by the horseshoe prior and lasso prior are plotted in figure ((ref)). The focal point of this paper, the horseshoe prior of carvalho2010horseshoe, employs two half Cauchy distributions for $\lambda$ and $\tau$:
which by the change of variables theorem imply a Beta(0.5,0.5) distribution on the shrinkage factors carvalho2010horseshoe. When $\tau$ and $\lambda$ are strongly identified, this prior results in complete or no shrinkage for each coefficient in the limit, as can be visually confirmed from figure ((ref)). Intuitively, this shape is induced by the Cauchy distribution having most mass on zero with fat enough tails to allow signals to escape shrinkage.
The lasso prior of park2008bayesian cast into the global local form, instead is defined as\footnote{Note here that that an auxiliary variable u is introduced which after integration yields the desired double exponential Laplace density on the coefficient vector (see park2008bayesian)}:
which, as can be seen from figure ((ref)), has the unfortunate trait that large signals can escape shrinkage, but noise variables are not shrunk aggressively enough. This will result in too little shrinkage in large dimensional problems with many noise variables.
In order to generalize the horseshoe prior to the ALL in ((ref)), it needs to be formulated under the assumption of independence between the $\beta$ and $\sigma$ prior so that the posterior takes a conditionally normal form kozumi2011gibbs. While prior independence between regression coefficients and error variance might seem like a strong assumption, moran2018variance have shown that in high-dimensional settings, the independence assumption aids inference of the error variance. This is due to the fact that conjugate priors act mathematically as additional observations which artificially bias the error variances downwards when K$>>$T.
The general independent global-local prior takes the following hierarchical form:
As recommended by gelman2006prior, we select a weakly informative inverse Gamma distribution as the prior for the error variance $\sigma^2$ instead of non-informative priors and two independently distributed half Cauchy distributions on the positive support for the scale parameters of the $\beta$ prior:
Under the above priors, the posteriors take the following form:
where $\overline{V}= (X'UX+\underline{V})^{-1}, U=diag(\frac{1}{\tau^2z_i\sigma})$ and $\underline{V}=\nu^2diag(\lambda^2_1,\cdots,\lambda^2_K)$, and $\overline{\beta}$ is defined as $\overline{\beta}=(X'UX+\overline{V}^{-1})^{-1}(X'Uy+\overline{V}^{-1}\underline{\beta})$ with $\underline{\beta}=\boldsymbol{0}_K$. The conditional posterior of the scale parameter is given by
and $\overline{a}=\underline{a}+\frac{3T}{2}$, $\overline{b}=\underline{b}+\sum_{t=1}^T\frac{(y_t-x_t'\beta-\theta z_t))^2}{2z_t+\tau^2}+\sum_{t=1}^Tz_t$.
Due to the assumption of independence on the scales $(\lambda,\nu)$, it is straightforward to show that their posteriors follow independent Cauchy distributions. Since the Cauchy has no defined moments which would enable sampling, the literature has proposed Gibbs samplers which rely either on slice sampling polson2014bayesian or mixture representations makalic2015simple. Unfortunately, both rely on conjugate formulations. We use the fact that the posterior distribution of $\lambda_j$ conditional on $\nu$ remains independent of all other parameters by assumption, to formulate a block slice sampling algorithm for for $\lambda=(\lambda_1,\cdots,\lambda_K)'$ akin to polson2014bayesian which is given in (ref). Finally, since the latent $z_t$ are sampled independently, the conditional posterior follows the reciprocal of the inverse Gaussian:
Where I-G stands for the inverse Gaussian distribution with location and rate parameters respectively, $\overline{c_t}=\frac{\sqrt{\theta^2+2\tau^2}}{|y_t-x_t'\beta|}$ and $\overline{d_t}=\frac{\theta^2+2\tau^2}{\sigma\tau^2}$. \\
With these conditional posteriors at hand, we utilize a standard Gibbs sampler. The dynamics of the Markov chain $\{(\beta_m,\sigma_m,\lambda^2_m,\nu^2_m,z_{m})\}_{m=0}^{\infty}$ are implicitly defined through the following steps
khare2012geometric show that the Markov chain of this sampler is geometrically ergodic and also valid in K$>>$T settings which gives theoretical justification to apply this sampler to high dimensional settings. However, a computational bottleneck is present in very high dimensions in evaluating the $K\times K$ dimensional inverse for the conditional posterior of $\beta$. Cholesky decomposition based methods will generally be of order $O(K^3)$. Taking into consideration that in quantile settings, one is usually interested in obtaining more than one expected quantile, this can result in prohibitively long computation times. We therefore provide a more efficient sampling algorithm for $\beta$ which leverages data augmentation similar to the algorithm developed by bhattacharya2016fast which is of order $O(T^2K)$ and especially beneficial in high dimensional settings.
As derived above, using the scale mixture representation in ((ref)), the conditional posterior of $\beta$ given all other parameters can be written as:
Suppose, we want to sample from $N_K(\mu, \Sigma)$, where
Assume $D \in \mathbb{R}^{K \times K}$ is a positive definitive matrix, $\boldsymbol{\phi} \in \mathbb{R}^{T \times K}$, and $\alpha \in \mathbb{R}^{T \times 1}$. Then ((ref)) is a special case of of ((ref)) when setting $\Phi=\sqrt{U}X$, $D=\Lambda_{*}$ and $\alpha=\sqrt{U}y$. An exact algorithm to sample from ((ref)) is thus given by:
\newmdtheoremenv{theo}{Algorithm}
Proof. Using the Sherman-Morrison-Woodbury identity, $\mu=D\Phi'(\Phi D \Phi' + I_T)^{-1}\alpha$. Plugging in 2. into 3., we obtain $\theta = u + D\Phi'(\Phi D \Phi'+I_T)^{-1}(\alpha-\xi)$. Since by definition $\xi \sim N(0,\Phi D \Phi + I_K)$, $\theta$ follows a normal distribution with mean $D \Phi'(\Phi D \Phi'+ I_K)\alpha=\mu$. As $cov(u,\xi)=D \Phi'$, it follows that $cov(\theta)=D-D \Phi'(\Phi D \Phi' + I_K)^{-1}\Phi D$ which by the Sherman-Morrison-Woodbury identity is equal to $\Sigma$. More details are provided in ((ref)). The provided algorithm is not specific to the horseshoe prior and follows through for any prior of the form in ((ref)). The computational advantage provided in algorithm 1 compared to Cholesky based decompositions is that we can cheaply sample from $(u,\xi)'$ which via data augmentation yields samples from the desired distributions.
In order to verify the theoretical advantages of the HS-BQR over the exponential based quantile regression priors laid out above, we conduct a variety of high dimensional Monte Carlo simulations that test the priors' ability to adapt to different degree's of sparsity and error distribtions in the data generating processes. We consider 3 variants of the original lasso prior which have been adapted to the Bayesian Quantile regression:
Three sample sizes are considered: $T_i \in \{200,500,1000\}$ \footnote{$T_3$ is only considered for the HS-BQR due to the prohibitively long computation times of the lasso based Gibbs samplers.}. In total 100 Monte Carlo datasets were generated\footnote{except for $T=500$ for the block case where only 20 Monte Carlo experiments were done due to the time it takes to run the estimator on such large dimensions.} for which the last 100 observations are constructed to be the same for each $T_i$ in order to make forecast errors comparable. The remainder of the observations are used as training samples to retrieve the mean posterior $\hat{\beta}(p)$ vector to calculate bias\footnote{Alternatively, one could also use the MAP (mean-absolute-posteriori) estimate of the regression posterior as the point estimate. This might seem more natural when comparing Bayesian quantile regression methods to frequentist quantile estimators due to their equivalence as discussed in kozumi2011gibbs. We found that since the conditional posteriors are normal, there is no practical difference between the posterior mean and MAP.}.
We consider 12 designs in total which vary along two different dimensions: the degree of sparsity and the error generating process. We test the following sparsity patterns:
Notice that for $T_1$ there are always more coefficients than observations.
Consider a linear model as in ((ref)). To retrieve the true quantile regression coefficients, $\beta(p)$, we make use of koenker2005's alternative representation of the quantile regression:
where $u_t$ is assumed to be i.i.d. having some CDF, $F$. The dimensionality of $\vartheta$ is K$\times$1 and determines which covariates have non constant quantile functions. This can be seen from the solution for $\beta(p)$ to equation ((ref)):
Hence, the true $\beta(p)$ profile of a quantile regression model has a random coefficient model interpretation, where the vector of coefficients can be decomposed into a fixed plus a random component. In particular, the random component depends on the inverse CDF of the error, $F^{-1}(p)$. One can therefore think of $\vartheta$ as determining which variable is correlated with the error, where by default the first entry, $\vartheta_0$, is set to 1. This entails that location effects will always be present.\footnote{While it is possible for $\vartheta$ to take on any value, for simplicity we assume that the elements of $\vartheta$ only to take on the values \{0, 1\}.}
From a frequentist' perspective Equation ((ref)) is our oracle estimator for $\beta(p)$ for a given quantile $p$, which, given that the ALL approximation in equation ((ref)) holds, can be compared to the mean of the posterior of equation ((ref)) kozumi2011gibbs. With this in mind, it is trivial to calculate the true $\beta$'s for the error generating processes considered.
The second dimension along which the DGPs differ is in their error process. The proposed DGPs can be grouped into two broad cases: (1) i.i.d. errors ($y_1$ and $y_2$); and (2) heteroskedastic errors ($y_3$ and $y_4$). In $y_1$, we assume that the error distribution follows a standard normal distribution and in $y_2$, the error has student-t distributed errors with 3 degrees of freedom. For the other cases, we assume simple heteroskedasticity caused by correlation between the second covariate (whose coefficient we denote as $\beta_1$) and $\epsilon$. Lastly, $y_4$ can be thought of as containing a mixture between a uniform and a standard normal error distribution. In all simulations, the design matrix is simulated using a multivariate normal distribution with mean 0 and a covariance matrix with its $(i,j)^{th}$ element defined as $0.5^{|i-j|}$.
Relating the assumed error processes to the random coefficient representation ((ref)), it is clear that, under i.i.d. errors, only the constant has a non-constant quantile function caused by $F^{-1}$ (hereinafter called location shifters). Under the heteroskedastic designs, apart from the constant, $\beta_1$ will have a non-constant quantile function as well. Hence, $\beta_1$ in $y_3$ is determined by $F^{-1}_{N(0,1)}$ across p, and $\beta_1$ in $y_4$ follows $F^{-1}_{U(0,2)}$, i.e., increases linearly with p. The simulation designs (and the true quantile functions) are summarized in table ((ref)).
We evaluate the performance of the estimators in terms of bias in the coefficients and forecast error. Using the true quantile profile in $\beta(p)$ in ((ref)), we calculate root mean coefficient bias (RMCB) and root mean squared forecast error (RMSFE) as:
where $iter$ is the number of Monte Carlo experiments. For most cases $iter=100$, except for Block $T_2$, where it is set to 20.\footnote{The only estimator where there is a deviation from this is the BALQR where the variance covariance matrix of the posterior coefficients was not invertible for some of the cases. This is indicative that the BALQR prior did not shrink enough}
The bias results for the three designs (sparse, dense, block) across a selection of quantiles are presented in table ((ref)) and the results of the forecast performance are presented in table ((ref)). To shed light on whether the estimators capture the variable's quantile function appropriately, we additionally show plots for variables with non constant quantile curves for each quantile. The HS-BQR's plots are presented in figure ((ref)). The line in the plots shows the average, while the shaded region highlights the 95% coverage of $\beta$ values across the Monte Carlo runs.
Table ((ref)) shows that the HS-BQR performs competitively compared to the considered estimators in all i.i.d designs regardless of what type of sparsity structure is considered. In particular, for the sparse case the HS-BQR provides the lowest coefficient bias for both $y_1$ and $y_2$ for all quantiles. The forecast results from table ((ref)) corroborate these findings with the HS-BQR providing the lowest root mean squared forecast errors among the estimators considered.
The HS-BQR's performance is competitive for the dense and block cases as well, as can be seen in table ((ref)), however falls slightly short for the central quantiles to the BQRENET in the dense and to the BALQR in the block cases for $T_1$. Forecast errors in table ((ref)) confirm these results. This coheres with the theoretical properties of the priors. The ridge component in the BQRENET provides better inference for dense designs, while the BALQR benefits in block structures from adaptive shrinkage without having to identify a global shrinkage parameter.
Figures ((ref)) and ((ref)) show the performance of the estimators at different parts of the block design for $T_1$ and $T_2$ respectively. It reveals how the HS-BQR does extremely well in the sparse regions of the data for $y_1$ and $y_2$ while not being able to replicate this performance in the dense regions of the data for $T_1$. This is not to say that it performs poorly: while the HS-BQR yields higher average bias than the competing estimators, this is not statistically different from the bias of the other estimators. When more data are introduced in $T_2$, the difference in bias for the dence parts become even smaller among the different priors, while the sparse parts are estimated with considerably more accuracy for the HS-BQR.
Generally, as more data are introduced, the performance across the estimators converge to similar bias and forecast results, which confirms asymptotic validity of the priors and their samplers. An exception is presented by the BALQR which seems to fare worse with more data for the block design.
Both the normally distributed $y_1$ and t-distributed $y_2$ showcase a situation where the extreme quantiles (0.1 and 0.9) have higher bias than the central quantile (0.5) for all the estimators considered. This is a common finding in quantile regressions which is on account of more extreme quantiles being "data sparse" as a few observations get large weights. While it is expected that there is a U-shape in the coefficient bias as we move across the quantiles, the slope of this shape is not uniform across the estimators. In particular, it can be seen in table ((ref)) that the HS-BQR's bias does not increase as much as the other estimators.\footnote{Apart from the HS-BQR in the block design of $T_1$, where the estimators have lower coefficient bias and forecast error for its extreme low quantiles than its central quantiles.} Similarly, extreme quantiles generally tend to have higher forecast errors for all estimators, but the HS-BQR's extreme quantiles don't suffer as much as it's competition as shown in table ((ref)). This property cannot be overstated, as quantile regression is often employed for extreme quantiles. The only case where the HS-BQR's extreme quantiles performance are less accurate is for the dense design of $T_2$, where the BQRENET's performance does not suffer as much as the HS-BQR's when considering the extreme quantiles.
Figure ((ref)) underpins the findings of the tables: the HS-BQR captures the normal inverse CDF shape for $y_1$ and inverse t-distribution for $y_2$ very well in the sparse design for all $T_i$, however in the dense design, it only identifies location shift's for the more extreme quantiles for $T_1$. Nevertheless, this property is fixed when more data is available. The figure also highlights how the HS-BQR struggles the most with block designs: It only captures the quantile profiles correctly for $T_3$. This finding underpins, that in designs with unmodeled block structures and, hence, badly identified global shrinkage, quantile effects might be shrunk away. Implementation of group-level shrinkage along with prior information about the sparsity pattern in the data might be able to alleviate this problem, which we leave for future research.
As with the homoskedastic DGPs, we see that for all estimators, the error rate increases when moving away from the central quantiles and that coefficient bias as well as forecast accuracy worsens for dense and block designs compared to the sparse design. Further, the bias and forecast results in tables ((ref)) and ((ref)) show that the HS-BQR provides competitive performance to the alternative estimators, where it consistently outperforms the other priors for $y_4$ in sparse designs.\footnote{The LBQR does surprisingly poorly in the sparse $T_1$ heteroskedastic cases. This is on account of the estimator completely missing the quantile profiles for both $y_3$ and $y_4$ (see respective figures in the appendix).} Similar to the previous discussion, the HS-BQR stands out in that it provides consistently more stable inference of extreme quantiles independent of the sparsity structure, with the exception of $T_2$ dense.
In dense designs, as for the homoskedastic simulation results, the BQRENET aided by the ridge component in the prior, provides lower coefficient bias and forecast error, than the HS-BQR, whereas in block DGPs, the BALQR outperforms the HS-BQR for $y_3$ in $T_1$ but not in $T_2$. A different picture emerges for $y_4$. Here, the HS-BQR's performance in coefficient bias is only rivaled by the BQRENET for both $T_1$ and $T_2$ for both dense and block cases.
Consulting figures ((ref)) and ((ref)) for $y_3$ and $y_4$ shows how the HS-BQR performs particularly well in the sparse regions of the data, as was the case for $y_1$ and $y_2$. Just like in the homoskedastic designs, the HS-BQR is not able to replicate its performance in the dense regions of the data for $T_1$, but it does not do much worse than the competing estimators. Similarly, the HS-BQR's bias for sparse parts of the block DGP's are far smaller while for the dense parts, it's on par with the other estimators.
The plots in figure ((ref)) provide another explanation as to why the HS-BQR's forecast performance is much better for the block case of $y_4$, which is that it captures some aspects of the quantile function for $\beta_0$, even for the smallest data setting $T_1$. The plots in figure ((ref)) and ((ref)) also highlight why the estimators have lackluster performance for $y_3$ and $y_4$ for $T_1$ even for the sparse designs: The estimators have difficulties identifying the quantile profiles of $\beta_0$ and $\beta_1$ simultaneously. This deficiency is amended with more data as shown by the plots for $T_2$ and $T_3$: The HS-BQR captures the quantile profiles for both the sparse and dense DGPs, however, its performance on the block design only gets better for $T_3$. This shows the scale at which the methods require data to identify the correct quantile profiles of the variables when the DGP contains mixed sparsity structures. This shows the scale at which the methods require data to identify the correct quantile profiles of the variables when the DGP contains mixed sparsity structures.
The simulations have shown that the HS-BQR provides competitive results but also that all quantile methods under consideration have difficulty simultaneously identifying the true regressors and partialling out the location ($\beta_0(p)$) and scale ($\beta_1(p)$) effects in high dimensional setting especially when data are not abundant.
We now compare the HS-BQR to the same set of competing estimators as above in estimating forecast densities of US quarterly GDP growth as well as its down- and upside risks, commonly referred to as GaR. The need for GaR was highlighted by the global financial crisis which showed how downside risks, so the lower quantiles of the density of GDP growth, evolve with the state of credit and financial market adrian2019vulnerable,prasad2019growth. Quantifying this vulnerability is of key interest of policymakers, as it is a well-known that recessions caused by financial crises are often more severe than ordinary recessions jorda2015leveraged.
Unlike much of the previous GaR literature which focuses on GDP growth density forecasts based on only one indicator of financial distress, we apply the HS-BQR to forecasting the entire conditional GDP density using the McCracken database, a large macro economic data set. The ability to produce well calibrated density forecasts in the face of large data contexts is important for nowcasting applications, in which the information flow is necessarily high-dimensional, or variable selection of large amount of competing uncertainty indexes. The latter purpose has been suggested by adams2020forecasting and figueres2020vulnerable who have argued that is not a-priori clear which index of market frictions impacts GDP growth the most.
The mccracken2020fred database\footnote{\url{https://research.stlouisfed.org/econ/mccracken/fred-databases/}} consists of 248 variables (including GDP) from 1959 Q1 at a quarterly frequency and is updated monthly. We take the quarter-on-quarter growth rate of annualized real GDP as our dependent variable and all others as independent covariates. These variables include a wide variety of macroeconomic effects which cover real, financial as well as national accounts data. Since not all variables start at 1959 Q1, for the growth at risk application, only variables that are available from 1970 Q1 were considered which gives 229 explanatory variables.
To obtain the forecasts, we use the general linear model:
for $t=1,\cdots,T-h$, where h refers to the forecast horizon. We consider one- to four-quarter ahead forecast horizons ($h=1,\cdots,4$). Using the quantile setup, forecasts from each quantile are denoted as $y_{T+h|T}^p$. Note, that these h-step-ahead forecasts are equivalent to the h-step-ahead $p^{th}$ Value-at-Risk. Forecasts are computed on a rolling basis where the initial in-sample period uses the first 50 observations of the sample, which makes for 149-h rolling forecast windows. We estimate a grid of 19 equidistant quantiles to construct the predictive density $p(\hat{y}_{T+h|T})$ via a normal kernel \footnote{Alternatively, one could follow the popular density construction approach by adrian2019vulnerable who fit their quantiles to a skewed t-distribution. We argue when discussing the results that this approach is less flexible than the proposed approach.}.
Forecast densities are evaluated along Kolmogorov-Smirnov (KS) statistics based on (unsorted) Probability Integral Transforms (PIT) and average log-scores \footnote{There are a plethora of tests to evaluate distributions based on QQ-plot of the PIT. The choice of the KS was based solely on its simplicity to compute and any other test would suffice for evaluation.}. The PIT is often used when evaluating density forecasts and provides a measure of calibration which is independent of the econometricians loss function. In particular, the PIT is the corresponding CDF of the density function evaluated at the actual observation of the out-of-sample periods, $y_{t+h}$:
The estimated predictive density is consistent with the true density when the CDF of $g_{t+h}$ form a 45 degree line diebold1998vevaluating, i.e forms the CDF of a uniform distribution. Deviation from uniformity is tested via the Kolmogorov-Smirnov test.
Secondly we compare density fit via average log-scores. Log-scores provide a strictly proper scoring rule in the sense of gneiting2007strictly and take into account location, skewness and kurtosis of the forecast distribution gelman2013bayesian. Since quantile crossing may lead to nonsensical density forecasts, before calculating the log-scores we sort the estimated quantiles and perform kernel smoothing to obtain $p(y_{t+h})$. Average log-scores are then calculated as follows:
We break from the forecast density literature a bit, by not exclusively focusing on testing the whole density, but also evaluating specific quantiles' performance as well. To appraise the HS-BQR compared to the alternative estimators, the pseudo $R^2$ for the quantiles is computed\footnote{Since Growth-at-Risk is meant to be a VaR of growth, utilizing tests designed to test the adequacy of VaR models is a natural extension for evaluation. Two popular tests to verify the performance of a specific quantile are the DQ test of engle2004caviar and the VQR test of gaglianone2011evaluating. These tests provide a principled way of testing the null hypothesis of the selected quantile being correct. However, they do not offer a comparative measure as to how much better the proposed method provides better fit for a specific quantile.}, following koenker1999goodness. The pseudo $R^2$ of the following regression is obtained from:
where $V_{t+h}(p)$ is the fitted value of of the estimator for the $p^{th}$ quantile. Running the regression in equation ((ref)) for the $p^{th}$ quantile gives an intuitive test for the ability of the estimated fitted value to capture the dynamics we are interested in. In particular the pseudo $R^2$ is calculated the following way:
where $RASW$ is the residual absolute sum of weighted differences, so the residuals of equation ((ref)) and $TASW$ is the total absolute sum of weighted differences, so the residuals of equation ((ref)), where $\beta_1$ is constrained to 0. In essence, the pseudo $R^2$ shows how much information $V_t(p)$ adds to the regression compared to a Quantile regression with only a constant.
To gain a visual understanding of how the forecast densities perform over time, figure ((ref)) plots in its left panel the one-step-ahead forecast densities of the HS-BQR and the right panel shows all other competing estimators. The figure highlights that the HS-BQR provides better calibration especially in the beginning period of the forecast evaluation during which upper, lower and middle quantiles span a reasonable range of values despite the relative scarcity of observations to number of covariates. It is clear from the right panel that precisely in the early forecast periods, the lasso based priors offer too little regularization, yielding far too extreme upper and lower quantile growth forecasts. As more data comes in, the right panel shows that the extreme variability of the lasso based estimators decreases somewhat over time, but the HS-BQR provides good fit throughout the entire evaluation period. Forecast densities for two-three- and four- quarter ahead forecast densities are presented in section B and confirm these findings.
The visual inspection is corroborated by the more formal PIT based KS statistics and average forecast log-scores in table ((ref)): the KS statistics show that the HS-BQR is the only estimator to provide forecasts densities whose PIT are statistically indistinguishable from a uniform CDF at the 10% significance level, and whose log-scores are highest for all but the 2-quarter ahead horizon. As expected, the test statistics as well as the QQ-plots of the PITs plotted in section B, indicate that as the forecast horizon increases to 2- and 3- quarters, density calibration deteriorates somewhat for all estimators. Contrary to carriero2020nowcasting and mazzi2019nowcasting, however, we find that density fit increases again at the 4th horizon, which suggests that the HS-BQR is useful not only for short-term density forecasts, but also for medium-term forecasts.
An additional feature of the HS-BQR forecasts is that they exhibit limited quantile crossing problem. Ideally, we want the forecasted quantiles to be monotonically increasing. When this monotonicity is violated, our estimated density is invalid. The HS-BQR's forecasted quantiles exhibit very little quantile crossing, especially when comparing it to the alternative estimators. In fact, in the one-step-ahead case, the HS-BQR is the only estimator that yields non-crossing quantiles.
To quantify the relative performance of the estimators in capturing tail risks, we show in the third panel of table ((ref)) estimates of the pseudo $R^2$ which are calculated as in equation ((ref)) for the extreme and middle quantiles. It is apparent that not only does the HS-BQR provide better quantile fit at all shown quantiles, but that the largest margin (compared to the other estimators) is at the lowest and highest quantiles at all horizons which echos the results from the simulations. This is corroborated by the PIT graphs, which show that the HS-BQR's tail quantiles are consistently the closest to the ideal 45-degree line.
The proposed estimator also provides competitive point forecasts which are shown in the right panel of table ((ref)) for the 50th quantile. As shown by he1990tail, median quantile forecasts are more robust to outliers than conditional mean forecasts. Table ((ref)) clearly shows that the HS-BQR offers sizable improvements in root-mean-squared-forecast-error over the competing quantile models of 25%-66%, which are all statistically significant as per the diebold2002comparing test.
To showcase how these improvements translate to actual events of importance to policymakers, we plotted density forecasts at all horizons right before NBER marked recession or trough dates. We concentrate on the quarters before the height of each individual crisis, as the recent growth-at-risk literature highlights the usefulness of quantile methods to detect vulnerabilities to parts of the economy before these vulnerabilities materialize adrian2019vulnerable. Representative for all other pre-crisis period shown in appendix B, figure ((ref)) shows forecast densities for Q2 2008. The actual realization is marked by a vertical grey line. Two points emerge from this graph: the HS-BQR provides the largest mass at the actual realization of growth (which translates to the highest density fit for this realization as measured by the log-score) and it provides a bi-modal distribution which yields a policy relevant characterisation of forecasted risk. The second mode hovers over negative growth outcomes, thereby giving a clear indication of risks of a recession. Compared to the HS-BQR, the competing quantile methods do provide mass on negative growth outcomes which is corroborated by carriero2020nowcasting,mazzi2019nowcasting however, provide little, or close to no mass on the actual realization. In fact, consulting figure ((ref)), one can see that the lasso based BQR methods throughout the entire forecast evaluation period provide mass on negative growth outcomes, in other words forecast positive probability of recessions. This is less confidence inspiring than the forecast densities of the HS-BQR which are more conservative with mass on negative growth outcomes. To argue that this is not an artifact of the kernel smoothing, we provide forecast densities for relatively 'tranquil' economic times, namely 2005Q1, in Part B ((ref),(ref), (ref),(ref)). For these forecast densities, the HS-BQR combines to a unimodal, non-skewed, normal looking forecast density with high mass on the realization. This highlights an advantage of quantile smoothing compared to adrian2019vulnerable approach of fitting the quantiles to a t-distribution. By smoothing the 19 forecasted quantiles via a kernel we impose no restrictions on the number of modes of degree of skewness of the combined density. Finally, to compare the utility of the HS-BQR approach to a widely used forecast density constructed by survey expectations, we plotted the Survey of Professional Forecasters distribution (SPF) into the same density graphs \footnote{Smooth densities have been estimated based on a normal kernel over 19 equidistant quantiles of the survey.}. From figure ((ref)), one can see that the HS-BQR not only outperforms the SPF but provides a better indication of the looming recession indicated by larger mass on negative growth outcomes.
The ability to produce well calibrated density forecasts and, especially accurate downside risk measures in the face of large data contexts makes the HS-BQR a powerful tool for nowcasting applications or variable selection of large amount of competing uncertainty indexes.
In this paper, we have extended the widely popular horseshoe prior of carvalho2010horseshoe to the Bayesian quantile regression and provided a new algorithm to sample the shrinkage coefficients via slice sampling for the independent prior and a fast sampling algorithm that speeds up computation significantly in high dimensions.
In our simulations, we considered a variety of sparse, dense and block designs with different error distributions which revealed three points about the HS-BQR. First, the HS-BQR provides better or comparable performance in terms of both coefficient bias and forecast risk where best performance can be expected for sparse designs. Second, the aggressive shrinkage profile of the HS-BQR leads to especially good performance in tail estimation (0.1 and 0.9). Finally, an issue that all BQR methods share is simultaneously identifying the correct location and scale effects in high-dimensional setting.
Our empirical application shows that the HS-BQR provides considerable gains in calibration, density fit and even point estimates compared to double exponential based priors at all horizons, especially so at short, h=1. and medium term, h=4, horizons. Local measures of fit confirmed that HS-BQR's fitted quantiles provide the best goodness of fit. The HS-BQR proved especially useful right before NBER marked recession and trough dates, providing forecast densities foreshadowing crises. This shows that the HS-BQR is an adequate method to give credible Value-at-Risk estimates. We expect therefore that the HSBQR performs well in nowcating settings such as carriero2020nowcasting,mazzi2019nowcasting which we leave for future research.
The results show that the HS-BQR is a competitive estimator for which especially good behaviour can be expected in sparse designs with few observations. However, there are multiple fronts on which the proposed HS-BQR can be improved upon. For instance, the simulations highlighted that in dense and block designs, the HS prior tends to shrink the constant too aggressively. Hence, extensions which allow for differing shrinkage terms for subsets of the regressors might be able to alleviate this problem. Extensions to the HS-BQR should also address the problems of simultaneously estimation location and scale effects as this is needed to attain oracle properties in quantile regression.