Bayesian Smoothed Quantile Regression

• Bayesian Smoothed Quantile Regression (BSQR) is a family of methods that smooths the check loss or coefficient trajectories to produce differentiable likelihoods and stable quantile estimates.
• BSQR employs kernel-smoothed losses, continuous shrinkage priors, and Gaussian-process adjustments to address the limitations of traditional asymmetric Laplace-based quantile regression.
• Empirical applications of BSQR highlight enhanced effective sample sizes, reduced bias at extreme quantiles, and robust, coherent inference in various complex data scenarios.

Searching arXiv for Bayesian Smoothed Quantile Regression and closely related Bayesian quantile smoothing papers. Bayesian Smoothed Quantile Regression (BSQR) denotes a family of Bayesian quantile-regression procedures in which some form of smoothing is introduced into quantile inference, posterior computation, or cross-quantile regularization. In the most explicit recent usage, BSQR replaces the non-differentiable check loss of standard quantile regression with a kernel-smoothed loss, producing a continuously differentiable likelihood that supports Hamiltonian Monte Carlo and yields a consistent posterior for the target quantile parameter (Liu et al., 3 Aug 2025). In related literature, however, “smoothed” has also referred to continuous shrinkage priors for sparse quantile regression, Gaussian-process adjustment across quantile levels to enforce noncrossing, trend-filtering priors on discrete differences of a latent quantile curve, Bayesian empirical-likelihood coupling across quantiles, and simultaneous spline-based modeling of the full quantile process (Li et al., 2021).

1. Conceptual scope and relation to standard Bayesian quantile regression

The common starting point is the linear conditional quantile model

$Q_Y(\tau\mid x)=x^\prime\beta(\tau),$

with the classical estimator obtained by minimizing the check loss

$$\sum_{i=1}^n \rho_\tau(y_i-x_i^\prime \beta), \qquad \rho_\tau(u)=u\big(\tau-\mathbb{I}(u<0)\big).$$

Standard Bayesian quantile regression typically replaces this optimization problem with an asymmetric Laplace distribution (ALD) working likelihood,

$$f(y \mid \mu,\sigma,\tau) = \frac{\tau(1-\tau)}{\sigma} \exp\left\{ -\rho_\tau\!\left(\frac{y-\mu}{\sigma}\right) \right\},$$

where $\mu=x^\prime\beta(\tau)$ targets the $\tau$-th conditional quantile (Santos et al., 2016).

This ALD formulation is foundational for much of the BSQR literature, but it also motivates smoothing. One line of work emphasizes that the classical quantile estimator does not adapt to sparsity and that its asymptotic variance depends on the unknown conditional density at zero, which must otherwise be estimated nonparametrically (Li et al., 2021). Another line argues that ALD-based Bayesian quantile regression suffers from two fundamental limitations: the non-differentiability of the check loss precludes gradient-based MCMC methods, and the posterior mean provides biased quantile estimates (Liu et al., 3 Aug 2025).

The scope of BSQR is therefore broader than a single algorithm. In some papers it means smoothing the loss; in others it means smoothing across quantile levels, smoothing coefficient behavior through continuous shrinkage, or smoothing a latent quantile trend. This suggests that BSQR is best understood as a family of Bayesian quantile methods whose defining feature is regularized, cross-quantile, or differentiable quantile inference rather than a single universally fixed hierarchy.

2. Kernel-smoothed BSQR as a differentiable likelihood framework

The most direct formulation of BSQR replaces the check loss $\rho_\tau$ with a kernel-smoothed version,

$L_h(e;\tau) = (\rho_\tau * K_h)(e),$

where $K_h(v)=h^{-1}K(v/h)$ and $h>0$ is a bandwidth. The resulting smoothed empirical objective is

$$\widehat R_h(\boldsymbol{b};\tau,h) = \frac1n\sum_{i=1}^n L_h(e_i(\boldsymbol{b});\tau), \qquad e_i(\boldsymbol{b})=y_i-\boldsymbol{x}_i^\prime\boldsymbol{b}.$$

Its score is

$$\sum_{i=1}^n \rho_\tau(y_i-x_i^\prime \beta), \qquad \rho_\tau(u)=u\big(\tau-\mathbb{I}(u<0)\big).$$0

