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.
94,794 characters · 25 sections · 110 citation commands
Joint Quantile Shrinkage: A State-Space Approach toward Non-Crossing Bayesian Quantile Models
\thispagestyle{empty}
Quantile regression estimates the conditional quantile function of a response variable given a set of covariates. It is a powerful tool for inference on the relationship between response and covariates, especially in the presence of non-linearity in the covariates' impacts across the distribution of the response. While independent estimation of quantiles has been the norm in the applied literature koenker2005, the presence of non-monotonically increasing fitted quantile functions, referred to as quantile crossing, remains largely unaddressed. The probability of observing the crossing problem increases with the number of conditional quantiles estimated and the dimensionality of the covariate set wang2024composite. Many constrained optimisation solutions have been suggested in the literature, yet full probabilistic inference remains a challenge, particularly when retaining the assumption of non-parametrically modelling the error distribution. In this paper, we propose a general framework to joint quantile regression where we use the connection between a constrained multiple quantile objective function and an implied negative log posterior to motivate a novel prior structure that penalises crossing of fitted quantiles. We name the suggested prior approach the quantile-varying-parameter ($\operatorname{QVP}$) prior due to the connection to time-varying-parameters ($\mathrm{TVP}$) models popular in the state space literature. Compared to previous approaches, crossing is penalised via the structure of the prior instead of the structure of the likelihood. In particular, in this paper we:
The simulated and real data experiments confirm superior inference and prediction performance with our joint prior approach compared to commonly used Bayesian and frequentist models. For the real-world data application, we extend the methods to the Quantile vector-autoregressive (QVAR) model with Euro Area data presented in chavleishvili2024forecasting.
In section (ref) we begin by discussing previous approaches to quantile regression as background to this work. In Section (ref) we show that one can view a probabilistic generalisation of the non-crossing constraint as a prior that penalises differences across coefficients of quantiles. (ref), shows that the pseudo-likelihood can be derived from the objective function of interest. We show that this likelihood has a convenient representation as a mixture of normals. Section (ref) derives an efficient posterior sampling algorithm. Section (ref) presents an alternative representation of the $\operatorname{QVP}$ model that offers improved sampling efficiency and shrinkage properties when the data imply low amount of quantile variation. In Section (ref), we discuss post-processing methods for achieving exact sparsity of the $\operatorname{QVP}$ parameter posteriors for improved inference when sparsity in the quantile coefficient vector is suspected. In Section (ref), we investigate the theoretical shrinkage properties of the $\operatorname{QVP}$ prior. We investigate finite sample performance in Section (ref). We apply the methods presented to multivariate target data, where we estimate quantile vector autoregressive models ($\operatorname{QVAR}$) with real world data in Section (ref). We conclude in Section (ref).
yu2001bayesian established that the commonly used tick-loss quantile regression objective function implies an asymmetric-Laplace distribution ($\mathcal{ALD}$) kotz2001asymmetric as a likelihood function. Treated as a working-likelihood,\footnote{Inference on a set of quantile regression coefficients $\beta_q$ for a quantile index $q$, where the percentile $\tau_q \in (0,1)$, is asymptotically equivalent to frequentist treatment of the quantile objective, when treating the $\mathcal{ALD}$ as a working likelihood under conditions discussed in sriram2013posterior.} Bayesian inference on quantile regression models has become ubiquitous. The more recent focus being on priors for high dimensional problems kohns2024horseshoe,alhamzawi2015model,li2010bayesian. Similar to the literature on normal observation models, those priors are designed to heavily shrink coefficient of noise variables to zero polson2010shrink. Yet these approaches assume that quantile functions are independent, and therefore do not address the problem of crossing of conditional quantiles.
Many probabilistic methods have been put forward to address the issue of crossing fitted quantile functions. These generally fall into the class of semi-parametric reich2011bayesian,reich2012spatiotemporal,reich2013bayesian,kottas2001bayesian,yang2017joint, fully non-parametric scaccia2003bayesian,taddy2010bayesian, empirical likelihood lancaster2010bayesian,yang2012bayesian,yang2015quantile as well as two-step methods reich2013bayesian,rodrigues2017regression.
A common perspective taken in the semi-parametric quantile literature is to centre the quantile process for a set of finite number of percentiles on a fully parametric model which is linked piece-wise via some valid quantile function, such as the normal quantile function reich2013bayesian. This approach has some notable drawbacks. One such drawback is that the likelihood may not be available in closed form reich2011bayesian necessitating approximation methods. An additional drawback is that conditional posteriors are not available in closed form, thus prohibiting efficient updating via Gibbs MCMC methods reich2012spatiotemporal,reich2013bayesian. Computational complexity is also the bottleneck for the empirical likelihood and full non-parametric methods, even in moderate dimensions rodrigues2017regression. rodrigues2017regression propose a computationally more convenient approach, in which the first stage of estimating individual quantiles are estimated with an $\mathcal{ALD}$ likelihood. These are combined with a Gaussian process into a valid joint density in a secon step. Their approach maintains valid frequentist coverage. Yet quality of inference heavily depends on the first stage, and joint estimation in the first step has been shown to significantly improve inference even for individual quantiles bondell2010noncrossing.
Closer to our approach are the methods presented in wu2021bayesian as well as wang2024composite where quantiles are estimated jointly, and information is shared via a prior on the differences. While wu2021bayesian take a moment-based approach following chernozhukov2003mcmc, wang2024composite model the smoothness of neighbouring quantile curves via basis function expansions. In contrast, our proposed $\operatorname{QVP}$ approach connects the logic of difference penalisation of linear quantile models to the non-crossing constrained objective function presented in bondell2010noncrossing. Additionally, the $\operatorname{QVP}$ framework allows for efficient posterior computation via Gibbs sampling due to the availability of standard conditional posteriors. While we maintain the assumption that the quantile function is linear in parameters, the method can be extended to use smoothing splines as in bondell2010noncrossing.
Estimation frameworks that only allow location shifts of the quantile function are referred to as composite quantile regression ($\operatorname{CQR}$) models zou2008composite. Here only the parameter on the intercept identifies differences across quantile functions which leads to non-crossing fitted quantiles. In contrast to the above methods, the $\operatorname{QVP}$ allows for a unifying framework in which the process that models quantile variation is centred on the coefficient vector implied by the $\operatorname{CQR}$ model.
The frequentist literature has proposed many solutions to the quantile crossing problem, where one of the simplest solutions is to sort the fitted quantiles post-estimation chernozhukov2010quantile. In this paper, interest resides in both prediction and inference on quantile regression coefficients. We therefore do not consider sorting any further. Most closely related to the QVP framework in the frequentist literature is szendrei2023fused who show that, one can re-formulate the exact non-crossing objective function of bondell2010noncrossing as a fused lasso model with a particular formulation of the penalisation constant on the fused term.
The symbol $\sim$ is used both as a sampling statement as well as signifying a variable's probability density function or a likelihood function, synonymously written as $p\left(\right)$. This makes the notation close to probabilistic programming languages such as Stan carpenter_stan_2017. Denote by $\mathcal{S}$ a non-empty sample space on which the $\sigma$-algebra $\mathcal{M}$ is defined. Then, by $P\left(X|Y\right)$ we denote the probability of $X=x$ given $Y=y$, where $(X,Y) \subseteq D$. We suppress the differentiation between random and ordinary bound variables for readability. We refer to all unspecified parameters in definitions of conditional probability density functions by $\vartheta$.
Let $X = (X_1,\dotsc,X_K)^T$ be a set of covariates which are related to a response vector $y$. Let $D \subseteq \mathbbm{R}^K$ be a closed convex polytope represented as the convex hull of $\mathcal{T}$ points in $K$ dimensions. We are interested in a regression at $\mathcal{Q}$ quantile levels $0<\tau_1<\dotsc<\tau_{\mathcal{Q}}<1$, where $\mathcal{Q}$ is a finite integer. Denote by $Q_{\tau_q}(X,\beta_q)$ the $q$\textsuperscript{th} quantile of $y$ as a function of $X$, given some set of quantile-specific regression coefficients $\beta_q$, such that $P\left(y\leq Q_{\tau_q}(X,\beta_q) \vert X, \beta_q \right) = \tau_q$. The classic solution to quantile function estimation approach is to presume quantile functions to be independent, regardless of $\mathcal{Q}$. However, this neglects that the quantile coefficients, $\beta_q$, are often correlated across quantiles. Such information sharing can drastically improve inference, even if the modeller is only interested in inference on a single quantile bondell2010noncrossing,jiang2013interquantile. The key interest in this paper, is to implement information sharing across quantiles via priors on the differences of the quantile coefficients. Importantly, the quantile coefficient process is centred on a quantile invariant coefficient vector - akin to the composite quantile model zou2008composite. By doing so, as differences are adaptively shrunk to zero, the model reduces to the composite quantile regression model, which estimates parallel quantile functions. In this way the model penalises quantile crossing.
To motivate the functional form of the prior, consider the following objective function:
where $\rho_{\tau_q}(u) = u(\tau_{q}-I(u<0))$ is the tick-loss function and the constraint in Equation (ref) ensures monotonicity of the conditional quantile functions. Here, $\alpha_q\in \mathbbm{R}$ is an intercept specific to each quantile. $\beta_0 \in \mathbbm{R}^{K}$ are the quantile invariant coefficients, and $\beta_q \in \mathbbm{R}^{K}$ capture the variation in the effects of covariates across quantiles. Denote by $x_t$ the $t$-th row of X. If $\alpha$ is a monotone function of $\tau$, then the penalty term included, shrinks the model to uniformly parallel, or pair-wise parallel lines. szendrei2023fused show that when recasting the data domain to $x_t \in \left[-1,1\right]^{K}$,\footnote{Recasting the data domain to $x_t \in \left[-1,1\right]^{K}$, rather than $x_t \in \left[0,1\right]^{K}$ as in bondell2010noncrossing, has the advantage that negative ($\beta_q - \beta_{q-1}<0$) and positive ($\beta_q - \beta_{q-1}>0$) differences are treated symmetrically.} Objective Function (ref) with the constraints in Equation (ref) is equivalent to a fused lasso type quantile regression problem in which differences in $\beta$ across $q$ are penalised by the $L$-1 norm:
where $\lambda_q$ is a quantile specific shrinkage parameter, $\varsigma$ is a hyperparameter regulating the degree of “tightness” of the non-crossing constrains in szendrei2023fused and $\left\lVertx\right\rVert_j$ refers to the $j\textsuperscript{th}$ norm of x. Setting $\varsigma=1$ recovers the non-crossing constraints of bondell2010noncrossing.
In the following we offer a full probabilistic solution to this problem along with several improvements, where we allow the quantile specific penalties to be freely estimated. The imposed penalisation pushes the $\beta_q$ parameters toward the desired constrained space.
To visually motivate that this is the case, we plot in Figure (ref), the posterior means of the proposed $\mathrm{QVP}$ along-side the $\mathrm{BQR}$ model with uninformative priors which estimates quantiles independently, as well as the frequentist model of bondell2010noncrossing, $\mathrm{BRW}$. The grey area which indicates the space for which two adjacent quantile regression functions are non-crossing. We estimate the models for the median and $51^{\mathrm{st}}$ quantile. While the $\mathrm{BQR}$ leads to posterior point-estimates outside the space of non-crossing, the penalisation implied by the $\mathrm{QVP}$ , and naturally the $\mathrm{BRW}$ model, lead to coefficient estimates inside the non-crossing area. It would be expected that as quantiles are chosen to be further apart ($\Delta\alpha$ increases, volume of grey area increases), one would obtain estimates of all models in the non-crossing area. However, when the modeller wants to get an accurate a picture of coefficient-heterogeneity across quantiles, many quantiles are typically estimated which reduces a given pair of quantiles $\Delta\alpha$, thus decreasing the space in the parameter domain for which the quantile curves are non-crossing. We will show that the proposed $\mathrm{QVP}$ scales well with the amount of estimated quantiles.
Let $\vartheta$ contain parameters: $\{\alpha_q,\beta_0,\beta_q\}$. Assume prior beliefs about $\vartheta$ are represented by $p(\vartheta)$, then a valid and coherent update of $p(\vartheta)$, following the general belief updating framework of bissiri2016general, is to the posterior $p(\vartheta|X)$:
where $\ell\left(\vartheta,X\right)=\sum_{q=1}^{\mathcal{Q}}\sum_{t=1}^{\mathcal{T}} \rho_{\tau_q} \left( y_t - \alpha_q - x_t^T\beta_0 - x_t^T\beta_q\right)$.\footnote{The general updating of beliefs framework of bissiri2016general includes an unknown scalar multiplying the loss function in order to calibrate the amount of relative information from the data vis-a-vis the prior. We follow the recommendations of Bayesian quantile literature yu2001bayesian,li2010bayesian to assume that this scale is equal to 1. mclatchie2025predictive show that large enough data, the scale only marginally influences inference.} Exponentiation of the loss function is equivalent up to a constant of proportionality to the commonly employed asymmetric-Laplace ($\mathcal{ALD}$) working likelihood used for probabilistic quantile regression modelling:
Hence, minimising the expected multiple quantile loss is equivalent to maximising the combination of the individual $\mathcal{ALD}$ likelihoods of each quantile, where quantiles are assumed to be exchangeable, conditional on $\left(\alpha_q,\beta_0,\beta_q\right)$. For computational convenience, we make use of the fact that the $\mathcal{ALD}$ can be written as a mixture of normal distributions, with the scale parameter having an exponential distribution, following kozumi2011gibbs:
where $\theta_q = (1-2\tau_q)/(\tau_q(1-\tau_q))$, $\zeta_q^2 = 2/(\tau_q(1-\tau_q))$, and $\omega_{q,t}\sim \mathrm{exp}\left(\sigma^{y}_q\right)$. To vectorise across $q$, define $\boldsymbol{y} = \mathbbm{1}_{\mathcal{Q}} \otimes y$, $\boldsymbol{X} = \mathbbm{I}_{\mathcal{Q}} \otimes X$, where $\mathbbm{1}_{\mathcal{Q}}$ denotes a $\mathcal{Q}$-dimensional vector of ones, and let $\mathbbm{I}_Q$ be the identity matrix of dimension $Q\times Q$. Further, we denote the location adjustment due to data augementation by $\mu_{q,t} = \theta_q\omega_{q,t}$ . Denote $\omega_q = (\omega_{q,1},\dotsc,\omega_{q,\mathcal{T}})^T$, $\boldsymbol{\alpha} = (\alpha_1,\dotsc,\alpha_1,\dotsc,\alpha_\mathcal{Q},\dotsc,\alpha_\mathcal{Q}) \in \mathbbm{R}^{\mathcal{Q} \mathcal{T}}$, $\boldsymbol{\mu} = (\mu_{1,1},\dotsc,\mu_{1,\mathcal{T}},\dotsc,\mu_{\mathcal{Q},1},\dotsc,\mu_{\mathcal{Q},\mathcal{T}})^T$, then the joint likelihood implied is:
where $\boldsymbol{\beta} = (\beta_1^T,\dotsc,\beta_{\mathcal{Q}}^T)^T$ captures the heterogeneity induced by the covariates, and $\boldsymbol{\Omega} =\text{diag}(\omega_1,\dotsc,\omega_{\mathcal{Q}})$. $\mvn()$ stands for the multivariate normal distribution.
Consider again the penalised quantile objective of Equation (ref). A probabilistic generalisation can be found by the following discrete state space representation:
where $\beta_0$ is the quantile invariant vector and $\Sigma_q = \text{diag}(\sigma_{q,1}^2,\dotsc,\sigma_{q,K}^2)$. $\Sigma_q$ controls the variability of the coefficients between quantiles and therefore how strongly correlated the coefficients are.\footnote{Allowing for non-zero off diagonals would allow for any quantile coefficient to affect any of the other estimated quantiles directly. While relevant, we leave investigation of the properties of this modelling approach to future research.} Modelling $\epsilon_q^{\beta}$ with Laplacian densities would create the exact Bayesian equivalent to the penalisation implied by the lasso penalty in the objective function. Namely, when $\beta_q - \beta_{q-1} \sim \mathcal{MAL}\left(0, \Sigma_q \right)$, or equivilantly, $\beta_q \sim \mathcal{MAL}\left(\beta_{q-1}, \Sigma_q \right)$, where $\mathcal{MAL}\left( \right)$ stands for the multivariate Laplace distribution kotz2001asymmetric. A similar logic is presented in the derivation of the univariate lasso prior of park2008bayesian. We, however, follow the more recent literature on shrinkage priors that show superior shrinkage properties with normal kernels carvalho_handling_2009,piironen_hyperprior_2017. To adhere to the convention of the state-space literature, we will refer to Equation (ref) as the observation equation, and Equation (ref) as the state equation for quantile $\tau_q$ respectively.
The state vector $\beta_0$ plays a special role in this model setup. It is both the initialisation of the state process, and equivalent to the quantile invariant vector in Objective (ref). This can be easily verified by backward substitution of the state equation into the observation equation. In fact, when $\mathbf{\beta}=\boldsymbol{0}$, then $\beta_0$ is equivalent to the composite quantile regression vector, considered in zou2008composite. The significance of this for the shrinkage properties will be further investigated in Section (ref).
The state-space representation results in a particular structure to the joint prior. Write the state process in Equation $\ref{eq:state-equation-centred}$ stacked across $\mathcal{Q}$ in matrix form as:
where $\boldsymbol{\epsilon^{\beta}} \sim \mvn(0,\boldsymbol{\Sigma})$, $\boldsymbol{\Sigma} = \mathrm{diag}(\Sigma_1,\dotsc,\Sigma_{\mathcal{Q}})$, $\boldsymbol{\tilde{\beta}} = (\beta_0^T,0,\dotsc,0)^T$ and
$\boldsymbol{H}$ is a $\mathcal{Q}K \times \mathcal{Q}K$ difference matrix, which is invertible since $|\boldsymbol{H}| = 1$. It is straightforward to show that $\boldsymbol{H}^{-1}\boldsymbol{\tilde{\beta}} = \mathbbm{1}_{\mathcal{Q}} \otimes \beta_0$. Then, the joint prior for $\boldsymbol{\beta}$, conditional on $\left(\beta_0,\boldsymbol{\Sigma}\right)$ is
Due to the close connection to joints priors for time-varying parameter regression models in which $\beta$ is indexed by time, we call this prior the joint quantile-varying parameter ($\operatorname{QVP}$) prior. And by standard manipulations, the posterior is given by
where
and $\boldsymbol{y^{*}} = \boldsymbol{y} - \boldsymbol{\alpha} - \boldsymbol{\mu}$.
The priors on $\boldsymbol{\Sigma}$ determine the amount of quantile variation. We model explicitly: 1) adaptivity of shrinkage on the difference in coefficients across quantiles, specific to each covariate, 2) global regularisation of each $\beta_q-\beta_{q-1}$ difference vector. Therefore, it is natural to follow the global-local prior literature where the prior hierarchy on $\sigma_{q,j}$ employs a mixture of fat-tailed distributions with singularity at 0 to allow for both large and small changes in coefficients. A plethora of priors can be considered polson_half-cauchy_2012, however, we consider here the horseshoe prior carvalho_handling_2009, adapted to fused shrinkage. In particular, define $\sigma^2_{q,j} = \nu_q^2\lambda_{q,j}^2$,\footnote{Setting double-Laplace priors for $\lambda_{q,j}$, the exact Bayesian interpretation of the absolute deviation penalisation of the motivating objective function in Equation (ref) can be recovered.} then the fused horseshoe prior is
where $C_+()$ stands for the half Cauchy distribution. Following piironen_hyperprior_2017, we scale the global parameter by the number of data points. While $\nu_q$ controls overall approximate sparsity of differences in quantiles, the local scales $\lambda_{q,j}$ control the local adaptivity in how much the differences are shrunk dependent on the quantile level as well as the covariate.
We again set a horseshoe prior. Let $\Sigma_0 = \nu^2_0\text{diag}(\lambda^2_{0,1},\dotsc,\lambda^2_{0,K})$, then
This prior may also induce approximate sparsity in the quantile invariant vector.
We summarise the $\operatorname{QVP}$ prior model with the following definition:
where $p()$ stands for some probability density. With uninformative priors $\boldsymbol{\alpha} \propto 1$, any location shift in the quantile function of $\boldsymbol{\beta}$ is determined by the data only. Notice, that this deviates from incrementing the posterior with the difference prior on $alpha_q-\alpha_{q-1}$ in (ref), which would influence only the conditional posterior of the quantile specific global-shrinkage. For this and the additional reason that shrinking toward the quantile invariant vector favours shrinking toward parallel quantile curves, we opt for uninformative priors. However for completeness, we present simulation evidence in Appendix (ref) with the full difference prior on $\alpha_q-\alpha_{q-1}$. The results are virtually identical. Finally, we set the relatively uninformative prior of $\sigma^y_q \sim \text{\normalfont IG}\left(0.1,0.1\right)$, where $\text{\normalfont IG}\left(\underline{a},\underline{b}\right)$ stands for the inverse-Gamma distribution with rate $\underline{a}$ and scale $\underline{b}$.
The priors outlined above results in known conditional posterior distributions and therefore an efficient Gibbs sampling scheme which iterates through the following updates:
for $s = (1,\dotsc,S)$ until convergence. The conditional posteriors are given in Appendix (ref). The main computational bottleneck in sampling from Posterior (ref) is the inversion of the $\mathcal{Q}K \times \mathcal{Q} K$-dimensional full covariance matrix $K_\beta^{-1}$ which can easily become high-dimensional. Computing the Cholesky factor for this covariance matrix will involve $\mathcal{O}\left((\mathcal{Q} K)^3\right)$ operations. Notice that the precision matrix, $\text{K}_{\beta}$, on the other hand has a band structure, which will typically look like in Figure (ref). This motivates a more efficient sampling algorithm that utilises the sparse nature of the matrix. In particular, computing the Cholesky factor of the precision matrix only involves $\mathcal{O}\left(\mathcal{Q} K\right)$ operations, which can be sped up in practice with sparse matrix routines available in most programming languages.
Hence, to obtain draws from the conditional posterior of $\boldsymbol{\beta}$, we make use of the following steps:
where $\boldsymbol{\overline{\beta}^{(s)}}$ and $({\text{C}^{(s)}}^T)^{-1}$ can be found efficiently by solving linear equations. The rest of the sampling steps are standard and further explained in Appendix (ref).
A challenge within the multiple quantile estimation literature is the identification of quantile variation of coefficients. Viewed from the state-space representation in Equations (ref)-(ref), selection of quantile variation can be seen as a variance selection problem. This entails boundary estimation which is computationally difficult and known to lead to slow convergence of $\mathrm{MCMC}$ samplers fruhwirth2010stochastic, even with horseshoe priors bitto2019achieving. For this reason, we formulate the state-space in its non-centred form in the spirit of fruhwirth2010stochastic. This shifts the variance selection problem to a standard conjugate variable selection problem. To see this, re-write the State-equation (ref) as:
with initial condition $\tilde{\beta}_{0,j} = 0$. Using this transformation the Observation-equation (ref) can be equivalently written as:
Re-writing $x^{T}_t\text{diag}(\sigma_{q,1},\dotsc,\sigma_{q,K})\tilde{\beta}_q$ as $\tilde{x}_{q,t}^{T}\sigma_q$ where $\tilde{x}_{q,t} = \tilde{\beta}_q^{T}\text{diag}(x_t)$, the state standard deviations may be viewed simply as regression coefficients, motivating a shift of the domain of $\sigma_{q,j}$ from the positive only to the entire real line. Doing so avoids the boundary estimation issues and additionally results in conditionally conjugate posteriors, allowing for efficient Gibbs sampling. For simplicity, we employ horseshoe priors with a normal kernel for $\beta_0$ and $\sigma_q$. The prior for $\beta_0$ remains the same as in Equation (ref) and $(\sigma_1^T,\dotsc,\sigma_\mathcal{Q}^T)^T = \boldsymbol{\sigma}$ now takes the following form:
where $\tilde{\Sigma}_q = \tilde{\nu}_q^2\mathrm{diag}({\tilde{\lambda}_{q,1}}^2,\dotsc,{\tilde{\lambda}_{q,K}}^2)$ and $\boldsymbol{\tilde{\Sigma}} = \mathrm{diag}(\tilde{\Sigma}_1,\dotsc,\tilde{\Sigma}_{\mathcal{Q}})$. It can be shown, that a normal prior on the scale, $\sigma_q$, implies a generalised inverse-Gaussian prior on the variance whose properties for state-space models are studied in cadonna2020triple.\footnote{In fact, the horseshoe prior on the scale can be shown to be nested by the more general triple-gamma prior cadonna2020triple framework} This results in a higher concentration rate of the marginal prior on $\sigma_{q,j}$ near the origin and lower rate of tail-decay, which is desirable for variable selection type inference tasks polson_half-cauchy_2012.
In order to draw inference on the non-centred $\operatorname{QVP}$ model in Equation (ref), we can again make use of efficient updating via conditional posteriors:
for $s=(1,\dotsc,S)$ until convergence. Compared to the sampling steps in the centred QVP model, sampling individually the quantile invariant vector $\beta_0$ and associated hyper-parameters, adds two further sampling blocks: one for the $\boldsymbol{\sigma}$ and $\boldsymbol{\tilde{\Sigma}}$, and modified sampling steps for $\beta_0$, which we discuss in turn.
To update $\beta_0$, the relavant likelihood and prior contributions are proportional to:
where $\tilde{X}_q = X \odot \left(\mathbbm{1}_T\otimes \tilde{\beta}_q^T\right)$, and $\Omega_q = \text{diag}(\omega_q)$. Since conditional on $\beta_0$, all likelihood contributions across quantiles are exchangeable, the posterior of $\beta_0$ may be efficiently updated one quantile at a time.\footnote{This follows from basic probability theory in that $p(\theta|Y_1,Y_2)\propto p(Y_2|Y_1,\theta)\times p(\theta|Y_1)$} The conditional posterior for $\beta_0$ is thus normal:
where $K_{\beta_0} = (\sum_{q=1}^{\mathcal{Q}}X^T\boldsymbol{\Omega}^{-1}_qX) + \Sigma_0^{-1}$ and $\overline{\beta}_0 = K^{-1}_{\beta_0}(\sum_{q=1}^{\mathcal{Q}}X^T\Omega_q^{-1}(y - \alpha_q - \mu_q - \tilde{X}_q\sigma_q))$. See Appendix (ref) for further derivation of the posterior moments.
Due to the non-centred representation of the state-space, the prior for $\boldsymbol{\tilde{\beta}}$ simplifies to $\boldsymbol{\tilde{\beta}} \sim \mvn\left(0,\left(\boldsymbol{H}^T\boldsymbol{H}\right)^{-1}\right)$. Define $\boldsymbol{\check{X}} = \boldsymbol{X}\text{diag}(\boldsymbol{\sigma})$ and $\boldsymbol{\beta_0} = \mathbbm{1}_{\mathcal{Q}} \otimes \beta_0$, then the posterior for $\boldsymbol{\tilde{\beta}}$ is conditionally normal:
where $\tilde{K}_{\tilde{\beta}} = (\boldsymbol{\check{X}}^{T}\boldsymbol{\Omega}^{-1}\boldsymbol{\check{X}} + \boldsymbol{H}^T\boldsymbol{H})$ and $\boldsymbol{\tilde{\beta}} = \tilde{K}_{\tilde{\beta}}^{-1}(\boldsymbol{\check{X}}^{T}\boldsymbol{\Omega}^{-1}(\boldsymbol{Y}-\boldsymbol{\alpha} - \boldsymbol{\mu} - \boldsymbol{X}\boldsymbol{\beta}_0))$. This posterior retains its band-matrix structure as in Figure (ref), which makes for fast computation with any sparse matrix routine.
To sample the state standard deviations $\boldsymbol{\sigma}$, we can rely again on standard regression results. Let $\boldsymbol{\tilde{X}}$ contain the stacked $\tilde{X}_q$ matrices across all quantiles, then the posterior for $\boldsymbol{\sigma}$ is normal:
where $K_{\boldsymbol{\sigma}} = (\boldsymbol{\tilde{X}}^{T}\boldsymbol{\Omega}^{-1}\boldsymbol{\tilde{X}} + \boldsymbol{\tilde{\Sigma}}^{-1})$ and $\overline{\boldsymbol{\sigma}} = K_{\boldsymbol{\sigma}}^{-1} (\boldsymbol{\tilde{X}}^{T}\boldsymbol{\Omega}^{-1}(\boldsymbol{Y} - \boldsymbol{\alpha} - \boldsymbol{\mu} - \boldsymbol{X}\boldsymbol{\beta}_0))$.
Shift of the domain of $\sigma_{q,j}$ to the entire real line has advantages for computation but it introduces sign-unidentifiability. To aid mixing, we randomly permute the signs of $(\boldsymbol{\sigma},\boldsymbol{\tilde{\beta}})$ as proposed in bitto2019achieving.
Although horseshoe priors on the quantile invariant and difference vectors shrink toward sparsity, exact sparsity cannot be achieved in finite samples due to absolute continuity of the prior distributions carvalho_handling_2009. We suggest post-processing the posterior with a thresholding algorithm motivated from decision theory lindley1968choice to project the posterior onto a possibly sparse subset.\footnote{Simply calculating Bayesian p-values based on the posteriors for $\beta_0$ and $\boldsymbol{\sigma}$ might be misleading due the effect of correlation in the posteriors as well as potential multi-modality} On the one hand, this helps intuitively assessing the marginal importance of quantile invariant as well as quantile variant effects. piironen_sparsity_2017 show that post-estimation projection may also reduce variance in variable selection. On the other, it can sharpen inference of true zero effects, as well as lead to improved predictions huber2021inducing,kohns2025flexible. This literature has renewed attention with hahn2015decoupling for normal linear models which has been extended to Bayesian quantile regression in kohns2021decoupling, and feldman2023bayesian.
Define the linear predictor of interest for the quantile invariant effect as $X\beta_0$, then a possibly sparse vector, $\xi^0$, can be found by solving:
where $\upsilon^0_j$ is an adaptive penalty factor akin to adaptive lasso zou2006adaptive. Likewise, define the linear predictor that captures quantile variation as $\boldsymbol{\tilde{X}}\boldsymbol{\sigma}$, then a possibly sparse vector, $\boldsymbol{\xi^{\sigma}}$, can be found by solving:
These projections can be solved for each MCMC draw $s = 1,\dotsc,S$ in order to retrieve pseudo model-average posteriors bhattacharya2016fast. The level of sparsity in Equations (ref) and (ref) is determined by the penalties $(\upsilon^0,\upsilon^{\sigma})$. For computational convenience, the penalty term is set inversely proportional to the posterior draw of the coefficient following ray2018signal. We maintain the naming of the authors in describing this sparsification as algorithm as the $\mathrm{SAVS}$ algorithm. We refer to the $\operatorname{NC-QVP}$ with $\mathrm{SAVS}$ algorithm applied as the $\operatorname{NC-QVP}_\mathrm{SAVS}$ model.
In this section, we analyse the $\operatorname{QVP}$ prior in terms of its implications on the shrinkage scale space. It is common to analyse global-local priors in terms of their implied distribution on a scale that allows to gauge the shrinkage effect away from the maximum likelihood estimate polson_half-cauchy_2012. As is common in the literature, we analyse the shrinkage scale distribution for the normal model.\footnote{Derivations can be extended to the $\mathcal{ALD}$ using the shrinkage coefficient definitions in kohns2024horseshoe} Applied to an isotropic normal observation model, we define fused shrinkage akin to the $\operatorname{QVP}$ in the following model structure:
$K_{\beta_t} = (\nu^{-2}\Lambda_t^{-1} + \frac{1}{\sigma_y^2}x_t^Tx_t)$ and $\overline{\beta}_t = K^{-1}_{\beta_t}\left(\nu^{-2}\Lambda_t^{-1}\beta_{t-1} + \frac{1}{\sigma_y^2}x_t^Ty_t\right)$ where $\beta_t \in \mathbbm{R}^{K}$. Assume further that $x_t^Tx_t \approx \mathrm{diag}(1,\dotsc,1)$, then the posterior mean can be further decomposed as:
where $\kappa_{t,j} = \frac{1}{\nu^2\lambda_{t,j}^2\sigma^{-2}_y+1}$ and $\beta_{\mathrm{ML},t,j} = \left(x_{t,j}\right)^{-2}x_{t,j}y_t$ can be understood as a maximum likelihood esimate to the coefficient $\beta_{t,j}$. This decomposition shows that the conditional posterior mean is a convex combination of the prior coefficient $\beta_{t-1}$ and a maximum likelihood estimate. Therefore, as $\lambda^2_{t,j}\rightarrow0$, then $\kappa_{t,j}\rightarrow1$ and $(1-\kappa_{t,j})\rightarrow 0$. Hence, with strong shrinkage toward the origin, there is no further updating implied from the data, and the posterior concentrates on the previous coefficient $\beta_{t-1,j}$. On the other-hand when $\lambda^2_{t,j}\rightarrow\infty$, then $\kappa_{t,j} \rightarrow0$ and $(1-\kappa_{t,j})\rightarrow1$, so updating the posterior only happens with the data information at observation $t$, and the previous coefficient has no further influence.
See Appendix (ref) for further derivations. The implied prior distribution for $\kappa_{t,j}$ is shown in Figure (ref). Hence, shrinkage properties that are well understood for the prior on levels of coefficients also transfer to difference shrinkage. The Cauchy prior on $\lambda_{t,j}$ results in a continuous approximation to variable selection type behaviour: probability density is highest on strong or very little shrinkage respectively.
In this section we will test the performance of the proposed $\operatorname{QVP}$ priors via a simulation study following the data-generating processes of bondell2010noncrossing. We generate data from the following location scale heteroscedastic error model:
In total we consider the following 5 DGPs:
$\mathrm{DGP}$-1 and $\mathrm{DGP}$-2 are identical to the simulation study in bondell2010noncrossing. $\mathrm{DGP}$-1 features non-zero location effects of the covariates, with a shallow but continuous upward slope in the quantile coefficient profile. $\mathrm{DGP}$-2 adds sparsity to the location where true zero coefficients do not display quantile variation. $\mathrm{DGP}$-3 is defined by non-zero location effects of the covariates, where only the first three coefficients have a shallow quantile profile. $\mathrm{DGP}$-4 adds to $\mathrm{DGP}$-3 extra covariates with true zero location and quantile varying effects. Additionally, 4 covariates display what kohns2021decoupling call quantile specific sparsity patterns: location and quantile variation are zero for central quantiles, with large quantile variation in the extreme quantiles $\tau_q\geq0.9$ and $0.1\leq\tau_q$. $\mathrm{DGP}$-5 is generated with non-zero location effects, where for two covariates there is large continuous quantile variation. For ease of discussion below, we refer to coefficients of covariates with only location effects as quantile constant coefficients, to those coefficients with only a shallow profile $({\eta_1}_i = 0.1)$ as quantile varying coefficients, and those quantile varying coefficients with big jumps or large variation as extreme varying coefficients.
We simulate sample sizes $\mathcal{T}=\{100,300\}$, and estimate models for various number of total quantiles $\mathcal{Q}=\{9,19,49\}$. We do this to address the fact that the $\operatorname{QVP}$ prior's shrinkage depends on how many quantiles are being estimated and because the previous literature warns of increasing probability of crossing quantile curves with many estimated quantiles jiang2013interquantile. We allow the design matrix to be correlated with constant correlation of $\varDelta \in \{0, 0.5, 0.9\}$.\footnote{For $\varDelta>0$, the design is sampled form a normal distribution and then transformed using the min-max transformation on a column-wise basis.} We generate $N_{\mathrm{sim}} = 500$ simulated data sets. For posterior inference we obtain per model 25000 MCMC draws for 4 chains in parallel, of which we discard the first 10000 each as burnin.
Recovery of $\boldsymbol{\beta}$ is measured based on the root-mean-squared-error for the i\textsuperscript{th} DGP, by $\mathrm{RMSE}_i$:
where $\beta_{k,q}$ is the DGP's true coefficient at the $q$\textsuperscript{th} quantile and $\hat{\boldsymbol{\beta}}$ refers to the posterior mean defined as:
Denote the $q\textsuperscript{th}$ fitted quantile for observation $t$ based on the posterior mean as $\hat{Q}_{\tau_q,t}$. We evaluate predictions with the Quantile score ($\mathrm{QS}_{q,i}$) for the q\textsuperscript{th} quantile and i\textsuperscript{th} DGP:
The $\mathrm{QS}$ score belongs to the strictly proper scoring rules gneiting_strictly_2007 calculated for every simulation run based on 100 independently generated out of sample observations. The $\mathrm{QS}$ may be used to evaluate different combinations of quantiles to get an idea of the performance at different parts of the predictive density. To achieve this, we will follow gneiting2011comparing in calculating the quantile weighted $\mathrm{QS}$:
where $w_{\tau_q}$ denotes a weighting scheme. We consider four different weighting schemes: (a) $w_{\tau_q}^1=\frac{1}{Q}$ places equal weight on all quantiles, which is equivalent to taking an average of the weighted residuals; (b) $w_{\tau_q}^2={\tau_q}(1-{\tau_q})$ places more weight on central quantiles; (c) $w_{\tau_q}^3=(1-{\tau_q})^2$ places more weight on the left tail; and (d) $w_{\tau_q}^4={\tau_q}^2$ places more weight on the right tail.
Lastly, we measure the incidence of crossing. The crossing incidence is calculated by comparing the fitted quantiles with the sorted quantiles following the procedure of chernozhukov2010quantile:
where $\hat{Q}^{\mathrm{sort}}_{{\tau_q,t}}$ is the sorted predicted quantile. The crossing incidence measures the proportion of quantiles that need sorting after estimation to adhere to the property of monotone quantile functions. The lower the $\mathrm{Cross}$ value, the less quantiles need to be rearranged after estimation.
Next to the $\operatorname{QVP}$, $\operatorname{NC-QVP}$, $\operatorname{NC-QVP}_{\mathrm{SAVS}}$, we consider two approaches for independent Bayesian quantile regression, the $\operatorname{BQR}$ as presented in kozumi2011gibbs with flat priors and the $\operatorname{HSBQR}$ of kohns2024horseshoe which uses horseshoe priors on the quantile coefficients, respectively. Since both assume that the quantile functions are unrelated, we refer to these as independent quantile regression methods. We expect the independent quantile methods to do well for $\mathrm{DGP}$-4 where quantile specific sparsity is present. Additionally, to illustrate the benefit of estimating the quantile profile $\beta_q$ in addition to the composite quantile vector $\beta_0$, we also estimate a Bayesian composite quantile model which only models $\beta_0$, the $\mathrm{CQR}$ model. Lastly, since the QVP prior is motivated from the non-crossing quantile objective function of bondell2010noncrossing, we include their model too for comparison, denoted $\operatorname{BRW}$.
We summarise coefficient recovery in Figure (ref) for $\mathcal{T}=300$ and $\mathcal{Q}=19$ estimated quantiles. In Figures (ref) and (ref), we additionally show estimated coefficient profiles for a variable with quantile specific sparsity for $\mathrm{DGP}$-4 and large quantile variation for $\mathrm{DGP}$-5 respectively.\footnote{See Appendix (ref) for figures also for coefficients the other DGPs.} See Appendix (ref) for average performance for the central and extreme quantiles for the different permutations of $(T,\varDelta)$.
As expected, Figure (ref) shows that the $\operatorname{NC-QVP}$ and $\operatorname{CQR}$ models perform generally best for coefficients with low quantile variation or no quantile vairation at all (columns 1 and 2), particularly visible for DGPs 2-3. The $\mathrm{CQR}$ does well for these DGPs because estimating the true shallow quantile profile to be zero only causes small increases in $\mathrm{RMSE}$. $\operatorname{NC-QVP}$ and $\operatorname{NC-QVP}_{\mathrm{SAVS}}$ shrink low quantile variation heavily to zero, and therefore perform for these DGPs well, too. In fact the SAVS variant is on par with the $\mathrm{CQR}$ model. The $\operatorname{NC-QVP}$ offers not only a large improvement over methods that estimate the quantiles independently ($\operatorname{BQR}$,$\operatorname{HSBQR}$), but it also clearly outperforms $\mathrm{BRW}$. This is particularly visible for the quantile constant coefficients of the $\mathrm{DGP}$-3, in which the true constant coefficients are non-zero. In line with the motivation of the $\operatorname{QVP}$ prior's structure, we find that for $\mathrm{DGP}$-2-3, it generally performs on par with the $\mathrm{BRW}$.
For $\mathrm{DGP}$-4 in which there is pronounced quantile specific sparsity, we see that that models which do not explicitly model the difference in quantile coefficients ($\operatorname{BQR}$, $\operatorname{HSBQR}$, $\operatorname{CQR}$) may outperform the $\operatorname{QVP}$ models. However, Figure (ref) shows that these models also falsely shrink away the true profile to zero in the extreme tails. The $\operatorname{QVP}$ models, in contrast, correctly recover the coefficient profile, albeit with larger variance in the tails. With increasing number of observations, we would expect the reduction in posterior variance to lead also to superior $\mathrm{RMSE}$ for the $\operatorname{QVP}$ models over the independent quantile methods.
For coefficients that showcase large continuous variation across all quantiles, $\mathrm{DGP}$-4 and $\mathrm{DGP}$-5 (column 3 of Figure (ref)), we can see that recovery of the $\mathrm{QVP}$ models is competitive in terms of RMSE (Figure (ref)) and estimated coefficient profile (Figure (ref)). As expected, the $\mathrm{QVP}$ model tends to outperform the $\operatorname{NC-QVP}$ here since less shrinkage on the difference between quantiles is exerted in the centred formulation. The uncertainty bands compared to independent estimation of the quantiles, indicates large efficiency benefits to joint estimation.
These findings are robust to the number of data points (Appendix, Figure (ref)), magnitude of correlation between covariates (Appendix, Figures (ref) and (ref)), the number of quantiles estimated (Appendix, Figures (ref) and (ref)) and from the perspective of predictive performance (Appendix, Section (ref)). In fact, with more quantiles estimated, we find that there are further performance gains with the $\operatorname{QVP}$ compared to the $\operatorname{NC-QVP}$ in $\mathrm{DGP}$-5. Here, less shrinkage of the coefficient profile with the $\operatorname{QVP}$ over the $\operatorname{NC-QVP}$ allows for improved modelling in the tails with the finer resolution offered by modelling more quantiles.
Hence, in terms of parameter recovery, the $\operatorname{QVP}$ framework provides an excellent balance between the composite quantile model that only models location effects and the $\operatorname{BRW}$ model which strictly enforces non-crossing. The $\operatorname{NC-QVP}$ model in particular benefits from stronger between-quantile shrinkage when little quantile variation is present, yet also does not overtly shrink true large variation in quantile profiles. To choose in practice between the centred and non-centred formulation, we recommend using the $\operatorname{NC-QVP}$ as the baseline with which to test for the presence of no quantile variation with the $\mathrm{SAVS}$ algorithm. If significant quantile variation is present, then we recommend using the $\operatorname{QVP}$ model. \footnote{We leave investigation in terms of model selection properties for future research.} Despite the relatively low dimensionality of the covariate set, the DGPs show that the QVP methods are always preferable to independent quantile models.
The improvements in parameter recovery of the $\mathrm{QVP}$ models also translate to a lower incidence of crossing fitted quantiles as well as improved sampling efficiency. Table (ref) shows that independent of the simulation design, the $\mathrm{QVP}$ models almost completely eliminate crossing. Hence, free estimation of the implied non-crossing penalty parameter on the differences across quantiles (compare Equation (ref)), is sufficient to regularise the posterior toward the desired area of non-crossing. Independent quantile models struggle comparatively more with $\mathrm{DGP}$ 1 and $\mathrm{DGP}$ 2, where separation of quantile varying and non-varying coefficients is complicated by the relatively low amount of cross-quantile variation.
Figure (ref) shows that for all $\mathrm{DGP}$s, the $\operatorname{QVP}$ models' $\hat{R}$ and effective sample size, $\mathrm{N}_{\mathrm{eff}}$ of vehtari2021rank indicate good mixing, as well as large efficiency gains of the $\mathrm{MCMC}$ sampler over the $\operatorname{BQR}$ and $\operatorname{HSBQR}$ models. As expected, sampling is more efficient for the $\mathrm{QVP}$ model for $\mathrm{DGP}$-5 due to the large quantile variation.
For a real word data application, we apply the $\operatorname{QVP}$ prior to the quantile vector-autoregressive ($\operatorname{QVAR}$) model, as presented in chavleishvili2024forecasting. The $\operatorname{QVAR}$ generalises the autoregressive quantile model of koenker2006quantile to vector valued targets and is commonly used to examine the interactions of endogenous variables across their respective conditional distributions. $\operatorname{QVAR}$s represent an important policy tool for conducting stress-tests on financial systems, and more recently to quantify probability of tail events in macroeconomic time-series chavleishvili2023quantifying,chavleishvili2023measuring. The literature has proposed many solutions to the multiple quantile function estimation problem (see hallin2017multiple for an overview). Ambiguity arises because, unlike the univariate case, no single universally accepted mathematical framework for defining a quantile function in multiple dimensions is accepted. The approach proposed in wei2008approach is particularly convenient for economic models since the statistical identification assumption is also often defensible from the standpoint of economic theory chavleishvili2024forecasting.
The $\operatorname{QVAR}$ of order one, written $\operatorname{QVAR}(1)$, takes the following form for $t=2,\dotsc,\mathcal{T}$:
where $Y_t = (y_{1,t},\dotsc,y_{m,t})^T \in \mathbbm{R}^{m}$, $\mathcal{B}_{\tau_q} = ({b}_{1}, \dotsc, {b}_m) \in \mathbbm{R}^m $ is a vector of intercepts, $\mathcal{A}_{1,\tau_{q}} \in \mathbbm{R}^{m \times m}$ is the coefficient matrix on the lag vector, and $\mathcal{A}_{0,\tau_q} \in \mathbbm{R}^{m \times m}$ is a contemporaneous impact matrix. Denote the entry of the $i$\textsuperscript{ith} row and $j$\textsuperscript{th} column of $\mathcal{A}_{1,\tau_q}$ by $a_{i,j,1,q}$. wei2008approach show that under the assumption of a lower triangular structure to $\mathcal{A}_0$, one obtains valid multivariate quantile function estimates by estimating the system one equation at a time. With this, $\Psi_t$ is the information set at time $t$, and differs for each variable due to the lower triangular structure of $\mathcal{A}_{0,\tau_q}$: $\Psi_{1,t} = \left\{ Y_{t-1} , Y_{t-2} ,\dotsc \right\}$, $\Psi_{i,t} = \left\{ Y_{i-1,t},\Psi_{i-1,t} \right\}$ for $i = 2,\dotsc,m$.
Thus, following the steps in Sections (ref)-(ref), the probabilistic representation of the model for the $i$\textsuperscript{th} variable and $q$\textsuperscript{th} quantile of the $\operatorname{QVAR}(1)$ becomes:
where $\tilde{z}_{i,t}$ refers to the $i^{\mathrm{th}}$ information set $\Psi_{i,t}$ in vectorised form, $\text{\normalfont vec}\left(\Psi_{i,t}\right)$, and similarly $\tilde{a}_{i,q}$ vectorises the $i^{\mathrm{th}}$ row of coefficients. Note that this setup differs from the recent literature on multi-variate modelling of multiple quantiles using the multivariate $\mathcal{ALD}$ ($\mathcal{MALD}$). This is not further considered here since the $\mathcal{MALD}$ does not, without further modification, prohibit crossing of quantiles. And inference on the quantile covariance matrix is complicated due to its non-standard conditional posterior iacopini2023money.
Representing quantile function (ref) in the sample space of $Y_t$ is convenient for the subsequent forecasting and causal analysis. Define $U_t = (U_{1,t},\dotsc,U_{m,t}) \in (0,1)^m$, such that each element is distributed independently as a uniform distribution. wei2008approach show that if the joint distribution of $Y_t$ is absolutely continuous, there exists a one-to-one mapping between the sample space of $Y$ and the hyper-cube $\left(0,1 \right)^m$\footnote{This is known as the Rosenblatt transformation.}:
where, one can obtain the standard $\operatorname{VAR}$ like representation of the model by making the right hand side, a function of a constant and lags only:
where $v(U_t) \equiv \left( \mathbbm{I}_m - \mathcal{A}_0(U_t) \right)^{-1}\mathcal{B}(U_t)$ and $\mathcal{C}(U_t) \equiv \left( \mathbbm{I}_m - \mathcal{A}_0(U_t) \right)^{-1}\mathcal{A}_1(U_t)$. Set in Equation (ref) $U_t = \tau_q$ for $t=2,\dotsc,\mathcal{T}$.\footnote{These inverses can be shown to always exists whenever the $\operatorname{QVAR}$ coefficients imply a stationary process, which is equivalent to the stationary conditions on parameter matrices of standard $\operatorname{VAR}$ models.}
The data set for this application, obtained from chavleishvili2024forecasting on Euro Area macrofinancial data, contains the industrial production growth-rate $(\operatorname{IP})$, which measures real economic activity, and the composite indicator of systemic stress $(\operatorname{CISS})$, representing a measure of financial health of the Euro Area.\footnote{For more examples of the $\operatorname{CISS}$ used in Growth-at-Risk models, see figueres2020vulnerable, szendrei2023revisiting or varga2025non among others.} The time-series are plotted in Figure (ref). Data are monthly and available from January 1999 to July 2018. The goal of the study is to jointly forecast $\operatorname{IP}$ and $\operatorname{CISS}$ with the $\operatorname{QVAR}$ system and perform causal analysis of how the conditional distribution of $\operatorname{IP}$ responds to perturbations in the $\operatorname{CISS}$ indicator. chavleishvili2024forecasting give economic justification to a lower-triangular identification scheme where $\operatorname{IP}$ impacts $\operatorname{CISS}$ contemporaneously, but $\operatorname{IP}$ only impacts the $\operatorname{CISS}$ with a lag.\footnote{This identification assumption is implicit in adrian2019vulnerable and many more studies that followed. Lower‐triangular identification entails that shocks are identified through how they propagate dynamically through the system of equations. This approach is widely used in macroeconomics because it yields a unique decomposition of structural shocks for location scale $\operatorname{VAR}$ models, often matching intuitive causal narratives sims1980macroeconomics. }
We estimate the $\operatorname{QVAR}$ with the same set of priors as in Section (ref). We generate h-step-ahead quantile predictions $Y_{t+h}$ for $h = 1,3,6$, which are evaluated with the $\operatorname{QS}$ score presented in Section (ref). Then, we perform causal analysis common to the $\operatorname{VAR}$ applications in which we estimate the impulse response function of $\operatorname{IP}$ to a one-standard-deviation shock in the $\operatorname{CISS}$.
To retrieve quantile predictions, we follow the procedure given in chavleishvili2024forecasting, however using the MCMC draws where appropriate. We summarise this in Algorithm (ref).
Forecasts are produced on an expanding in-sample time-window $t = 1\dotsc,\mathcal{T}_r$, with initial in-sample period $\mathcal{T}_{\mathrm{start}}=96$ and $\mathcal{T}_{\mathrm{end}}=224$. Quantile weighted $\operatorname{QS}$ scores for the $\operatorname{QVAR}$s are shown in Table (ref). Similar to the results in Section (ref), we show forecasting performance relative to the $\operatorname{BQR}$-$\operatorname{QVAR}$ model.
Several clear patterns emerge. First, joint‐quantile models consistently outperform independent quantile models ($\operatorname{BQR}$, $\operatorname{HSBQR}$), confirming the benefit of `partially pooling' information across quantiles observed in Section $\ref{sec:simulation}$. Second, while $\operatorname{HSBQR}$ often excels at central quantiles, its tail performance is much worse than the $\operatorname{QVP}$ models, especially at longer forecast horizons. This is particularly visible for the $\operatorname{IP}$ variable. Third, the $\operatorname{QVP}$ prior variants can even outperform the $\operatorname{BRW}$ model, particularly in the tails. Here, integration over the parameter space with the Bayesian models induces more smoothness of the quantile function of the coefficients (see Section (ref)) and therefore reduces variance in the tails. Finally, forecast accuracy gets worse as the horizon increases which is a consequence of model parsimony and accumulation of forecast uncertainty. However, the loss in forecasts accuracy is less pronounced for joint estimation frameworks, particularly for the $\operatorname{IP}$ equation.
In terms of crossing of the estimated (in-sample) quantiles, we find, similar to the Section (ref), that the $\operatorname{QVP}$ priors almost completely eradicate the issue of crossing quantile curves (see Table (ref)).
$\operatorname{QIRF}$s measure the distributional causal effect of an unexpected change in $Y_{i,t}$ on $Y_{j,t}$ over some horizon $h = 1,\dotsc,W$. Such analyses are of interest to policy institutions like central banks, who tailor their policy instruments to the likelihood the potential paths real economic output take in relation to movements of the financial sector. Common to location-scale $\operatorname{VAR}$ analysis, the $\operatorname{QVAR}$ approach allows the analysis of such dynamic responses along particular points of the conditional distributions, in particular the tails.
Denote by $Y_t^{*}$ a hypothetical `shocked' vector $Y_t$, to which $\iota \in \mathbbm{R}^m$ is added. As in chavleishvili2024forecasting, we define $\iota$ as a vector of zeros, except entry $i$ which is equal to the shock's magnitude. The quantile impulse-response function at time point $t$ is defined as $\delta_{t}(U_t) \equiv Y_t^{*} - Y_t$. chavleishvili2024forecasting show that under the triangular identification scheme from above, this reduces to $\mathcal{D}(U_t)\iota$, where $\mathcal{D}(U_t) = \left( \mathbbm{I}_m - \mathcal{A}_0(U_t) \right)^{-1}$. The impulse response function for $h$-steps ahead is then given by: $\delta_{t+h}\left( U_{t+h} \vert U_t,\dotsc, U_{t+h-1}\right) = \prod_{g=1}^h\mathcal{C}(U_{t+g})\mathcal{D}(U_t)\iota$. We define the quantile response levels of interest of $\operatorname{IP}$ to be equal to $U_{\operatorname{IP},t+g} = \tau = (0.05,0.1,\dotsc,0.95) \forall g$ while that of $\operatorname{CISS}$ are fixed to the median level, $U_{\operatorname{CISS},t+g} = 0.5 \forall g$. Note, however, that for estimation of the coefficient-posteriors, it remains that both equations are estimated for all $\tau$.\footnote{We keep the quantile levels for the $\operatorname{QIRF}$ constant respective to each equation of the $\operatorname{QVAR}$. In principle, the quantile indices must stay constant across the horizons of the $\operatorname{QIRF}$.} The initial shock $\iota$ is set equal to a one standard deviation of the residuals at the median quantile of the $\operatorname{CISS}$ equation chavleishvili2024forecasting. The previous literature finds that shocks to financial conditions have a more pronounced negative impact on the left tail of real economic output adrian2019vulnerable,chavleishvili2024forecasting. It is therefore expected that the $\operatorname{QIRF}$s show a pronounced negative impact at low quantile levels that even out to 0 at high quantile levels.
Figure (ref) shows the impulse response functions for the various quantile models for $h=1,\dotsc,30$.
As expected, Figure (ref) shows generally that all models estimate a pronounced negative impact of a shock to $\operatorname{CISS}$ to the left tail of $\operatorname{IP}$ with impulse response functions petering out to zero as the quantile level increases. While the $\operatorname{BQR}$ and $\operatorname{HSBQR}$ models exhibit a notably sharp dip in the lower quantile only, the $\operatorname{QVP}$ models estimate a more gentle slope along the quantile levels. This is due to the $\operatorname{QVP}$ priors leading to smoother posterior quantile coefficient profiles - even compared to the $\operatorname{BRW}$ model.\footnote{Posteriors of the coefficients are shown in the appendix in figures (ref), (ref), and (ref).}
Compared to chavleishvili2024forecasting, we find that joint estimation of quantiles leads to significant differences in the $\operatorname{QIRF}$s, particularly at the median. Figure (ref) shows the $\operatorname{QIRF}$ at selected quantiles for better visibility (including posterior uncertainty intervals). While the $\operatorname{BQR}$ model, in line with chavleishvili2024forecasting, predicts no significant impact of the $\operatorname{CISS}$ shock at the median, the $\operatorname{QVP}$ model predicts a persistent negative one. Taken together, we confirm the previous literature's finding that post shock distributions of $\operatorname{IP}$ exhibit negative skewness, however, we also observe a significant downward location shift identified by the $\operatorname{QVP}$ models. This highlights how using joint estimation, along with a prior that regulates information `pooling', can have a significant influence on the inference drawn from these models.
In this paper, we defined a prior for multiple quantile regression in which information across quantiles is shared via an adaptive joint shrinkage prior. The structure of the prior is motivated from the penalised non-crossing objective function from bondell2010noncrossing, and is shown to imply a quantile state-space representation, named $\operatorname{QVP}$ model, where unknown states are equal to the quantile regression coefficients. This allows for the derivation of efficient sampling methods where the resultant triangular structure of the conditional posterior precision allows for fast computation. We extend the $\operatorname{QVP}$ framework to a non-centred formulation ($\operatorname{NC-QVP}$) as well as a post-estimation sparsification algorithm that allow for stronger shrinkage on state variability and sparsity, respectively. With this method, we were able to tackle the issue of quantile crossing through a structured prior that regularises toward the desired parameter sub-space, rather than modifying the likelihood.
A simulation exercise shows that the $\operatorname{QVP}$ priors result in far superior predictive performance and parameter recovery compared to Bayesian methods that estimate quantiles independently. Additionally, crossing is nearly completely eliminated with this approach. For low true quantile variation, the $\operatorname{NC-QVP}$ models can also offer large gains over frequentist methods that strictly enforce non-crossing.
In the empirical application of a $\operatorname{QVAR}$ on the Euro Area, following chavleishvili2024forecasting, we show the practical advantages of the $\operatorname{QVP}$ prior in modelling complex macroeconometric dynamics. We produce quantile forecasts as well as conduct a causal study of the effect of financial shocks to the distribution of industrial production, $\operatorname{IP}$. $\operatorname{QVP}$ models produce very competitive forecasts, often outperforming all models under comparison.
For the causal study, we generate impulse response functions ($\operatorname{QIRF}$) of $\operatorname{IP}$ in reaction to shocks to worsening financial conditions. We verify the finding that financial shocks exert markedly asymmetric and persistent effects across the conditional distribution of $\operatorname{IP}$. The $\operatorname{QVP}$ priors produce smoother $\operatorname{QIRF}$s with larger negative effects at the lower tails which are more persistent.
Despite these advantages, there are several avenues by which the method can be improved. First, the paper focuses on implementing the framework to linear quantile regression models. However, the method can be extended to nonlinear settings as well. This extension can enhance the applicability of the $\operatorname{QVP}$ prior framework. Second, we have exclusively focused on the horseshoe prior for modelling the differences across quantiles. The $\operatorname{QVP}$ framework can be used on various other types of shrinkage priors, such as the GIGG prior (see for example kohns2025flexible).
We acknowledge the computational resources provided by the Aalto Science-IT project, and the support of the Research Council of Finland Flagship programme: Finnish Center for Artificial Intelligence, Research Council of Finland project (340721), and the Finnish Foundation for Technology Promotion. We thank Gary Koop, Niko Hauzenberger, Ping Wu, Aristeidis Raftapostolos, and all the participants of the 2024 CFE conference, 2025 RSS conference for their feedback. We further thank the Scotland national rugby union teams for their continued effort both on and off the field. \onecolumn