Bayesian Smoothed Quantile Regression

Bayesian Smoothed Quantile Regression

\doublespacing

Introduction

Quantile regression (QR) Koenker1978 constitutes a cornerstone of modern econometrics, enabling robust characterization of distributional heterogeneity essential for financial risk management and policy evaluation Chernozhukov2013. Rigorous specification testing and inference for quantiles have been extensively studied in the frequentist literature Linton2005,Whang2006. In the Bayesian paradigm, QR is prized for coherent uncertainty quantification and flexibility in hierarchical and latent-variable modeling. In practice it is used for tail risk, inequality, and heterogeneous treatment effects, often with short time spans or moderately high-dimensional covariates. In such settings we find that standard Bayesian quantile regression (BQR) can be systematically misleading: Bayesian quantile regression based on the asymmetric Laplace distribution (BQR-ALD) produces severely biased and under-covering posteriors at extreme quantiles in small samples, penalized exponentially tilted empirical likelihood (PETEL) Tang2022 exhibits poor mixing once $p$ becomes even moderately large, and in systemic risk applications conventional BQR-ALD delivers unstable dynamic betas with erratic credible bands.

Traditionally, BQR relies on the ALD working likelihood Yu2001 and its scale-mixture representation for Gibbs sampling Kozumi2011. This device makes BQR practically accessible, but it faces two fundamental bottlenecks. First, as emphasized by \shortciteA{Gneiting2011} and \shortciteA{Sriram2013}, the ALD-based posterior is decision-theoretically misaligned: while asymptotically consistent, the posterior mean minimizes expected squared error rather than check loss, violating proper scoring rules and inducing severe finite-sample bias, especially in the tails. Second, the non-smooth check loss precludes state-of-the-art gradient-based Markov chain Monte Carlo (MCMC) algorithms such as Hamiltonian Monte Carlo (HMC) Neal2011, effectively confining practitioners to conditionally Gaussian augmentation schemes whose mixing deteriorates rapidly with dimension, hierarchical structure, or extreme quantiles.

Recent work seeks to move beyond these limitations via generalized Bayesian inference (GBI) Bissiri2016 and joint quantile estimation frameworks Yang2017. PETEL Tang2022 delivers a moment-based posterior but inherits the geometric fragility of empirical likelihood: posterior support can be empty or nearly empty in finite samples, mixing in high dimensions is often extremely poor, and extensions to hierarchical or latent-variable models are structurally cumbersome. Likelihood-free approaches such as the martingale posterior Fong2025 instead treat the loss as a scoring rule and optimize predictive performance without an explicit parametric likelihood. While powerful for forecasting, they obscure the direct covariate--quantile link that underpins structural interpretation, marginal effects, and policy counterfactuals. In short, existing approaches either retain the non-smooth check-loss geometry and forgo gradient-based computation, impose fragile empirical-likelihood constraints that scale poorly, or abandon parametric likelihoods and thereby weaken structural interpretability.

We propose Bayesian smoothed quantile regression (BSQR), a principled framework that bridges frequentist smoothing and robust Bayesian inference. Building on the “conquer” convolution-smoothing strategy Fernandes2021,He2023 and the quasi-Bayesian paradigm Chernozhukov2003, we construct a principled, fully differentiable pseudo-likelihood whose negative log-density coincides with a kernel-smoothed check loss. This formulation explicitly manages the classical bias--variance trade-off in smoothing Gozalo2000, restores decision-theoretic coherence with the underlying quantile loss, defines a valid error distribution on $\mathbb{R}$ that avoids convex-hull pathologies, and unlocks the substantial efficiency gains of HMC and the No-U-Turn sampler (NUTS) Hoffman2014 for BQR.

Our contributions advance the literature along three dimensions: methodological realignment, theoretical rigor, and computational scalability. Methodologically, we exploit a symmetrizing property of convolution smoothing to rectify the decision-theoretic flaw of standard BQR. Using the smoothed quantile regression (SQR) objective of \shortciteA{He2023}, we define a new error distribution whose negative log-density is proportional to the smoothed check loss and show that the pseudo-true parameter coincides with the minimizer of expected check loss. This aligns the posterior mean with the conditional quantile decision problem and eliminates the finite-sample bias induced by squared-error optimality in BQR-ALD. The smoothed loss admits closed-form first and second derivatives in residuals, enabling efficient gradient-based sampling for regression and scale parameters. The resulting pseudo-likelihood is strictly positive on $\mathbb{R}$, ensuring global posterior existence and avoiding empty-support phenomena that can afflict empirical-likelihood-based methods. In simulations this rectification reduces out-of-sample check loss by up to 40--50% relative to BQR-ALD and markedly improves tail behavior and coverage, especially for extreme quantiles and short samples.

Theoretically, we go beyond simple consistency and provide a comprehensive asymptotic and geometric analysis. We establish posterior consistency for the true linear conditional quantile under mild regularity and an undersmoothing regime for the bandwidth, despite using a smoothed pseudo-likelihood. We then prove a Bernstein--von Mises (BvM) theorem under possible misspecification, leveraging the $C^2$-smoothness of the loss to obtain a classical local asymptotic normality expansion and an explicit limiting covariance for the centered posterior, thereby justifying credible intervals for frequentist coverage. Complementing this, we derive a generalized Wilks-type result: by calibrating the scale parameter to satisfy a generalized information equality, BSQR credible sets achieve asymptotically correct frequentist coverage even under misspecification; both plug-in and fully Bayesian treatments of this scale are analyzed and shown to enhance robustness. Furthermore, we analyze posterior propriety under various priors and derive geometric insights linking kernel “peakedness” to posterior concentration, establishing tail equivalence with ALD posteriors for compact kernels and clarifying the bias--variance and robustness trade-offs inherent in smoothing.

Computationally, this work marks, to our knowledge, the first systematic deployment of gradient-based sampling for BQR with an explicit likelihood and BvM justification. Smoothing the objective renders the posterior jointly differentiable in regression and scale parameters, enabling joint HMC/NUTS updates in standard probabilistic programming platforms such as Stan Carpenter2017. We show that the normalizing constant of the smoothed likelihood is log-convex in the scale parameter, implying a log-concave conditional posterior under standard priors and numerically stable gradient evaluations. Extensive simulations and a systemic risk application demonstrate that BSQR substantially reduces out-of-sample prediction error by up to 50% relative to BQR-ALD and significantly enhances sampler efficiency (up to 80% in our applications), maintaining robust coverage in extreme tails. In higher-dimensional settings (e.g., $p=50$), it achieves effective sample sizes per second that are an order of magnitude larger than those of PETEL, whose chains often suffer near-complete mixing collapse due to the geometric fragility of empirical likelihood constraints. In our systemic risk study, BSQR produces smoother and more persistent dynamic betas with tighter yet well-calibrated credible bands, while matching or improving the forecasting performance of BQR-ALD.