with the closed-form identity

$$\sum_{i=1}^n \rho_\tau(y_i-x_i^\prime \beta), \qquad \rho_\tau(u)=u\big(\tau-\mathbb{I}(u<0)\big).$$1

so the gradient is analytically available (Liu et al., 3 Aug 2025).

BSQR then defines a smoothed quantile error density

$$\sum_{i=1}^n \rho_\tau(y_i-x_i^\prime \beta), \qquad \rho_\tau(u)=u\big(\tau-\mathbb{I}(u<0)\big).$$2

with normalizing constant

$$\sum_{i=1}^n \rho_\tau(y_i-x_i^\prime \beta), \qquad \rho_\tau(u)=u\big(\tau-\mathbb{I}(u<0)\big).$$3

For regression data, the likelihood becomes

$$\sum_{i=1}^n \rho_\tau(y_i-x_i^\prime \beta), \qquad \rho_\tau(u)=u\big(\tau-\mathbb{I}(u<0)\big).$$4

and the posterior is proportional to this likelihood times $$\sum_{i=1}^n \rho_\tau(y_i-x_i^\prime \beta), \qquad \rho_\tau(u)=u\big(\tau-\mathbb{I}(u<0)\big).$$5 (Liu et al., 3 Aug 2025).

The theoretical motivation is not merely computational. Under regularity conditions, the posterior concentrates around the true quantile parameter $$\sum_{i=1}^n \rho_\tau(y_i-x_i^\prime \beta), \qquad \rho_\tau(u)=u\big(\tau-\mathbb{I}(u<0)\big).$$6, thereby addressing the inconsistency and bias issue attributed to the ALD-based posterior mean. The same paper establishes posterior propriety under improper uniform priors on $$\sum_{i=1}^n \rho_\tau(y_i-x_i^\prime \beta), \qquad \rho_\tau(u)=u\big(\tau-\mathbb{I}(u<0)\big).$$7, Gaussian priors, and hierarchical Gaussian priors, subject to conditions involving $$\sum_{i=1}^n \rho_\tau(y_i-x_i^\prime \beta), \qquad \rho_\tau(u)=u\big(\tau-\mathbb{I}(u<0)\big).$$8 and the prior on $$\sum_{i=1}^n \rho_\tau(y_i-x_i^\prime \beta), \qquad \rho_\tau(u)=u\big(\tau-\mathbb{I}(u<0)\big).$$9. It also characterizes kernel-dependent curvature through the Hessian

$$f(y \mid \mu,\sigma,\tau) = \frac{\tau(1-\tau)}{\sigma} \exp\left\{ -\rho_\tau\!\left(\frac{y-\mu}{\sigma}\right) \right\},$$0

showing that more peaked kernels produce a more concentrated posterior locally. Compact-support kernels such as the uniform and triangular kernels are emphasized because they preserve ALD-like tail behavior while remaining computationally efficient (Liu et al., 3 Aug 2025).

3. Meanings of “smoothing” in the broader BSQR literature

The term “smoothing” is used in several distinct but related ways across Bayesian quantile-regression research.

