# Bayesian Estimation of Simultaneous Regression Quantiles Using Hamiltonian Monte Carlo
## 4 Numerical Results
### Univariate Linear Quantile  
Suppose the error terms \(\epsilon_{i=1,\dots,T}\stackrel{{ i.i.d}}{{\sim}}\text{ALD}(0,0.1,0.75)\). As the Bayesian approach requires a likelihood, we use the ALD to make the Bayesian inference, and our quantiles of interest are \(\tau=0.25\) and \(\tau=0.5\); hence, we use the relationship between the two ALDs mentioned in Section 2.  
We independently generate 300 observations, using R, issued from the model defined in Equation (5). Since we are using simulated data, the theoretical quantiles are \(\beta_{0,0.25}=0.56\), \(\beta_{0,0.5}=0.83\), and \(\beta_{1,0.25}=\beta_{1,0.5}=2\).  
To study the quality of the estimate of the quantiles, we calculate the root mean square error:  
\[RMSE=\sqrt{\sum_{i=1}^{n}(q_{\tau}-\hat{q}_{\tau})^{2}},\]  
where \(q_{\tau}\) is the theoretical quantile and \(\hat{q}_{\tau}\) is the posterior mean quantile regression.  
To assess the impact of data variance on estimation, we assumed that the error terms \(\epsilon_{i=1,\dots,T}\stackrel{{ i.i.d}}{{\sim}}\text{ALD}(0,0.3,0.75)\). The obtained results are presented in Table 1.  
We run our algorithm for 20,000 iterations, with half of them being burn-in. Gelman-Rubin diagnostic (R-hat) values are depicted in Figure 1, with \(\hat{R}\) approaching 1, indicating the convergence of multiple chains. The densities of the estimated \(\beta\)s are depicted in Figure 2, providing insight into the distribution patterns.  
Furthermore, Figure 3 illustrates that the autocorrelation of each chain converges to 0. While comparing the HMC method used in this paper with the Metropolis-Hastings within Gibbs algorithm used in [27], we find that using HMC exhibits lower autocorrelation at earlier lags (around 0 after 10 lags) compared to the Metropolis-Hastings within Gibbs algorithm, where autocorrelation remains around 0 only after 50 lags. This suggests that HMC achieves greater independence between samples sooner, indicating potentially more efficient sampling and better exploration of the distribution space compared to the algorithm used in [27].  
Additionally, the examination of the trace plots in Figure 4 reveals consistent and stable behavior in the MCMC chains, showcasing no discernible oscillation, drift, or other abnormalities. This pattern suggests a thorough exploration of the parameter space, indicating successful convergence to the target distribution.  
The small variance observed in Table 1 further supports the reliability of the results. Furthermore, the Monte Carlo standard error (MCSE) is 0 for all \(\beta\)s, indicating that the estimated values are likely close to the true values, as shown in Table 1. In addition, when the variance of the data is larger, the estimation worsens (higher RMSE and higher variance). This is normal, as a larger variance makes it harder to estimate the parameters.  
It should be noted that diagnostics can only reliably be used to determine a lack of convergence and not detect convergence [47]. Therefore, from the obtained results, the trace plots do not seem to indicate non-stationary behavior.  
Figure 1: \(\hat{R}\) for univariate model.  
Figure 3: Autocorrelation plots for various \(\beta\) values in the univariate model (\(\beta_{1}\) on the left, \(\beta_{2}\) in the middle, and \(\beta_{3}\) on the right), showing chain 1 (**top**), chain 2 (**middle**), and chain 3 (**bottom**).  
Figure 2: Histograms showcasing diverse \(\beta\) values for the univariate model.