Collectively, these innovations position BSQR as a rigorous bridge between frequentist accuracy and Bayesian probabilistic richness, retaining the structural interpretability prized in applied econometrics while delivering decision-theoretically aligned inference, classical large-sample guarantees, and modern computational efficiency. The complete source code is publicly available for reproducibility.\footnote{The replication package is available at \url{https://github.com/BeauquinLau/BSQR}.}

The paper is organized as follows. Section (ref) presents the BSQR model and smoothed likelihood. Sections (ref), (ref), and (ref) establish the theoretical core, covering posterior consistency, the Bernstein--von Mises theorem, and posterior propriety. Section (ref) analyzes kernel effects and guides kernel and bandwidth selection. Section (ref) details the HMC/NUTS computational framework. Sections (ref) and (ref) present simulation evidence and an empirical application to asymmetric systemic risk. Section (ref) concludes. All proofs, algorithms, technical derivations, and simulation tables are provided in the appendices.

The Bayesian smoothed quantile regression framework

This section develops the BSQR framework. We first review the fundamental limitations of standard BQR to motivate our approach, then lay out the technical foundations of loss smoothing, and finally construct our principled Bayesian model upon this foundation.

Motivation: The challenge of non-smoothness in Bayesian quantile inference

The linear QR model posits that for i.i.d. observations $(y_i, \boldsymbol{x}_i)$ from a joint distribution $(Y, \boldsymbol{X})$, the $\tau$-th conditional quantile of the response $Y$ is a linear function of covariates for a given $\tau \in (0,1)$:

where $Q_Y(\tau \mid \boldsymbol{x}_i) \coloneqq \inf \{q : F_{Y \mid \boldsymbol{X}}(q \mid \boldsymbol{x}_i) \geq \tau \}$ is defined by the conditional cumulative distribution function (CDF) $F_{Y \mid \boldsymbol{X}}$, $\boldsymbol{x}_i \in \mathbb{R}^d$ is a covariate vector including an intercept, and $\boldsymbol{\beta}(\tau) \in \mathbb{R}^d$ is the coefficient vector. This implies the error term $\varepsilon_i \coloneqq y_i - \boldsymbol{x}_i^\top\boldsymbol{\beta}(\tau)$ has a zero $\tau$-th conditional quantile, $Q_{\varepsilon}(\tau \mid \boldsymbol{x}_i) = 0$. The assumption of independence between $\varepsilon_i$ and $\boldsymbol{x}_i$ simplifies this to the unconditional requirement $F_{\varepsilon}(0) = \tau$.

The population coefficient vector $\boldsymbol{\beta}(\tau)$ is the minimizer of the expected check loss Koenker1978:

where $\zeta \coloneqq Y - \boldsymbol{X}^\top\boldsymbol{b}$ is the population-level residual for a candidate vector $\boldsymbol{b}$, $F_{\zeta}$ is its CDF, and $\rho_{\tau}(e) = e(\tau - \mathbb{I}(e<0))$ is the non-differentiable “pinball” loss. Correspondingly, the sample estimator $\hat{\boldsymbol{\beta}}(\tau)$ minimizes the empirical risk:

where $e_i(\boldsymbol{b}) \coloneqq y_i - \boldsymbol{x}_i^\top\boldsymbol{b}$.

The BQR paradigm connects to this objective via the ALD, whose probability density function (PDF) is $p_{\mathrm{ALD}}(\cdot; \theta, \tau) \propto \exp(-\theta\rho_\tau(\cdot))$ Yu2001. While this formulation ensures that its posterior mode, which we denote $\check{\boldsymbol{\beta}}(\tau)$, numerically coincides with the frequentist point estimator $\hat{\boldsymbol{\beta}}(\tau)$, this reliance on the non-smooth check loss imposes critical limitations. First, it precludes the use of modern gradient-based samplers (e.g., HMC). Second, it creates an inferential misalignment: the posterior mean is a biased estimator of the true conditional quantile and does not minimize the expected check loss, the canonical strictly proper scoring rule for quantile forecasts Gneiting2011, Sriram2013. These challenges motivate a fundamental reformulation of the Bayesian likelihood.

Technical foundation: Kernel smoothing of the check loss

To address the non-differentiability, we adapt kernel smoothing\,---\,a technique with a rich history in econometrics for enhancing estimation efficiency and utilizing parametric information Gozalo2000\,---\,to construct a differentiable parametric likelihood. Specifically, we adopt the convolution-type smoothing strategy recently advanced in the frequentist literature, particularly the “conquer” framework established by \shortciteA{Fernandes2021} and \shortciteA{He2023}. While these works utilized smoothing primarily to facilitate gradient-based M-estimation and bootstrap inference, we leverage their formulation as the structural basis for our Bayesian likelihood.

Following \shortciteA{He2023}, the smoothed check loss is constructed by replacing the indicator function in the standard check loss with a smooth approximation derived from a kernel. Formally, let $K(\cdot)$ be a symmetric, non-negative kernel density function integrating to one, and let $K_h(v) = h^{-1}K(v/h)$ be the scaled kernel with bandwidth $h>0$. The smoothed check loss function (see Figure (ref)) is defined as the convolution of the standard check loss $\rho_{\tau}(\cdot)$ and the kernel $K_h(\cdot)$:

Based on this smoothed loss, the smoothed quantile regression (SQR) objective function is defined as the average smoothed loss over the sample. The resulting SQR estimator, denoted as $\hat{\boldsymbol{\beta}}_h(\tau)$, is:

A distinct advantage of this formulation, which we exploit for HMC sampling, is its differentiability. The derivative of the smoothed loss, denoted as $\Psi_h(e; \tau)$, is the convolution of the check loss subderivative $\psi_\tau(e) = \tau - \mathbb{I}(e<0)$\footnote{Strictly speaking, $\psi_\tau(e) = \tau - \mathbb{I}(e<0)$ is a subderivative, as $\rho_\tau(e)$ is non-differentiable at $e=0$ due to the discontinuity in $\psi_\tau$. However, this non-differentiability does not affect the subsequent smoothing and computation, as the convolution with a sufficiently smooth kernel $K_h$ yields a differentiable $\Psi_h(\cdot;\tau)$.} with the kernel:

\footnotetext{The interchange of differentiation and integration is justified by the dominated convergence theorem. The partial derivative of the integrand with respect to $e$, namely $\psi_\tau(e-v)K_h(v)$, is bounded in absolute value by $\max(\tau, 1-\tau)K_h(v)$. This dominating function is integrable over $v \in \mathbb{R}$ since $K_h$ is a probability density, thus satisfying the conditions for the theorem.}By evaluating this convolution integral, we explicitly derive a closed-form expression involving the kernel's CDF, $F_K(\cdot)$:

This analytical gradient is fundamental for efficient HMC-based Bayesian inference, as detailed in Section (ref). Furthermore, to elucidate the geometry of the posterior distribution sampled by HMC, we explicitly derive the second derivative (the curvature). Differentiating the smoothed score function $\Psi_h(e; \tau)$ from Eq. (ref) with respect to $e$:

The derivative of $\psi_\tau(u) = \tau - \mathbb{I}(u<0)$ is $\psi_\tau'(u) = \delta(u)$, where $\delta(u)$ is the Dirac delta function.\footnote{Strictly speaking, the classical derivative of $\psi_\tau(u)$ does not exist at $u=0$ due to the discontinuity. However, in the distributional sense, it is given by $\psi_\tau'(u) = \delta(u)$. To see this, note that we can write $\psi_\tau(u) = \tau - 1 + \mathbb{I}(u \geq 0)$. The distributional derivative of the indicator function $\mathbb{I}(u \geq 0)$ is the Dirac delta $\delta(u)$, leading to $\psi_\tau'(u) = \delta(u)$. This does not affect the convolution, which yields a smooth second derivative.} Thus, we obtain

by utilizing the sifting property of the Dirac delta function, which states $\int_{-\infty}^{\infty} f(x)\delta(x-a)\,\mathrm{d} x = f(a)$.

Consequently, the Hessian of the sum of smoothed losses (the SQR objective) with respect to $\boldsymbol{\beta}$ takes the form of a weighted covariance matrix:

This explicit Hessian structure serves a dual purpose in our framework: it characterizes the posterior curvature essential for constructing the mass matrix in HMC, ensuring efficient exploration of the posterior parameter space, and it is central to our asymptotic analysis in Section (ref) as well as the discussion of kernel effects in Section (ref).

The BSQR likelihood and posterior

While previous works utilized the smoothed loss $L_h(\cdot;\tau)$ primarily as a computational surrogate for fast optimization or bootstrap inference He2023, constructing a principled Bayesian framework requires a valid probability density. Since the standard ALD-based likelihood relies on the non-smooth check loss $\rho_\tau(\cdot)$ and generally $L_h(\cdot;\tau) \neq \rho_\tau(\cdot)$, the ALD posterior does not target the smoothed estimator $\hat{\boldsymbol{\beta}}_h(\tau)$. To achieve both computational efficiency and inferential coherence, we formulate a new likelihood grounded in the smoothing geometry.

We propose a novel error distribution, denoted as $f_{\text{SQR}}$, constructed explicitly with a negative log-density proportional to the smoothed check loss. Unlike the ALD, our distribution is smooth and explicitly depends on the bandwidth $h$:

where $\theta > 0$ is a scale parameter. A critical component here, which is absent in the optimization-based “conquer” framework, is the normalizing constant $Z(\theta, \tau, h)$ required for $f_{\text{SQR}}$ to be a valid PDF:

This constant ensures that $f_{\text{SQR}}$ is a proper density, a feature essential for Bayesian model selection and the valid sampling of $\theta$.

For an observed dataset $\boldsymbol{y} = (y_1, \dots, y_n)^\top$ and $\boldsymbol{\mathcal{X}} = (\boldsymbol{x}_1^\top, \dots, \boldsymbol{x}_n^\top)^\top$, under the assumption of independence, the joint likelihood for parameters $\boldsymbol{\beta}$ and $\theta$ is:

The log-likelihood is thus:

By construction, maximizing this likelihood with respect to $\boldsymbol{\beta}$ is equivalent to minimizing the SQR objective (Eq. (ref)), establishing a principled link between our Bayesian model and the frequentist SQR estimator. With independent priors $\pi(\boldsymbol{\beta}, \theta) = \pi(\boldsymbol{\beta})\pi(\theta)$, the posterior is:

This formulation overcomes the limitations of the standard framework by providing: (1) differentiability for efficient HMC sampling (Section (ref)); (2) a consistent estimator that resolves the inferential bias of BQR-ALD (Section (ref)); and (3) a theoretically valid basis for uncertainty quantification, backed by Bernstein-von Mises theorems (Section (ref)) and rigorous propriety conditions (Section (ref)).

Posterior consistency and bias resolution

Posterior consistency is foundational. While standard ALD-based BQR is asymptotically consistent Sriram2013, its posterior mean suffers from critical “inferential misalignment” in finite samples by minimizing squared error rather than the target check loss Gneiting2011. Although joint estimation frameworks like \shortciteA{Tang2022} and \shortciteA{Yang2017} address structural issues, they typically entail heavy computational burdens and struggle with mixing in high-dimensional settings, impeding hierarchical extensions. In contrast, BSQR resolves this misalignment directly through loss smoothing. This section establishes the consistency of the BSQR posterior under the standard assumption of a correctly specified model, demonstrating that its symmetrizing property ensures inferences are correctly centered and unbiased. We extend this analysis to allow for model misspecification in Section (ref).

The proof of this theorem is provided in \hyperref[pf:consistency]{Appendix A.1}.

Theorem (ref) provides the theoretical cornerstone for BSQR. It demonstrates that our model retains the foundational property of posterior consistency while simultaneously resolving the computational and finite-sample limitations that plague the standard BQR-ALD. The proof hinges on the symmetrizing property of our smoothing procedure: convolving the asymmetric check loss with a symmetric kernel yields an objective function that is asymptotically unbiased in expectation, as shown by the pivotal result that $\mathbb{E}[\Psi_{h}(\varepsilon_{0i}; \tau)] \to 0$. This property ensures that, in large samples, BSQR inferences\,---\,such as posterior means and credible intervals\,---\,are correctly centered on the true parameters, thereby providing a foundation for the superior out-of-sample performance observed in our simulations (Section (ref)), particularly at extreme quantiles where the finite-sample bias of BQR-ALD is most severe.

Limiting distribution and frequentist validity

While Theorem (ref) establishes posterior consistency, rigorous inference requires characterizing the limiting posterior shape and ensuring its validity even under potential model misspecification\,---\,a central theme in recent GBI literature Matsubara2022, Tang2022. Unlike standard BQR, where the non-smooth check loss necessitates complex empirical process theory, BSQR's kernel-smoothed loss $L_h(\cdot; \tau)$ is twice continuously differentiable (guaranteed by the continuity of $K(\cdot)$ in Assumption A4). This regularity allows us to establish a BvM theorem via a classical local asymptotic normality (LAN) expansion, addressing both posterior existence and frequentist coverage.

Global existence and well-posedness

A limitation of standard empirical likelihood approaches is that the posterior support is restricted to the convex hull of the estimating equations. While recent methods such as PETEL Tang2022 introduce penalization to alleviate this emptiness problem, BSQR ensures global validity naturally without requiring auxiliary penalty terms or relaxation of moment conditions.

The detailed proof is provided in \hyperref[pf:existence]{Appendix A.2}.

Rate of convergence and bias control

Before establishing the asymptotic normality of the posterior, it is essential to characterize the convergence rate of the frequentist smoothed estimator $\hat{\boldsymbol{\beta}}_h(\tau)$. This rate determines the allowable range of bandwidths to ensure that the smoothing bias does not contaminate the asymptotic distribution.

The detailed proof is provided in \hyperref[pf:rate]{Appendix A.2}.

Bernstein-von Mises theorem under misspecification

We now derive the asymptotic normality of the posterior. While frequentist approaches often employ specification tests to validate quantile models Whang2006, our Bayesian framework rigorously accounts for potential misspecification via the generalized Bernstein-von Mises theorem. To ensure robustness, we do not assume the linear quantile model is correctly specified. Instead, following \shortciteA{Matsubara2022}, we define the target parameter $\boldsymbol{\beta}^*_h(\tau)$ as the minimizer of the expected smoothed risk (the pseudo-truth): \[ \boldsymbol{\beta}^*_h(\tau) \coloneqq \arg\min_{\boldsymbol{b}} \mathbb{E}[L_h(e_i(\boldsymbol{b}); \tau)]. \] Note that the smoothed loss $L_h(\cdot; \tau)$ is convex and converges pointwise to the check loss $\rho_\tau(\cdot)$ as $h \to 0$. Consequently, assuming the linear specification is valid (i.e., the true conditional quantile satisfies $Q_{Y}(\tau \mid \boldsymbol{x}_i) = \boldsymbol{x}_i^\top \boldsymbol{\beta}_0(\tau)$), standard M-estimation theory for convex objective functions ensures that the pseudo-truth $\boldsymbol{\beta}^*_h(\tau)$ converges to the true parameter $\boldsymbol{\beta}_0(\tau)$ Horowitz1998a. We define the population Hessian and score covariance matrices at this pseudo-truth:

The proof exploits the $C^2$ smoothness of $L_h(\cdot; \tau)$ to apply a Taylor expansion of the log-posterior, avoiding the empirical process theory required for non-smooth losses (see \hyperref[pf:bvm]{Appendix A.2}).

However, valid frequentist inference for $\boldsymbol{\beta}_0(\tau)$ faces two impediments. First, under model misspecification, standard quasi-Bayesian posteriors exhibit covariance mismatch Tang2022: the posterior covariance $(\theta \mathcal{H}(\boldsymbol{\beta}^*_h(\tau)))^{-1}$ generally diverges from the true “sandwich” variance $\mathcal{H}(\boldsymbol{\beta}^*_h(\tau))^{-1} \mathcal{J}(\boldsymbol{\beta}^*_h(\tau)) \mathcal{H}(\boldsymbol{\beta}^*_h(\tau))^{-1}$. Departing from the empirical likelihood remedy of \shortciteA{Tang2022}, we calibrate the scale parameter of the smoothed objective directly to align these variances, following the sandwich adjustment strategy of \shortciteA{Mueller2013}. Second, smoothing introduces deterministic bias. The following corollary establishes that combining an undersmoothing schedule (controlling bias) with a generalized information equality condition (calibrating variance) yields valid asymptotic frequentist coverage.

The detailed proof is provided in \hyperref[pf:cor_validity_true]{Appendix A.2}.

\noindentExample (Analytic calibration in standard cases). To build intuition for the abstract condition in Eq. (ref), consider the standard i.i.d. linear model defined in Eq. (ref), where errors satisfy $f_{\varepsilon}(0)>0$. In this scenario, as the bandwidth $h \to 0$, the smoothed population matrices approach the classical limits established in standard quantile regression theory Koenker2005: \[ \mathcal{H}(\boldsymbol{\beta}^*_h(\tau)) \to f_{\varepsilon}(0) \Sigma_X \quad \text{and} \quad \mathcal{J}(\boldsymbol{\beta}^*_h(\tau)) \to \tau(1-\tau) \Sigma_X, \] where $\Sigma_X = \mathbb{E}[\boldsymbol{x}\boldsymbol{x}^\top]$. We apply the trace calibration formula from Eq. (ref) to these limits. First, we compute the product matrix inside the trace: \[ \mathcal{J}(\boldsymbol{\beta}^*_h(\tau)) \left[\mathcal{H}(\boldsymbol{\beta}^*_h(\tau))\right]^{-1} \to (\tau(1-\tau)\Sigma_X) (f_{\varepsilon}(0)\Sigma_X)^{-1} = \frac{\tau(1-\tau)}{f_{\varepsilon}(0)} \boldsymbol{I}_d, \] where $\boldsymbol{I}_d$ is the $d \times d$ identity matrix. Its trace is simply $d \cdot \frac{\tau(1-\tau)}{f_{\varepsilon}(0)}$. Next, substituting this trace value into the definition of the target parameter yields: \[ \theta^* = \frac{d}{\mathrm{tr}\left( \mathcal{J}(\boldsymbol{\beta}^*_h(\tau)) \left[\mathcal{H}(\boldsymbol{\beta}^*_h(\tau))\right]^{-1} \right)} = \frac{d}{d \cdot \frac{\tau(1-\tau)}{f_{\varepsilon}(0)}} = \frac{f_{\varepsilon}(0)}{\tau(1-\tau)}. \] This result explicitly recovers the standard scalar correction factor found in classical quantile regression theory. It demonstrates that our generalized trace calibration target $\theta^*$ naturally simplifies to the correct theoretical scaling in well-specified linear models. Consequently, the plug-in estimator $\hat{\theta}$ defined in Footnote (ref) will consistently estimate this classical factor, ensuring valid coverage without requiring prior knowledge of the error density. The geometric intuition behind this variance rescaling is illustrated in Figure (ref).

Posterior propriety under various prior specifications

This section investigates posterior propriety for the BSQR model under common prior choices for the regression coefficients $\boldsymbol{\beta}$ and the scale parameter $\theta$. We analyze cases ranging from an improper uniform prior for $\boldsymbol{\beta}$ to hierarchical structures. Our analysis uses the BSQR likelihood from Eq. (ref) and the shorthand $S(\boldsymbol{\beta}; \tau, h) \coloneqq \sum_{i=1}^n L_h(e_i(\boldsymbol{\beta}); \tau)$ for the sum of smoothed losses.

Although Corollary (ref) establishes asymptotic validity based on a consistent estimator $\hat{\theta}$, finite-sample Bayesian practice treats $\theta$ as a random parameter to capture scale uncertainty. Establishing posterior propriety under this specification is therefore essential to validate the joint Markov chain Monte Carlo (MCMC) sampling (Section (ref)). We structure this analysis in three progressive steps: first deriving baseline conditions under an improper uniform prior (Section (ref)), extending them to a standard Gaussian prior (Theorem (ref)), and finally validating the hierarchical Gaussian specification (Corollary (ref)) used in our implementation.

Propriety with an improper uniform prior for \texorpdfstring{$\boldsymbol{\beta}$}{beta}

We begin by analyzing posterior propriety under an improper uniform prior, $\pi(\boldsymbol{\beta}) \propto 1$. This setting is particularly insightful, as it isolates the likelihood's contribution to integrability from the regularizing effects of a proper prior for $\boldsymbol{\beta}$. We first establish baseline conditions for a fixed $\theta$ and then generalize to settings where $\theta$ has its own prior. Theorem (ref), proved in \hyperref[pf:thm_propriety1]{Appendix A.3}, summarizes these results.

Propriety with proper priors for \texorpdfstring{$\boldsymbol{\beta}$}{beta}

We now transition to proper priors for $\boldsymbol{\beta}$, which can enhance regularization and simplify propriety arguments. A conventional choice, seen in many BQR methodologies Kozumi2011,Li2010, is a Gaussian prior, often with a large variance to create a weakly informative yet mathematically convenient specification. Theorem (ref), proved in \hyperref[pf:thm_propriety2]{Appendix A.3}, establishes the propriety conditions under this prior.

A natural extension treats the prior variance $\sigma_{\boldsymbol{\beta}}^2$ as random, allowing the data to inform its scale while increasing flexibility. Assigning $\sigma_{\boldsymbol{\beta}}^2$ an Inverse-Gamma hyperprior, favored for its positive support and tractability, leads to the propriety result in the following corollary. The detailed proof is provided in \hyperref[pf:cor_1]{Appendix A.3}.

Theoretical analysis of kernel selection

Having established the conditions for posterior propriety in Section (ref), we now investigate how the choice of the kernel function $K(\cdot)$ shapes the specific geometry of the BSQR posterior. We first establish that for any compact support kernel, the BSQR posterior for $\boldsymbol{\beta}$ is equivalent in its tail behavior to that of standard ALD-based BQR. Second, we show that more “peaked” kernels yield a more concentrated posterior. These results provide a theoretical justification for our method and a principle for kernel selection.

Implications of kernel selection

Building on Section (ref), we examine how the choice of kernel $K(\cdot)$ influences the BSQR posterior. Two properties are central: (i) the behavior of the normalizing constant $Z(\theta, \tau, h)$, which governs both the likelihood and MCMC sampling of $\theta$; and (ii) the relationship between the BSQR and standard ALD posteriors for compact support kernels, clarifying the role of smoothing in inference on $\boldsymbol{\beta}$.

The proof of this proposition is provided in \hyperref[pf:prop1]{Appendix A.4}.

A pertinent question is how the BSQR posterior, $\pi_{\mathrm{BSQR}}(\boldsymbol{\beta} \mid \theta, \boldsymbol{y}, \boldsymbol{\mathcal{X}};\tau,h)$, relates to the ALD-based posterior, $\pi_{\mathrm{ALD}}(\boldsymbol{\beta} \mid \theta, \boldsymbol{y}, \boldsymbol{\mathcal{X}}; \tau)$, particularly in their tail behaviors. For compact-support kernels (e.g., Uniform, Epanechnikov, Triangular, with support normalized to $[-1,1]$), the smoothed loss $L_h(e; \tau)$ coincides with the check loss $\rho_\tau(e)$\,---\,up to an additive constant independent of $e$\,---\,whenever $|e/h| > 1$, thereby localizing smoothing to residuals near zero. This motivates a formal investigation into whether the two posteriors are equivalent under such kernels.

We define $L_{\mathrm{ALD}}(\boldsymbol{y} \mid \boldsymbol{\mathcal{X}}, \boldsymbol{\beta}, \theta; \tau) \propto \exp\left( -\theta \sum_{i=1}^n \rho_\tau(e_i(\boldsymbol{\beta})) \right)$, with a normalizing constant for the ALD PDF independent of $\boldsymbol{\beta}$. The corresponding posterior is $\pi_{\mathrm{ALD}}(\boldsymbol{\beta} \mid \theta, \boldsymbol{y}, \boldsymbol{\mathcal{X}};\tau) \propto L_{\mathrm{ALD}}(\boldsymbol{y} \mid \boldsymbol{\mathcal{X}}, \boldsymbol{\beta}, \theta; \tau) \pi(\boldsymbol{\beta})$. The following theorem, proved in \hyperref[pf:thm4]{Appendix A.4}, establishes this equivalence.

Kernel effects on posterior concentration

Beyond tail behavior equivalence under compact-support kernels, the choice of kernel $K(\cdot)$ also influences the posterior's concentration. Intuitively, a more “peaked” kernel assigns greater weight to small residuals (scaled by $h$), penalizing deviations from zero more sharply and thus yielding a more concentrated posterior for $\boldsymbol{\beta}$. We formalize this by examining the Hessian of the negative log-likelihood component of the posterior.

Let $U_L(\boldsymbol{\beta}; \theta, h, K) = \theta \sum_{i=1}^n L_h(e_i(\boldsymbol{\beta}); \tau)$ be the primary component of the negative log-posterior (or potential energy function) that depends on $\boldsymbol{\beta}$ through the sum of smoothed losses, with its argument $K$ signifying dependence on the kernel function $K(\cdot)$. The full negative log-posterior is $U(\boldsymbol{\beta}) = U_L(\boldsymbol{\beta}; \theta, h, K) - \log \pi(\boldsymbol{\beta}) + C_{\theta,h}$, where $C_{\theta,h}$ collects terms not dependent on $\boldsymbol{\beta}$.

We focus on the Hessian of the likelihood component $U_L$, denoted as $\boldsymbol{H}_L(\boldsymbol{\beta}; K)$. Utilizing the rigorous derivation based on the Dirac delta function established in Section (ref) (Eq. (ref) and Eq. (ref)), we directly obtain:

A larger Hessian (in the positive definite sense) at the posterior mode $\check{\boldsymbol{\beta}}(\tau)$ suggests a more sharply peaked posterior and, via Laplace approximation, a smaller posterior covariance.

As discussed in Section (ref), the true errors $\varepsilon_{0i} = y_i - \boldsymbol{x}_i^\top\boldsymbol{\beta}_0(\tau)$ are assumed to be independent and identically distributed following a common density $f_{\varepsilon_0}(\cdot)$, where $\boldsymbol{\beta}_0(\tau)$ represents the true value of $\boldsymbol{\beta}(\tau)$. Define $s_K(h) \coloneqq \mathbb{E}_{\varepsilon_0 \sim f_{\varepsilon_0}}\left[K\left(\frac{\varepsilon_0}{h}\right)\right]$. Let $\mathcal{H}_L(\boldsymbol{\beta}_0(\tau); K) \coloneqq \mathbb{E}_{\varepsilon_{0i}, \boldsymbol{x}_i} [\boldsymbol{H}_L(\boldsymbol{\beta}_0(\tau); K)]$ denote the expected Hessian of the negative log-likelihood component, evaluated at the true parameter $\boldsymbol{\beta}_0(\tau)$. Note that under the i.i.d.\ assumption, this relates to the population Hessian of the loss defined in Eq. (ref) by a scaling factor: $\mathcal{H}_L(\boldsymbol{\beta}_0(\tau); K) = n \theta \mathcal{H}(\boldsymbol{\beta}_0(\tau))$. This leads to the following theorem with the proof provided in \hyperref[pf:thm5]{Appendix A.4}.

Bayesian inference via MCMC

Bayesian inference for the BSQR parameters\,---\,regression coefficients $\boldsymbol{\beta}$ and scale $\theta$\,---\,is challenging, particularly when targeting uncertainty quantification. Standard BQR exploits the ALD’s scale-mixture-of-normals form Kozumi2011, Yu2001, enabling efficient Gibbs updates. However, the smoothed loss $L_h(\cdot; \tau)$ in our framework disrupts conditional conjugacy, rendering closed-form Gibbs sampling intractable. While Metropolis-Hastings within Gibbs is a viable alternative, it often suffers from slow mixing and requires careful tuning.

To overcome these limitations, we adopt the quasi-Bayesian framework of \shortciteA{Chernozhukov2003} for defining the posterior, and implement computation via Hamiltonian Monte Carlo (HMC) Duane1987, Neal2011 using the No-U-Turn sampler (NUTS) Hoffman2014 as implemented in Stan Carpenter2017. Unlike component-wise sampling schemes, NUTS updates the joint parameter space $(\boldsymbol{\beta}, \theta)$ simultaneously using Hamiltonian dynamics, leveraging gradient information to explore the high-dimensional posterior robustly.

Joint Hamiltonian Monte Carlo sampling

We target the joint posterior distribution $\pi(\boldsymbol{\beta}, \theta \mid \boldsymbol{y}, \boldsymbol{\mathcal{X}};\tau,h)$. The complete sampling procedure is summarized in Algorithm (ref) of \hyperref[sec:AppendixB]{Appendix B}. The HMC sampler augments the parameter space with auxiliary momentum variables.

To map the positively constrained scale parameter $\theta$ onto the unconstrained Euclidean space required by HMC, we utilize the logarithmic transformation (consistent with Stan's automatic parameter unconstraining mechanism). Let $\boldsymbol{q} = (\boldsymbol{\beta}^\top, \log\theta)^\top$ denote the extended state vector and $\boldsymbol{p}$ denote the corresponding momentum vector. The system's dynamics are governed by the Hamiltonian function: \[ H(\boldsymbol{q}, \boldsymbol{p}) = U(\boldsymbol{q} \mid \boldsymbol{y}, \boldsymbol{\mathcal{X}}) + T(\boldsymbol{p}), \] where $T(\boldsymbol{p}) = \frac{1}{2}\boldsymbol{p}^\top \boldsymbol{M}^{-1} \boldsymbol{p}$ is the kinetic energy with mass matrix $\boldsymbol{M}$ (typically diagonal and adapted during warmup), and $U(\boldsymbol{q} \mid \boldsymbol{y}, \boldsymbol{\mathcal{X}})$ is the potential energy. Since the sampling is performed on the unconstrained scale $\boldsymbol{q}$ (involving $\log \theta$), $U(\boldsymbol{q}\mid \boldsymbol{y}, \boldsymbol{\mathcal{X}})$ must include the Jacobian adjustment to ensure the correct marginal posterior density for $\theta$. Thus, the potential energy is defined as:

The evolution of the system over fictitious time $s$ follows Hamilton’s equations of motion Duane1987, Neal2011:

These equations preserve phase-space volume and are reversible. In practice, they are integrated using the leapfrog method (or its variant in NUTS). A single step of size $\epsilon$ updates the joint state $(\boldsymbol{q}, \boldsymbol{p})$ via:

Crucially, simulating these dynamics requires computing the gradients of the potential energy $\nabla_{\boldsymbol{q}} U(\boldsymbol{q} \mid \boldsymbol{y}, \boldsymbol{\mathcal{X}})$, which involve partial derivatives with respect to both $\boldsymbol{\beta}$ and $\theta$. For notational simplicity, we express these gradients in terms of the constrained parameters $(\boldsymbol{\beta}, \theta)$; however, in the actual HMC implementation, sampling is performed on the unconstrained state $\boldsymbol{q}=(\boldsymbol{\beta}^\top, \log\theta)^\top$, with the log-Jacobian determinant absorbed into $U(\boldsymbol{q} \mid \boldsymbol{y}, \boldsymbol{\mathcal{X}})$. For the regression coefficients $\boldsymbol{\beta}$, the gradient is derived analytically using the smoothed loss properties established in Section (ref): \[ \nabla_{\boldsymbol{\beta}} U(\boldsymbol{\beta}, \theta \mid \boldsymbol{y}, \boldsymbol{\mathcal{X}}) = \theta \sum_{i=1}^n \nabla_{\boldsymbol{\beta}} L_h(e_i(\boldsymbol{\beta}); \tau) - \nabla_{\boldsymbol{\beta}} \log \pi(\boldsymbol{\beta}), \] where $\nabla_{\boldsymbol{\beta}} L_h(e_i;\tau) = -\Psi_h(e_i;\tau)\,\boldsymbol{x}_i$. Using the closed-form expression derived in Eq. (ref), the scalar score component is simply:

where $K(\cdot)$ is the kernel density and $F_K$ its CDF. We provide the specific forms of $\Psi_h(\cdot;\tau)$ and $L_h(\cdot;\tau)$ for various kernels in Sections (ref)--(ref).

For the scale parameter $\theta$, the gradient involves the derivative of the log-normalizing constant $\nabla_\theta \log Z(\theta, \tau, h)$. While this term involves an intractable integral that typically necessitates gradient-free methods like Metropolis-Hastings, we leverage Stan's automatic differentiation capability. By implementing $Z(\theta, \tau, h)$ via numerical quadrature (e.g., integrate_1d), the system automatically computes exact gradients through the integration operator, enabling fully gradient-based NUTS sampling for $\theta$ as well.

Kernel-specific smoothed gradients

Efficient HMC sampling requires evaluating the smoothed loss $L_h(\cdot;\tau)$ and score $\Psi_h(\cdot;\tau)$ at each leapfrog step. Below, we present the closed-form expressions derived for common kernels; for brevity, detailed step-by-step derivations are provided in \hyperref[sec:AppendixC]{Appendix C}. These derivations confirm that bounded-support kernels (Uniform, Triangular, Epanechnikov) yield piecewise polynomial potentials. Since HMC dynamics rely primarily on gradient information (requiring only $C^1$ continuity), these kernels remain computationally efficient and robust despite minor discontinuities in higher-order derivatives.

Gaussian kernel

The standard Gaussian kernel is $K(v) = \phi(v) = (2\pi)^{-1/2} \exp(-v^2/2)$, with CDF $F_K(u) = \Phi(u)$. Applying Eq. (ref) yields:

The corresponding smoothed loss function $L_h(\cdot; \tau)$, obtained by integrating $\Psi_h(\cdot; \tau)$ and consistent with $\rho_\tau(\cdot)$ for small $h$, is given by (e.g., \shortciteNP{Horowitz1998, Koenker2005}):

Uniform kernel

The standard Uniform kernel is $K(v) = \frac{1}{2}$ for $v \in [-1, 1]$ and $0$ otherwise. Its CDF is $F_K(u) = 0$ for $u < -1$, $F_K(u) = \frac{u+1}{2}$ for $-1 \le u \le 1$, and $F_K(u) = 1$ for $u > 1$. From Eq. (ref), $\Psi_h(e; \tau)$ is:

The function $L_h(e; \tau)$ is derived by integrating $\Psi_h(e; \tau)$ and imposing continuity with $\rho_\tau(e)$ for $|e/h| \ge 1$. The resulting expression is:

Epanechnikov kernel

The standard Epanechnikov kernel is $K(v) = \frac{3}{4}(1-v^2)$ for $v \in [-1, 1]$ and $0$ otherwise. Its CDF, $F_K(u)$, is $0$ for $u < -1$, $\frac{3}{4}u - \frac{1}{4}u^3 + \frac{1}{2}$ for $-1 \le u \le 1$, and $1$ for $u > 1$. This yields $\Psi_h(\cdot; \tau)$ as:

Integration of $\Psi_h(e; \tau)$ and matching boundary conditions with $\rho_\tau(e)$ for $|e/h| \ge 1$ leads to $L_h(\cdot; \tau)$:

Triangular kernel

The standard Triangular kernel is $K(v) = 1-|v|$ for $v \in [-1, 1]$ and $0$ otherwise. The CDF, $F_K(u)$, is $0$ for $u < -1$; $\frac{1}{2}(1+u)^2$ for $-1 \le u < 0$; $1 - \frac{1}{2}(1-u)^2$ for $0 \le u \le 1$; and $1$ for $u > 1$. Consequently, $\Psi_h(\cdot; \tau)$ is:

The corresponding $L_h(e; \tau)$ is obtained by piecewise integration of $\Psi_h(e; \tau)$, ensuring continuity with $\rho_\tau(e)$ for $|e/h| \ge 1$ and at $e=0$:

Computational implementation details

The algorithm is implemented in Stan (code available in supplementary materials). To ensure numerical stability and computational efficiency in the NUTS, we employ two key strategies corresponding to specific blocks in our Stan code:

This fully gradient-based approach, combined with these numerical stabilizations, eliminates the need for manual tuning of proposal variances inherent in Metropolis-Hastings steps and ensures robust mixing.

Simulation

We evaluate the estimation accuracy, inferential validity, and scalability of the BSQR framework across three regimes: (1) general performance in standard settings (Section (ref)); (2) inferential validity under extreme sparsity (Section (ref)); and (3) high-dimensional scalability (Section (ref)). Unless otherwise noted in the specialized regimes of Sections (ref) and (ref), synthetic data are generated via the linear specification $y_i = \boldsymbol{x}_i^\top\boldsymbol{\beta} + u_i$, utilizing four error distributions to assess robustness: (1) $\mathcal{N}(0, 1)$; (2) $t(3)$; (3) a mixture $0.2\mathcal{N}(0, 3) + 0.8\mathcal{N}(0, 4)$; and (4) heteroscedastic normal errors $u_i \mid \boldsymbol{x}_i \sim \mathcal{N}(0, \sigma_i^2)$ with $\sigma_i = \exp(-0.25 + 0.5 x_{i1})$, which violates the linear conditional quantile assumption.

Benchmarking against standard quantile regression (StdQR) (quantreg) and BQR-ALD (brms), we implement BSQR in Stan Carpenter2017 using the NUTS Hoffman2014. We run two parallel chains of 4000 iterations (2000 warmup) to update the joint space $(\boldsymbol{\beta}, \theta)$, leveraging Stan's automatic differentiation with a hybrid strategy for the normalizing constant $Z(\theta, \tau, h)$: employing asymptotic approximation for sharply peaked integrands (large $\theta$) and robust numerical quadrature (integrate_1d) otherwise. The smoothing bandwidth $h$, balancing check loss fidelity and posterior stability, is selected via 5-fold cross-validation from a set of candidates scaled from Silverman's rule of thumb by factors $\{0.5, 0.75, 1.0, 1.5, 2.0\}$; fully Bayesian learning of $h$ is computationally prohibitive due to the intractability of $Z(\theta, \tau, h)$, though sensitivity is analyzed in Section (ref).

Priors are set to $\boldsymbol{\beta} \sim \mathcal{N}(\boldsymbol{0}, 1000\boldsymbol{I})$ and $\theta \sim \text{Gamma}(0.01, 0.01)$. While Corollary (ref) implies fixing $\theta$ for asymptotic coverage, we treat it as random to integrate out scale uncertainty, avoiding the noise of finite-sample Hessian/Jacobian estimation while empirically satisfying the generalized information equality. Point estimation accuracy is measured by MSE: $\|\hat{\boldsymbol{\beta}}(\tau) - \boldsymbol{\beta}_0(\tau)\|_2^2$; MAE: $\|\hat{\boldsymbol{\beta}}(\tau) - \boldsymbol{\beta}_0(\tau)\|_1 / d$; WMSE: $(\hat{\boldsymbol{\beta}}(\tau) - \boldsymbol{\beta}_0(\tau))^\top \boldsymbol{\Sigma}_X (\hat{\boldsymbol{\beta}}(\tau) - \boldsymbol{\beta}_0(\tau))$, with $\hat{\boldsymbol{\beta}}(\tau)$ the posterior mean. Predictive accuracy is the average check loss on the test set: $N_{\text{test}}^{-1} \sum_{i=1}^{N_{\text{test}}} \rho_{\tau}(y_i^{\text{test}} - \boldsymbol{x}_i^{\text{test}^\top} \hat{\boldsymbol{\beta}}(\tau))$. For Bayesian methods, we also report the empirical coverage and average width of 95% credible intervals, computation time (seconds), and MCMC diagnostics, including the maximum potential scale reduction factor ($\widehat{R}_{\max}$) and the minimum bulk effective sample size (ESS$_{\min}$) across all $\boldsymbol{\beta}$ coefficients.

General performance in standard settings

We first evaluate performance in standard sample sizes with $N_{\text{train}}=200$ and $N_{\text{test}}=1000$. Covariate vectors $\boldsymbol{x}_i \sim \mathcal{N}(\boldsymbol{0}, \boldsymbol{\Sigma}_X)$ have an autoregressive structure $(\boldsymbol{\Sigma}_X)_{jk} = \rho^{|j-k|}$ with $\rho=0.5$, as in \shortciteA{Fan2001}. We consider two specific designs: a sparse, high-dimensional setting ($d=20$) with $\boldsymbol{\beta}_0 = (3, 1.5, 0, 0, 2, 0, \dots, 0)^\top$, and a dense, lower-dimensional setting ($d=8$) with $\boldsymbol{\beta}_0 = (0.85, \dots, 0.85)^\top$. Simulations are run for $\tau \in \{0.25, 0.5, 0.75\}$ over $M=200$ independent replications.

Results (detailed in \hyperref[tab:sim_results_all_kernels]{Appendix D.1}) validate the proposed BSQR framework, showing consistent superiority over BQR-ALD across multiple evaluation dimensions. BSQR generally excels in estimation accuracy and prediction, achieving lower MSE, MAE, and WMSE for coefficient estimates in most scenarios across error distributions, though BQR-ALD is occasionally comparable or slightly better in specific cases. Critically, BSQR resolves BQR-ALD's predictive bias: out-of-sample check loss is reduced by 40--50% at $\tau=0.25$ and $\tau=0.75$, matching or even slightly surpassing the frequentist StdQR benchmark, and performances are comparable at $\tau=0.5$, confirming that smoothing restores the link between model parameters and the true conditional quantiles. For inference, BSQR delivers more reliable uncertainty quantification, with credible interval coverage often closer to 0.95; BQR-ALD's narrower intervals frequently under-cover, while Uniform and Triangular kernels balance coverage and width effectively. Computationally, the smoothed posterior yields 20--40% higher minimum bulk ESS (e.g., 2700--3200 vs. 2100--2400 for BQR-ALD), ensuring reliable estimates, and bounded-support kernels (Uniform, Triangular, Epanechnikov) maintain near-zero divergent transitions (consistently $<2$); Uniform (BSQR-U) is the fastest for large datasets\,---\,often surpassing BQR-ALD\,---\,while Triangular (BSQR-T) offers a favorable trade-off between speed, accuracy, and inference, Epanechnikov (BSQR-E) is slower with wider intervals, and Gaussian (BSQR-G) is computationally intensive and least stable. BSQR improves BQR by eliminating bias, improving computational efficiency, and favoring simple bounded-support Uniform and Triangular kernels for optimal trade-offs in speed, accuracy, and inference, making them the preferred choice for applied research.

Inferential validity under extreme sparsity

To rigorously assess finite-sample robustness, we conduct a “stress test” characterized by extreme data sparsity ($n=30, p=6$) and heteroscedasticity. This regime severely challenges the asymptotic justifications of semi-parametric estimators. Diverging from the general simulation settings in Section (ref), we adopt a specialized data generating process to isolate the effects of sparsity. Specifically, covariates are generated from independent standard normal distributions, and errors follow a linear heteroscedastic specification $u_i \sim \mathcal{N}(0, (1 + 0.5|x_{i1}|)^2)$, with the true coefficient vector set to $\boldsymbol{\beta}_0 = (1, 1, 0, 0, 0, 0)^\top$.

Furthermore, regarding computational implementation, given the extremely small sample size, cross-validation is prone to instability. Thus, distinct from the grid-search strategy employed in Section (ref), here we employ a fixed bandwidth determined by Silverman's rule of thumb for all BSQR variants. This constraint ensures that performance differences are driven solely by the kernel's geometric properties rather than bandwidth variability. We benchmark BSQR against the PETEL method of \shortciteA{Tang2022}, focusing on the critical trade-off between estimation precision and inferential validity at extreme quantiles ($\tau \in \{0.05, 0.95\}$) versus the median. The detailed numerical results are provided in Table (ref) in \hyperref[sec:AppendixD]{Appendix D.2}.

Table (ref) reveals a stark divergence in inferential reliability. While all methods successfully converged numerically in all replications (0% failure rate), the seemingly lower MSE of PETEL appears to be driven by an artificial reduction in variance at the cost of significant bias and under-coverage. Specifically, PETEL's squared bias is 3--8 times larger than that of BSQR with compact kernels (e.g., 91.4 vs. 27.0 at $\tau=0.05$). Consequently, its frequentist coverage drops to 83--84%, significantly below the nominal 95% level. This pattern suggests that the empirical likelihood surface becomes overly constrained in data-sparse regions, yielding estimates that are precise but inaccurate (“false precision”).

Conversely, BSQR offers two distinct advantages depending on the kernel choice. First, compact kernels (Uniform, Triangular, Epanechnikov) act as effective regularizers, reducing the squared bias by nearly an order of magnitude compared to PETEL (e.g., 6.26 vs. 49.7 at the median). Second, for valid uncertainty quantification, the Gaussian kernel functions as a robust inferential safety net. The infinite support of the Gaussian kernel provides non-zero gradient signals that guide the sampler even where empirical likelihoods face “empty set” constraints Koenker2005. Despite incurring a higher MSE and Interval Score (IS) due to the larger bandwidth required to bridge sparse data points, it is the only method that successfully maintains valid coverage ($>97\%$) at the tails. This reflects a theoretically justified trade-off: in regimes of critical information scarcity, ensuring conservative validity via smoothing is prioritized over the risk of over-confident, biased inference.

Scalability in high dimensions

To probe the computational boundaries of the competing methods, we simulate a high-dimensional regime ($n=300, p=50$). This setting exacerbates the “curse of dimensionality,” creating a parameter space where likelihood-free or moment-based methods often struggle due to complex posterior geometry. To explicitly isolate the computational cost of high dimensionality from statistical complications (such as multicollinearity or hyperparameter tuning), we adopt a simplified simulation protocols distinct from Section (ref). Specifically, regressors are generated from independent standard normal distributions, errors follow $\mathcal{N}(0,1)$, and the true coefficient vector is sparse with $\boldsymbol{\beta}_0 = (1, 1, 1, 1, 1, 0, \dots, 0)^\top$. Furthermore, we employ a fixed bandwidth derived from Silverman's rule for BSQR to strictly measure the sampling efficiency (ESS/sec) of the HMC kernel itself, excluding the overhead of cross-validation. We contrast the sampling efficiency of BSQR (utilizing the Gaussian kernel and HMC) with PETEL (utilizing random walk Metropolis), summarizing the results in Table (ref).

Table (ref) documents a catastrophic breakdown in the efficacy of the projected empirical likelihood approach. The PETEL algorithm exhibits a phenomenon of complete mixing collapse: despite successful numerical execution, the ESS is indistinguishable from zero, and the convergence diagnostics ($\widehat{R}_{\max}$) are undefined due to vanishing chain variance. This stagnation confirms that the Random Walk Metropolis sampler is ill-suited for the complex, constrained geometry of the empirical likelihood surface in 50 dimensions. Lacking gradient guidance, the blind proposal mechanism yields a near-zero acceptance probability, leaving the chain trapped at its starting value. Crucially, PETEL incurs a higher computational cost (10.20s vs. 6.15s) merely to reject every proposal. Consequently, the seemingly comparable MSE of PETEL ($6.44\times 10^{-3}$) is solely an artifact of the informative initialization provided by the frequentist estimator, rather than evidence of successful Bayesian learning.

In sharp contrast, the BSQR framework demonstrates HMC dominance in high-dimensional exploration. By smoothing the check loss, BSQR provides the necessary gradient information that allows the Hamiltonian Monte Carlo sampler to traverse the parameter space efficiently, overcoming the geometric bottlenecks identified by \shortciteA{Betancourt2017}. The sampler achieves near-perfect mixing ($\text{ESS} \approx 795$ out of 1000 iterations) and rapid convergence ($\widehat{R}_{\max} = 1.01$). This results in an effectively infinite efficiency gap: while BSQR generates approximately 129 independent posterior samples per second, PETEL generates none. These findings empirically establish that the smoothed likelihood formulation is not merely a theoretical convenience but a computational prerequisite for scalable Bayesian inference in high-dimensional quantile regression.

Empirical analysis: Asymmetric systemic risk exposure in the post-COVID era

We apply the BSQR framework to assess the asymmetric dynamics of systemic risk for a globally systemically important financial institution (G-SIFI). Quantile-based dependence measures, such as the quantilogram Linton2007, have highlighted the critical importance of tail dependency in financial econometrics. Building on this perspective, our analysis models the asymmetric link between daily stock returns of JPMorgan Chase & Co. (JPM), the largest U.S. bank by assets, and the S&P 500 Index (SP500),\footnote{Data for JPM and the S&P 500 are from \url{https://finance.yahoo.com/quote/JPM/} and \url{https://finance.yahoo.com/quote/^GSPC/}, respectively.} using daily log-returns from January 1, 2017, to January 1, 2025. Prices are converted to continuously compounded returns ($r_t = \ln(P_t/P_{t-1})$), yielding a rich time series whose summary statistics appear in Table (ref). The pronounced leptokurtosis in both series strongly motivates a quantile-based approach.

We specify a dynamic capital asset pricing model for quantile level $\tau$. Let $r_{jpm, t}$ and $r_{m, t}$ denote the log-returns of JPM and the SP500 at time $t$, respectively. The conditional quantile function is specified as:

where $\beta(\tau)$ measures systemic risk exposure. To capture asymmetry, we examine the downside beta ($\beta(0.05)$), quantifying exposure during severe downturns, and the upside beta ($\beta(0.95)$), measuring participation in market rallies. Using a rolling-window estimation with a one-year window (252 trading days) advanced monthly (21 days), we compare BSQR (Uniform and Triangular kernels) with the BQR-ALD benchmark across inferential stability, predictive accuracy, and sampler efficiency.

Inferential stability and economic insights from asymmetric betas

Methodologically, Figure (ref) demonstrates that BSQR provides substantially more stable and interpretable parameter estimates. In both panels, the BQR-ALD benchmark produces visibly wider and more erratic credible intervals (light gray shaded area) plagued by high-frequency noise that can lead to over-interpretation of statistically insignificant fluctuations. In contrast, BSQR’s inherent kernel smoothing attenuates this noise, yielding tighter credible intervals (orange and blue shaded areas) and smoother posterior mean paths. This enables a more reliable quantification of uncertainty, clearly distinguishing genuine shifts in the risk profile from mere sampling variation.

Economically, BSQR’s inferential clarity reveals a multi-phase and profoundly asymmetric response of JPM’s systemic risk to the COVID-19 shock. After an initial, transient symmetric spike in both betas during the early-2020 market panic, a more dominant and divergent dynamic quickly emerged. The downside beta ($\beta(0.05)$), previously stable around a pre-crisis baseline of roughly 0.97, decoupled from market downturns, plummeting to a trough of 0.53 by mid-2021 in a sustained “defensive decoupling” effect before settling into a new, lower regime averaging 0.92. This distinct regime shift highlights the model's capacity to capture rapid structural changes in systemic risk exposure, consistent with the structural break phenomena analyzed in quantile change-point literature Lee2018. In stark contrast, the upside beta ($\beta(0.95)$), originating from a similar 0.96 baseline, surged to a peak of 1.38 during the stimulus-fueled recovery, followed by a structural re-pricing that established a new, markedly lower equilibrium averaging 0.75. The superior stability of BSQR is what makes this nuanced dynamic sequence\,---\,a transient symmetric spike followed by a powerful asymmetric divergence\,---\,unambiguously clear.

Predictive accuracy and sampler efficiency

Beyond inferential quality, we assess BSQR's out-of-sample predictive accuracy and computational efficiency. In terms of predictive accuracy, measured by the average one-day-ahead check loss, Table (ref) shows that BSQR is highly competitive. For the downside quantile level ($\tau=0.05$), the check loss of BSQR-Uniform (0.002998) and BSQR-Triangular (0.003031) is economically indistinguishable from that of BQR-ALD (0.003027). Notably, for the upside quantile level ($\tau=0.95$), both BSQR methods (0.001338 and 0.001349) outperform the BQR-ALD benchmark (0.001363). This is a powerful result: BSQR’s substantial gains in inferential stability are achieved without sacrificing\,---\,and in the upside case, even enhancing\,---\,predictive accuracy.

Computationally, BSQR demonstrates superior sampler efficiency, with the Uniform kernel emerging as the dominant specification (Table (ref)). The average minimum bulk effective sample size (Avg. Min. ESS) for BSQR-Uniform (2554.72) is approximately 83% higher than that of BQR-ALD (1399.66). Crucially, this improvement is achieved with a lower computational time (0.69s vs. 0.77s), indicating that the combination of the Uniform kernel and HMC is strictly more efficient than the Gibbs sampling benchmark. While the BSQR-Triangular method incurs a higher computational cost (4.23s) due to the complexity of its smooth approximation, it similarly delivers a 79% increase in ESS, offering a robust alternative. Overall, BSQR-Uniform provides the optimal balance, delivering higher predictive accuracy and nearly double the sampling efficiency in less time.

\FloatBarrier

Robustness to bandwidth selection

To examine the influence of smoothing bandwidth (Section (ref)), we conducted a sensitivity analysis of the downside beta ($\beta(0.05)$) for a representative 252-day rolling window ending July 7, 2020, during significant market turmoil. Using the BSQR-Uniform method, we compared three bandwidths: the cross-validated choice ($h_{CV} \approx 1.08$), a smaller undersmoothed value ($h \approx 0.54$), and a larger oversmoothed one ($h \approx 2.16$). As shown in Figure (ref), the posterior location of the systemic risk parameter is sensitive to $h$, with the mean shifting from about 1.36 to 1.42 as smoothing increases, underscoring the need for a principled, data-driven bandwidth choice. Yet, despite this point-estimate sensitivity, the economic conclusion holds: all posterior distributions concentrate on values well above 1, confirming high systemic risk in the crisis period regardless of bandwidth. This analysis both quantifies model sensitivity and highlights the value of cross-validation as an objective, replicable criterion for balancing fit and smoothing.

Our empirical application shows BSQR to be a substantial methodological advance: it delivers more reliable and stable parameter estimates while maintaining competitive predictive performance. This enhanced inferential fidelity supports a more credible, nuanced understanding of economic phenomena, offering a powerful tool for researchers and practitioners prioritizing robustness.

\FloatBarrier

Discussion

This paper introduces Bayesian smoothed quantile regression (BSQR), a principled framework that reconciles the long-standing tension between computational efficiency and inferential validity in Bayesian quantile analysis. By constructing a fully differentiable pseudo-likelihood via kernel smoothing\,---\,a classical device for managing the bias--variance trade-off Gozalo2000\,---\,BSQR delivers decision-theoretically aligned inference for conditional quantiles, rectifies the finite-sample bias documented for BQR-ALD Gneiting2011, and enables efficient gradient-based MCMC such as HMC. Our theoretical analysis further establishes posterior consistency, Bernstein--von Mises (BvM) limits, and global posterior propriety, and it derives a generalized calibration condition for the scale parameter that guarantees valid frequentist coverage.

A distinct methodological advantage of BSQR over empirical-likelihood-based frameworks Tang2022 lies in its modularity and tractability. Whereas empirical likelihood methods often face severe geometric and computational hurdles in the presence of latent variables, many moment conditions, or high-dimensional constraints, BSQR provides an explicit, globally supported pseudo-likelihood. This “plug-and-play” structure facilitates seamless embedding into hierarchical and latent-variable architectures, while preserving a clear decision-theoretic interpretation. In particular, the differentiable formulation is well suited for extensions to endogeneity via control-function approaches Lee2007 and to semi- and nonparametric specifications (e.g., splines or Gaussian processes), directly addressing the limitations of purely linear quantile models.

Future research avenues are rich. Beyond a fully Bayesian treatment of bandwidth selection Silverman1986, the BSQR framework naturally extends to multivariate quantiles Hallin2010 and to dynamic models with time-varying parameters Gerlach2011, where smoothing-based pseudo-likelihoods can again regularize non-smooth objectives without sacrificing interpretability. Moreover, the coherent posterior quantification delivered by BSQR is directly applicable to probabilistic assessments of stochastic dominance, a central theme in modern econometrics Linton2005, allowing researchers to test dominance orders without resorting to complex resampling schemes. In this sense, BSQR is not merely a computational refinement but a structured way of combining frequentist smoothing precision with Bayesian probabilistic richness, opening a flexible platform for future developments in distributional and risk-oriented econometrics.

\phantomsection

Declaration of competing interest

\bookmark[ level=section, dest=\@currentHref,

] {Declaration of competing interest} The authors declare that they have no conflict of interest.

\phantomsection

Author contributions

\bookmark[ level=section, dest=\@currentHref,

] {Author contributions}

Bingqi Liu conceived the study, developed the BSQR framework, derived its theoretical properties, wrote the software, and conducted all numerical experiments. He also wrote the initial draft of the manuscript. Kangqiang Li provided assistance with the proofs of asymptotic posterior consistency and kernel effects on posterior concentration. Tianxiao Pang supervised the research and provided critical feedback on the manuscript. All authors reviewed and approved the final manuscript.

\phantomsection

Data and code availability

\bookmark[ level=section, dest=\@currentHref,

] {Data and code availability} The source code and data required to replicate all numerical results in this paper are publicly available on GitHub at the following repository: \url{https://github.com/BeauquinLau/BSQR}.

\phantomsection