Framework	Smoothing object	Main purpose
"Posterior Inference for Quantile Regression: Adaptation to Sparsity" (Li et al., 2021)	Continuous shrinkage on coefficients	Adaptation to sparsity without hard selection
"Bayesian empirical likelihood for quantile regression" (Yang et al., 2012)	Prior coupling across quantiles	Borrow strength across percentile levels
"Regression Adjustment for Noncrossing Bayesian Quantile Regression" (Rodrigues et al., 2015)	Gaussian-process smoothing over $$f(y \mid \mu,\sigma,\tau) = \frac{\tau(1-\tau)}{\sigma} \exp\left\{ -\rho_\tau\!\left(\frac{y-\mu}{\sigma}\right) \right\},$$1	Noncrossing and cross-quantile borrowing
"Fast and Locally Adaptive Bayesian Quantile Smoothing using Calibrated Variational Approximations" (Onizuka et al., 2022)	Shrinkage on discrete differences of a latent trend	Locally adaptive quantile trend estimation
"Bayesian Non-parametric Simultaneous Quantile Regression for Complete and Grid Data" (Das et al., 2016)	Tensor-product B-spline expansion of $$f(y \mid \mu,\sigma,\tau) = \frac{\tau(1-\tau)}{\sigma} \exp\left\{ -\rho_\tau\!\left(\frac{y-\mu}{\sigma}\right) \right\},$$2 or $$f(y \mid \mu,\sigma,\tau) = \frac{\tau(1-\tau)}{\sigma} \exp\left\{ -\rho_\tau\!\left(\frac{y-\mu}{\sigma}\right) \right\},$$3	Simultaneous noncrossing estimation of the whole quantile function
"Bayesian Quantile Regression with Subset Selection: A Decision Analysis Perspective" (Feldman et al., 2023)	Smoothing induced by a single Bayesian model for $$f(y \mid \mu,\sigma,\tau) = \frac{\tau(1-\tau)}{\sigma} \exp\left\{ -\rho_\tau\!\left(\frac{y-\mu}{\sigma}\right) \right\},$$4	Quantile extraction, uncertainty quantification, and subset selection from one posterior

In sparse linear quantile regression, the “smoothed” aspect is explicitly identified with continuous shrinkage priors that separate active and inactive coefficients without a discontinuous model-selection step. The adaptive lasso-type prior

$$f(y \mid \mu,\sigma,\tau) = \frac{\tau(1-\tau)}{\sigma} \exp\left\{ -\rho_\tau\!\left(\frac{y-\mu}{\sigma}\right) \right\},$$5

and the clipped-absolute prior

$$f(y \mid \mu,\sigma,\tau) = \frac{\tau(1-\tau)}{\sigma} \exp\left\{ -\rho_\tau\!\left(\frac{y-\mu}{\sigma}\right) \right\},$$6

induce sparsity through continuous posterior shrinkage rather than dichotomous variable selection (Li et al., 2021).

In cross-quantile smoothing, the empirical-likelihood literature uses informative and shrinking priors on $$f(y \mid \mu,\sigma,\tau) = \frac{\tau(1-\tau)}{\sigma} \exp\left\{ -\rho_\tau\!\left(\frac{y-\mu}{\sigma}\right) \right\},$$7 to encode commonality across quantiles, particularly in data-sparse tail regions (Yang et al., 2012). The Gaussian-process adjustment literature fits separate Bayesian quantile models at a grid of quantile levels and then smooths the induced quantiles over the quantile index $$f(y \mid \mu,\sigma,\tau) = \frac{\tau(1-\tau)}{\sigma} \exp\left\{ -\rho_\tau\!\left(\frac{y-\mu}{\sigma}\right) \right\},$$8 using a squared-exponential kernel

$$f(y \mid \mu,\sigma,\tau) = \frac{\tau(1-\tau)}{\sigma} \exp\left\{ -\rho_\tau\!\left(\frac{y-\mu}{\sigma}\right) \right\},$$9

thereby borrowing strength from neighboring quantiles and monotonizing the fitted quantile function (Rodrigues et al., 2015).

In nonparametric smoothing, simultaneous quantile regression expands the entire conditional quantile process as

$\mu=x^\prime\beta(\tau)$0

with the coefficient functions $\mu=x^\prime\beta(\tau)$1 themselves expanded in tensor-product B-splines over predictors. This imposes monotonicity and prevents quantile crossing by construction (Das et al., 2016).

4. Computational architectures

A major computational backbone for Bayesian quantile regression is the location-scale mixture representation of the ALD. In one standard form,

$\mu=x^\prime\beta(\tau)$2

with

$\mu=x^\prime\beta(\tau)$3

and latent $\mu=x^\prime\beta(\tau)$4 exponentially distributed a priori with mean $\mu=x^\prime\beta(\tau)$5. In regression form, this yields normal full conditionals for $\mu=x^\prime\beta(\tau)$6, inverse-gamma updating for $\mu=x^\prime\beta(\tau)$7 under an inverse-gamma prior, and generalized inverse Gaussian full conditionals for the latent variables $\mu=x^\prime\beta(\tau)$8 (Santos et al., 2016).

Kernel-smoothed BSQR changes this architecture. Because the loss is continuously differentiable, the posterior gradient with respect to $\mu=x^\prime\beta(\tau)$9 exists in closed form, and the paper uses a hybrid Metropolis-within-Gibbs strategy: Hamiltonian Monte Carlo or NUTS for $\tau$0, and Metropolis-Hastings for the scale parameter $\tau$1, since the normalizing constant $\tau$2 complicates direct gradient-based updates for $\tau$3 (Liu et al., 3 Aug 2025).

Quantile trend-filtering work combines the ALD mixture representation with mean-field variational Bayes. The factorization

$\tau$4

produces Gaussian, generalized inverse Gaussian, and inverse-gamma variational factors, followed by a residual-bootstrap calibration step

$\tau$5

to restore nominal frequentist coverage of credible intervals (Onizuka et al., 2022).

The Gaussian-process adjustment approach is modular: first-stage ALD-based Bayesian quantile models are fit independently at multiple quantiles, and the second stage has a closed-form GP posterior for the adjusted quantile. A further reduction shows that fitting the GP only to induced posterior means gives the same adjusted mean and variance as using all $\tau$6 MCMC draws, which substantially lowers the second-stage computational burden (Rodrigues et al., 2015).

Bayesian empirical likelihood uses Metropolis-Hastings directly on

$\tau$7

because evaluating the empirical likelihood ratio $\tau$8 is manageable even though global maximization is difficult (Yang et al., 2012). Simultaneous spline-based methods instead employ Block Metropolis-Hastings over simplex-constrained spline-coefficient blocks, with a warm start obtained by the derivative-free optimization algorithm GCDVSMS (Das et al., 2016).

5. Statistical properties and inferential goals

The inferential goals attached to BSQR differ across formulations, but several recurrent themes are visible: sparsity adaptation, noncrossing, locally adaptive regularization, improved uncertainty quantification, and valid posterior asymptotics.

In sparse quantile regression, continuous shrinkage priors yield two-rate posterior contraction under $\tau$9: active coefficients contract at the usual $\rho_\tau$0 rate, while inactive coefficients contract at $\rho_\tau$1. The posterior factorizes asymptotically into a Gaussian oracle component for active coefficients and a shrinkage-dominated component for inactive ones. The resulting posterior mean has oracle asymptotic efficiency for the active coefficients and super-efficiency for the inactive coefficients, while a variance-corrected posterior interval uses an extra factor $\rho_\tau$2 to account for shrinkage asymmetry (Li et al., 2021).

In Bayesian empirical likelihood, the posterior is asymptotically normal for fixed priors, centered at the maximum empirical likelihood estimator, and its variance approaches that of the MELE. Under shrinking priors that favor similar slopes across quantiles, the posterior center becomes a precision-weighted compromise between prior information and empirical-likelihood information. A common-slope corollary shows that, when the true slopes are equal across several quantiles, this shrinkage can reduce variance without introducing first-order bias at the scale relevant for inference (Yang et al., 2012).

In Gaussian-process adjustment for noncrossing quantile regression, monotonicity is obtained by bandwidth control. For any finite grid $\rho_\tau$3, the paper states that there exists a bandwidth $\rho_\tau$4 such that the adjusted estimates are ordered. The practical rule is to choose the smallest $\rho_\tau$5 that eliminates crossing everywhere, creating a continuum between separate Bayesian quantile fits and a fully smoothed noncrossing quantile surface (Rodrigues et al., 2015).

In locally adaptive quantile smoothing, shrinkage on $\rho_\tau$6-st differences of a latent trend allows the posterior to favor piecewise-polynomial structure while preserving abrupt jumps under Laplace-type or horseshoe-type priors. Because mean-field variational Bayes tends to produce overly narrow intervals, the calibrated variational framework inflates posterior covariance until empirical coverage matches the target level, especially at extreme quantiles (Onizuka et al., 2022).

In the formal kernel-smoothed BSQR framework, posterior consistency is established under regularity conditions on the error density, covariates, kernel, and prior support. The normalizing constant $\rho_\tau$7 is shown to be strictly decreasing in $\rho_\tau$8, with $\rho_\tau$9 convex in $L_h(e;\tau) = (\rho_\tau * K_h)(e),$0, and compact-support kernels yield BSQR posteriors that are sandwiched between constant multiples of the ALD posterior. This preserves similar tail behavior while modifying the local geometry of the likelihood near zero residuals (Liu et al., 3 Aug 2025).

6. Empirical behavior and application domains

The empirical literature presents BSQR-style methods as especially useful when quantile estimation is difficult because of sparsity, extreme $L_h(e;\tau) = (\rho_\tau * K_h)(e),$1, nonstationary structure, or incoherence across separately fitted quantiles.

In the kernel-smoothed BSQR formulation, simulations report up to a 50% reduction in predictive check loss at extreme quantiles relative to ALD-based methods and 20–40% improvements in effective sample size, with compact-support kernels—especially uniform and triangular—emerging as the most effective trade-offs between concentration and computation. In a rolling quantile-CAPM application to JPMorgan Chase returns against the S&P 500 during 2017–2025, BSQR produces smoother and more stable $L_h(e;\tau) = (\rho_\tau * K_h)(e),$2 and $L_h(e;\tau) = (\rho_\tau * K_h)(e),$3 paths than ALD-based Bayesian quantile regression, and the average minimum ESS is roughly 80% higher for BSQR in that application (Liu et al., 3 Aug 2025).

Locally adaptive Bayesian quantile smoothing via trend filtering is designed for nonstationary quantile trends, piecewise-constant and piecewise-linear signals, varying smoothness, and misspecified noise. The calibrated variational procedure is reported to be computationally much more efficient than the Gibbs sampler and to provide stable inference results, especially for high and low quantiles (Onizuka et al., 2022).

Related Bayesian quantile diagnostics also remain relevant to BSQR’s broader aims. In ALD-based Bayesian quantile regression, posterior learning of the scale parameter $L_h(e;\tau) = (\rho_\tau * K_h)(e),$4 is argued to be more principled than fixing it, because the ALD variance structure is quantile dependent and $L_h(e;\tau) = (\rho_\tau * K_h)(e),$5 is U-shaped. The posterior distributions of the latent variables $L_h(e;\tau) = (\rho_\tau * K_h)(e),$6 can also be used to detect quantile-specific outliers. In an application to Gini indexes in Brazilian states, the Federal District observations stand out in the lower tail and Santa Catarina in 2010 stands out in the upper tail, illustrating that outlyingness can itself be quantile dependent (Santos et al., 2016).

The application range of BSQR-style ideas is broad. Decision-theoretic posterior quantile projection has been used to produce quantile-specific subset selection and variable importance in educational outcomes in North Carolina, with gains in quantile estimation accuracy, inference, and variable selection over several Bayesian and frequentist competitors (Feldman et al., 2023). Bayesian quantile models with structured smoothing have been used for risk margin estimation in non-life insurance, including AL-based dynamic mean, variance, and skewness formulations for high-quantile reserve and capital analysis (Dong et al., 2014). For longitudinal count data, an artificial smoothing of counts through jittering and transformation enables ALD-based Bayesian quantile regression with random effects and Gibbs sampling, extending quantile-specific inference beyond conditional means (Jantre, 2022). Simultaneous spline-based quantile regression has been applied to North Atlantic hurricane intensity and US household income, where full noncrossing quantile surfaces reveal temporal changes in the upper tail that would be difficult to capture with separately estimated quantiles (Das et al., 2016).

A persistent misconception is that BSQR is simply synonymous with Bayesian quantile regression using the ALD. The literature does not support that equivalence. Some papers treat ALD-based Bayesian quantile regression as the baseline to be corrected, smoothed, or regularized; others use ALD only as a computational surrogate; and some avoid a parametric likelihood altogether through empirical likelihood or posterior decision analysis. Another misconception is that smoothing always refers to smoothing over covariates. In fact, the smoothing target may be the loss function, the coefficient process over $L_h(e;\tau) = (\rho_\tau * K_h)(e),$7, the discrete differences of a latent trend, the quantile function over the quantile index, or the posterior summaries extracted from a single Bayesian model. The resulting field is therefore heterogeneous, but its unifying concern is stable, coherent, and computationally tractable Bayesian inference for conditional quantiles.