2009-2
Gibbs Sampling Methods for Bayesian Quantile Regression
Hideo Kozumi Genya Kobayashi

Gibbs Sampling Methods for Bayesian Quantile Regression
Hideo Kozumi∗,a , Genya KobayashiaaGraduate School of Business Administration, Kobe University, 2-1 Rokko, Kobe 657-8501, Japan

Abstract
This paper considers quantile regression models using an asymmetric Laplace distribution from a Bayesianpoint of view. We develop a simple and eﬃcient Gibbs sampling algorithm for ﬁtting the quantile regression model based on a location-scale mixture representation of the asymmetric Laplace distribution. It isshown that the resulting Gibbs sampler can be accomplished by sampling from either normal or generalized inverse Gaussian distribution. We also discuss some possible extensions of our approach, including
the incorporation of a scale parameter, the use of double exponential prior, and a Bayesian analysis ofTobit quantile regression. The proposed methods are illustrated by both simulated and real data.
Key words: Asymmetric Laplace distribution, Bayesian quantile regression, double exponential prior,generalized inverse Gaussian distribution, Gibbs sampler, Tobit quantile regression

1. Introduction
Since the seminal work of Koenker and Bassett (1978), quantile regression has received increasing
attention both from a theoretical and from an empirical viewpoint. It is a statistical procedure based onminimizing sums of asymmetrically weighted absolute residuals and can be used to explore the relationship between quantiles of the response distribution and available covariates. Since a set of quantiles often
provides more complete description of the response distribution than the mean, quantile regression oﬀers
a practically important alternative to classical mean regression. There exists a large literature on quantile
regression methods and we refer to Yu et al. (2003) and Koenker (2005) for an overview.
Let yi be a response variable and xi a k × 1 vector of covariates for the i-th observation. Let q p (xi )denote the p-th (0 < p < 1) quantile regression function of yi given xi . Suppose that the relationshipbetween q p (xi ) and xi can be modeled as q p (xi ) = x′i β p , where β p is a vector of unknown parameters ofinterest. Then we consider the quantile regression model given byyi = x′i β p + ϵi ,(i = 1, . . . , n),where ϵi is the error term whose distribution (with density, say, f p (·)) is restricted to have the p-th quantile∫0
equal to zero, that is, −∞ f p (ϵi )dϵi = p.

∗Corresponding author
Email addresses: kozumi@kobe-u.ac.jp (Hideo Kozumi), genyako@hotmail.co.jp (Genya Kobayashi)
January 30, 2009

The error density f p (·) is often left unspeciﬁed in the classical literature. Thus, quantile regression
estimation for β p proceeds by minimizingn∑()ρ p yi − x′i β p ,

(1)

i=1

where ρ p (·) is the check (or loss) function deﬁned byρ p (u) = u {p − I(u < 0)} ,(2)

and I(·) denotes the usual indicator function. Since, however, the check function is not diﬀerentiable
at zero, we cannot derive explicit solutions to the minimization problem. Therefore, linear programing
methods are commonly applied to obtain quantile regression estimates for β p (see, e.g., Koenker andPark (1996) and Portnoy and Koenker (1997)).
From a Bayesian point of view, Walker and Mallick (1999), Kottas and Gelfand (2001), and Hanson
and Johnson (2002) considered median regression, which is a special case of quantile regression with
p = 0.5, and discussed nonparametric modeling for the error distribution based on either Pólya tree
or Dirichlet process priors. Regarding general quantile regression, Yu and Moyeed (2001) proposed a
Bayesian modeling approach by noting that minimizing (1) is equivalent to maximizing a likelihood
function under the asymmetric Laplace error distribution (see also Tsionas (2003)).
As discussed in Yu and Moyeed (2001), the use of the asymmetric Laplace distribution for the error terms provides a natural way to deal with the Bayesian quantile regression problem. However, theresulting posterior density for β p is not analytically tractable due to the complexity of the likelihood
function. Therefore, Yu and Moyeed (2001) considered Markov chain Monte Carlo (MCMC) methods
for posterior inference. Speciﬁcally, they used a random walk Metropolis algorithm with a Gaussian proposal density centered at the current parameter value. Although the random walk sampler is a convenient
choice to generate candidate values, the corresponding acceptance probability depends on the value of p
through the likelihood function. As a result, tuning parameters of proposals such as a proposal step sizeneed to be adjusted so as to attain some appropriate acceptance rates for each value of p, and this limits
the applicability of the random walk sampler in practice.
This paper considers Bayesian quantile regression models using the asymmetric Laplace distribution
and proposes MCMC methods that are not only computationally eﬃcient but also easy to implement.
In particular, we develop the Gibbs sampling algorithms based on a location-scale mixture representation of the asymmetric Laplace distribution. It is shown that the mixture representation provides fully
tractable conditional posterior densities and considerably simpliﬁes the existing estimation proceduresfor quantile regression models. To our knowledge, no previous work has employed a location-scale mixture representation to analyze quantile regression models. Furthermore, we show that our approach canreadily incorporate a scale parameter and can be directly extended to Tobit quantile regression.
The rest of the paper is organized as follows. In section 2 we present a mixture representation of anasymmetric Laplace distribution and derive full conditional densities for parameters. Section 3 discusses2some possible extensions of our approach, including the incorporation of a scale parameter, the use of
double exponential prior for regression coeﬃcients, and a Bayesian analysis of Tobit quantile regression
models. Some numerical results for both simulated and real data are presented in Section 4, and someconclusions are given in the ﬁnal section.
2. Posterior inference
2.1. Mixture representation
Following Yu and Moyeed (2001), we consider the linear model given byyi = x′i β p + ϵi ,(i = 1, . . . , n),and assume that ϵi has the asymmetric Laplace distribution with density{}f p (ϵi ) = p(1 − p) exp −ρ p (ϵi ) ,(3)

where ρ p (·) is deﬁned in (2). The parameter p determines the skewness of distribution and the p-th
quantile of this distribution is zero. It is also known that the mean and variance of the asymmetric
Laplace distribution with density (3) are given, respectively, byE(ϵi ) =1 − 2pp(1 − p)and Var(ϵi ) =1 − 2p + 2p2.
p2 (1 − p)2

Some other properties of the asymmetric Laplace distribution can be found in Yu and Zhang (2005).
As shown in Kotz et al. (1998), the asymmetric Laplace distribution has various mixture representations. For example, if ξ and η are independent and identical standard exponential distributions,ξ/p − η/(1 − p) has the asymmetric Laplace distribution. To develop Gibbs sampling algorithms for thequantile regression model, we utilize a mixture representation based on exponential and normal distributions, which is given by the following proposition.
Proposition. Let z be an standard exponential variable and u a standard normal variable. If a randomvariable ϵ follows the asymmetric Laplace distribution with density (3), then we can represent ϵ as alocation-scale mixture of normals given by√ϵ = θz + τ zu,
where
θ=

1 − 2pp(1 − p)and τ2 =2.
p(1 − p)

Proof. See Appendix.
From this result, the response yi can be equivalently rewritten as√
yi = x′i β p + θzi + τ zi ui ,3

(4)

where zi ∼ E(1) and ui ∼ N(0, 1) are mutually independent, and E(ψ) denotes an exponential distribution
with mean ψ. As the conditional distribution of yi given zi is normal with mean x′i β p + θzi and varianceτ2 zi , the joint density of y = (y1 , . . . , yn )′ is given by

 n
 n


∑ (yi − x′ β p − θzi )2 ∏ − 1 




i
2
f (y|β p , z) ∝ zi  exp 
,
−



2


2τ zi
i=1

(5)

i=1

where z = (z1 , . . . , zn )′ .
2.2. Gibbs sampler
To proceed a Bayesian analysis, we assume the priorβ p ∼ N(β p0 , B p0 ),(6)

where β p0 and B p0 are the prior mean and covariance of β p , respectively. As proved in Yu and Moyeed (2001), all posterior moments of β p exist under the normal prior (6). Using data augmentation, a
Gibbs sampling algorithm for the quantile regression model is constructed by sampling β p and z fromtheir full conditional distributions. Since (4) is a normal linear regression model conditionally on zi , it is
not diﬃcult to derive the full conditional density of β p given by
β p |y, z ∼ N(β̂ p , B̂ p ),
where
B̂−1
p =

n
∑
xi x′
i

τ2 zi
i=1

+ B−1
p0 ,

(7)

 n





∑ xi (yi − θzi )−1β̂ p = B̂ p 
+
B
β
.

p0
p0


2


τ zi
i=1

From (5) together with a standard exponential density, the full conditional distribution of zi is proportional to
z−1/2
exp
i

{

)}
1 ( 2 −12
− δ̂i zi + γ̂i zi ,2(8)where δ̂2i = (yi − x′i β)2 /τ2 and γ̂i2 = 2 + θ2 /τ2 . Since (8) is the kernel of a generalized inverse Gaussian
distribution, we have)1
zi |y, β p ∼ GIG , δ̂i , γ̂i ,
2
(

(9)

where the probability density function of GIG(ν, a, b) is given by{}(b/a)ν ν−1
1 2 −12
f (x|ν, a, b) =x exp − (a x + b x) , x > 0, −∞ < ν < ∞, a, b ≥ 0,
2Kν (ab)2
and Kν (·) is a modiﬁed Bessel function of the third kind (see Barndorﬀ-Nielsen and Shephard (2001)).
There exist eﬃcient algorithms to simulate from a generalized inverse Gaussian distribution (see, e.g.,Dagpunar (1989) and Hörmann et al. (2004)), so that our Gibbs sampler deﬁned in (7) and (9) can beeasily applied to quantile regression estimation.
We note that Tsionas (2003) also developed a Gibbs sampling algorithm for the quantile regression
model by employing a diﬀerent representation of the asymmetric Laplace distribution. Although it does4

not require Metropolis–Hastings steps, the Gibbs sampler proposed by Tsionas (2003) needs to update
each element of β p separately. Thus his sampling algorithm may produce highly correlated draws andbecome less eﬃcient than our algorithm (see Liu et al. (1994)). Furthermore, Tsionas (2003) mentionedthat his Gibbs sampler is complicated and can be somewhat slow when the number of observations islarge.
3. Some extensions
3.1. Inference with scale parameterIn the previous section, we have considered the quantile regression model without taking into accounta scale parameter. If one may be interested to introduce a scale parameter σ > 0 into the model, ourapproach can incorporate it by rewriting (4) as√
yi = x′i β p + σθzi + στ zi ui .

(10)

However this expression is not convenient to develop Gibbs samplers as the scale parameter appears in
the conditional mean of yi . Therefore, we reparameterize (10) as√
yi = x′i β p + θvi + τ σvi ui ,where vi = σzi . To complete the model speciﬁcation, we assume that β p ∼ N(β p0 , B p0 ) and σ ∼IG(n0 /2, s0 /2), where IG(a, b) denotes an inverse Gamma distribution with parameters a and b.
We now need to sample β p , v = (v1 , . . . , vn )′ and σ from their conditional distributions. The usual
Bayesian calculations show that the full conditional density of β p is given by
β p |y, v, σ ∼ N(β̃ p , B̃ p ),where B̃−1
p =

(11)

{
}
xi x′i−1 and β̃ = B̂ ∑n xi (yi −θvi ) + B−1 β+
B. Similarly to the previous section,
p
2
2
p
p0
i=1 τ σvi
i=1 τ σvi
p0
p0

∑n

we can easily obtain that()
1vi |y, β p , σ ∼ GIG , δ̃i , γ̃i ,2(12)where δ̃2i = (yi − x′i β)2 /τ2 σ, and γ̃i2 = 2/σ + θ2 /τ2 σ. By noting that vi ∼ E(σ), the full conditional
density of σ is proportional to( ) n20 + 32 n+1nn∑∑ 1 (yi − x′i β p − θvi )2 


1
s


0
exp − 
v
+
+
,

i

2

σ
σ
2τ2 vii=1
i=1
so that we have( ñ s̃ )(13)σ|y, β p , v ∼ IG , ,2 2∑∑
where ñ = n0 + 3n and s̃ = s0 + 2 ni=1 vi + ni=1 (yi − x′i β p − θvi )2 /τ2 vi . Consequently, the introductionof scale parameter does not cause any diﬃculties in our Gibbs sampling algorithm.

53.2. Double exponential priorInstead of a normal prior, we consider an important alternative, that is, a double exponential prior forβ p . The density of double exponential prior is given byπ(β p ) =k∏
j=1π(β p j ) ∝

k
∏

(
)
exp −λ0 |β p j − β p j0 | ,j=1

where β p j is the jth element of β p , and β p j0 and λ0 are hyperparameters. Using this prior distribution, Yu and Stander (2005) showed that all posterior moments of β p exist for Tobit quantile regression
models. Recently, Park and Casella (2008) considered this prior in the context of Lasso estimation (Tibshirani (1996)) and discussed the choice of hyperparameters.
As shown in Park and Casella (2008), the double exponential density can be expressed as
 2 

∫ ∞
2
 λ0 ω j 1 (β p j − β p j0 ) − dω jexp −exp
π(β p j ) =
√



2ω j
2 
2πω j
0
where ω j has an exponential distribution with mean 2/λ2 . This suggests the following hierarchicalrepresentation of the prior:
β p |ω ∼ N(β p0 , Ω),ωi ∼ E(2/λ20 ),where ω = (ω1 , . . . , ωk )′ and Ω is a diagonal matrix with the jth element ω j . It follows from thisspeciﬁcation that the full conditional distributions of β p , σ and v are the same as those under the normal
prior with B p0 replaced with Ω. Also, the full conditional density of ω j is proportional to
[
}]
1{
− 12
2 −1
2
ω j exp − (β p j − β p j0 ) ω j + λ0 ω j ,2which implies that)1
ω j |y, β p , σ, v ∼ GIG , |β p j − β p j0 |, λ0 .
2(

3.3. Tobit quantile regression
Tobit quantile regression models have received much attention in the classical literature (see, e.g.,
Powell (1986), Hahn (1995), Buchinsky and Hahn (1998) and Bilias et al. (2000)). Yu and Stander (2005)
proposed a Bayesian framework for Tobit quantile regression based on the asymmetric Laplace distribution. Here we show that our methodology is directly extended to the analysis of Tobit quantile regression
models.
As in the standard Tobit model, we assume that the response yi is generated according toy∗i if y∗i > 0,yi = 
0 if y∗i ≤ 0,y∗i = x′i β p + ϵi ,6(14)

where y∗i is a latent variable. To develop a Tobit quantile regression model, we again assume that ϵi has
the asymmetric Laplace distribution and rewrite (14) as√y∗i = x′i β p + θvi + τ σvi ui ,(15)

where ui ∼ N(0, 1) and vi ∼ E(σ).
Since the model (15) has a regression form conditionally on vi , the method developed by Chib (1992)
can be applied to the sampling of y∗i , that is,y∗i |y, β p , v, σ ∼ yi I(yi > 0) + T N (−∞,0] (x′i β p + θvi , τ2 σvi )I(yi = 0),
where T N (a,b] (µ, σ2 ) denotes a normal distribution with mean µ and variance σ2 truncated on the interval(a, b]. Moreover, by assuming that β p ∼ N(β p0 , B p0 ) and σ ∼ IG(n0 /2, s0 /2), the full conditional
distributions of β p , v and σ can be easily obtained from (11)–(13) by replacing yi everywhere with y∗i .
Thus our approach oﬀers a simple sampling method for Tobit quantile regression models compared withthat of Yu and Stander (2005), which relies on Metropolis-Hastings algorithms and requires to control
tuning parameters for each value of p.
4. Numerical examples
We present three examples to illustrate the Gibbs sampling methods developed in Sections 2 and 3.
The ﬁrst example uses a simulated data set and the second is based on the real data analyzed by Wanget al. (1998). While the ﬁrst and second examples estimate the quantile regression models, the third
example considers the Tobit quantile regression model using the data from Mroz (1987).
4.1. Simulated data
A data set of n = 100 observations was generated from the modely = β1 + β2 x2 + β3 x3 + β4 x4 + β5 x5 + ϵ,ϵ ∼ N(0, 1),

where each covariate was simulated from a uniform distribution on (−2, 2) and all the β j were set to one.
The normal prior (6) was used for β p and the hyperparameters were chosen as β p0 = 0 and B p0 = 100I.
We ﬁtted the quantile regression model without a scale parameter for p = 0.05, 0.5 and 0.95.
To assess the sampling eﬃciency of our algorithm, we estimated ineﬃciency factors for β p by running the Gibbs sampler for 15000 iterations with an initial burn-in of 5000 iterations. The ineﬃciency
factor is deﬁned as a ratio of the numerical variance of the sample mean from the Markov chain to thevariance from independent draws (see Chib (2001)). The results are summarized in Table 1 together withthose obtained from the random walk sampler of Yu and Moyeed (2001). When applying the random
walk sampler, we adjusted proposal step sizes in order that the ineﬃciency factors achieve the lowest
values.
It is observed that the ineﬃciency factors for the Gibbs sampler are smaller than those for the random
walk sampler in all the cases. This indicates that our Gibbs sampler is more eﬃcient than the random
7

walk sampler. It is also found that both the Gibbs sampler and the random walk algorithm attain thesmallest ineﬃciency factors in the case of p = 0.5.
4.2. Patent data
As an real data example, we consider the data examined by Wang et al. (1998). The data set contains
information about the number of patent applications from 70 pharmaceutical and biomedical companies in 1976. A more detailed explanation of the data can be found in Hall et al. (1988). Following
Tsionas (2003), we analyze the relationship between patents and research and development (R&D) by
estimating the modellog(1 + N) = β1 + β2 log(RD) + β3 log(RD)2 + β4 log( RD )+ ϵ,
S ALE

where N is the number of patent applications, RD is R&D spending, and RD/S ALE is the ratio of
R&D to sales. In this study, we included a scale parameter σ and assumed that β p ∼ N(0, 100I) andσ ∼ IG(5/2, 0.1/2). Using the prior speciﬁcation, we ran the Gibbs sampler for 20000 iterations anddiscarded the ﬁrst 50000 iterations as a burn-in period.
Figure 1 shows the sample autocorrelation functions of β p for p = 0.5. It is observed that the
autocorrelations from the Gibbs sampler generally decline to zero by lag 20. Although diﬀerent priors
are employed in Tsionas (2003), a comparison with the results shown in Figure 3 of Tsionas (2003)
shows that the autocorrelations from our Gibbs sampler decay more quickly than those from the Gibbs
sampler of Tsionas (2003). This is because the Gibbs sampler of Tsionas (2003) updates each elementof β p separately while our algorithm draws all the elements jointly.
4.3. Labor supply data
To illustrate our method for Tobit quantile regression models, we consider the data set from Mroz (1987),
which is comprised of 753 married women between the ages of 30 and 60. The response variable is givenby the total number of hours the wife worked for a wage outside the home during year 1975 and measured in 100 hours. Of the 753 women in the sample, 325 of the women worked zero hours and thecorresponding responses are treated as left censored at zero. As the explanatory variables, we include a
constant term, income which is not due to the wife (nwi f einc), years of education (educ), years of work
experience (exper), wife’s age (age), the number of children under 6 years old (kidslt6), and the numberof children over 6 years old (kidsge6).
For the regression coeﬃcients β p , we considered both the normal and the double exponential prior,and speciﬁed the following hyperparameters: β p0 = 0, B p0 = 100I and λ0 = 0.14. Note that theprior variances are the same under both the priors. Furthermore, we incorporated a scale parameter σ
and assumed that σ ∼ IG(10/2, 0.02/2). All results are based on an MCMC sample of 10000 draws
obtained after a burn-in of 5000 iterations.
Table 2 shows the posterior estimates of the parameters. From the table, we see that both the priors
yield very similar posterior estimates and notice that posterior means of some parameters (e.g., the coeﬃcients on educ and kidslt6) are diﬀerent across the diﬀerent values of p. To study more closely the8relationships between the quantiles and the posterior estimates, Figure 2 plots the posterior estimates ofβ p against various values of p under the normal prior. We observe from the ﬁgure that nwi f einc, kidslt6
and age are negatively associated with the wife’s work in the lower and middle quantiles, but nwi f eincand kidslt6 become less related in the higher quantiles. On the contrary, educ and exper have positiveeﬀects on the wife’s working time in the lower and middle quantiles, but they are also less related in thehigher quantiles. Finally, the results for expersq and kidsge6 show almost ﬂat paths compared with theother covariates, indicating that they have constant eﬀects on the wife’s propensity to work.
5. Conclusions
We have developed a Gibbs sampling method for quantile regression models based on the locationscale mixture representation of the asymmetric Laplace distribution. The proposed Gibbs sampling algorithm is easy to implement in practice since one cycle of the algorithm can be accomplished by thesimulation from a normal and a generalized inverse Gaussian distribution. We have also discussed someextensions of our approach, which are inference with scale parameters, the use of double exponential
priors, and Tobit quantile regression analysis. Finally, we have provided illustrations with simulated andreal data sets and shown the superiority of our Gibbs sampler to the existing sampling methods.
The mixture representation utilized in this paper allows us to express a quantile regression model as anormal regression model. Therefore, our approach can be further extended to more complicated models
such as nonlinear models, and this is left for future research.
Acknowledgments
This research is partly supported by Grants-in-Aid for Scientiﬁc Research from the Ministry of Education, Culture, Sports, Science and Technology of Japan. The computational results are obtained using
Ox version 5.10 (Doornik (2007)).
[2009.1.30. 906]9

Appendix
Suppose that a random variable ϵ has the probability density function given by}{f p (ϵ) = p(1 − p) exp −ρ p (ϵ) ,where ρ p (ϵ) = ϵ{p − I(ϵ < 0)}. Then the characteristic function of ϵ is derived asϕ(t) = E[eitϵ ]∫ ∞=eitϵ f p (ϵ) dϵ
−∞
∫ 0
∫ ∞
itϵ+(1−p)ϵ=p(1 − p)edϵ +p(1 − p)eitϵ−pϵ dϵ
−∞
0
}
{
1
1
+
= p(1 − p)it + (1 − p) p − it{}−1
1
1 − 2p
2
=
t −i
t+1
,
p(1 − p)p(1 − p)where i2 = −1.
Let z be an standard exponential variable and u a standard normal variable. The characteristic func√tion of ϵ ′ = θz + τ zu can be expressed as[√ ]φ(t) = E eit(θz+τ zu)
∫ ∞
[ √ ]
eitθz E eitτ zu e−z dz.
=
0

Since u follows a standard normal distribution, we have
[ √ ]
1 2 2
E eitτ zu = e− 2 t τ z .
Thus, the characteristic function of ϵ ′ is obtained as
∫ ∞
1 2 2
φ(t) =e−z(1+ 2 t τ −iθt) dz
(

0

)−1
1 22
= τ t − iθt + 1 .
2
1−2p2
These two characteristic functions are equivalent when θ = p(1−p)
and τ2 = p(1−p).

10References[1] Barndorﬀ-Nielsen, O.E. and Shephard, N. (2001). “Non–Gaussian Ornstein-Uhlenbeck-basedmodels and some of their uses in ﬁnancial economics,” Journal of the Royal Statistical Society
B63, 167–241.
[2] Bilias, Y., Chen, S., and Ying, Z. (2000). “Simple resampling methods for censored regression
quantiles,” Journal of Econometrics 68, 303–338.
[3] Buchinsky, M. and Hahn, J. (1998). “An alternative estimator for censored quantile regression,”
Econometrica 66, 653–671.
[4] Chib, S. (1992). “Bayes inference in the Tobit censored regression model,” Journal of Econometrics 51, 79–99.
[5] Chib, S. (2001). “Markov chain Monte Carlo methods: computation and inference,” in J.J. Heckman and E. Leamer (eds.), Handbook of Econometrics 5, 3569–3649. Amsterdam: North Holland.
[6] Dagnapur, J. (1989). “An easily implemented generalized inverse Gaussian generator,” Communications in Statistics–Simulation 18, 703–710.
[7] Doornik, J.A. (2007). Object-Oriented Matrix Programming Using Ox (3rd ed.). London: Timberlake Consultants Press and Oxford.
[8] Hahn, J. (1995). “Bootstrapping quantile regression estimators,” Econometric Theory 11, 105–
121.
[9] Hall, B. H., Cummins, C., Laderman, E., and Mundy, J. (1988). ”The R&D master ﬁle documentation,” Technical Working Paper 72. National Bureau of Economic Research, Cambridge, MA.
[10] Hanson, T. and Johnson, W.O. (2002). “Modeling regression error with a mixture of Pólya trees,”
Journal of the American Statistical Association 97, 1020–1033.
[11] Hörmann, W., Leydold, J., and Derﬂinger, G. (2004). Automatic Nonuniform Random Variate
Generation. Berlin: Springer.
[12] Koenker, R. (2005). Quantile Regression. New York: Cambridge University Press.
[13] Koenker, R. and Bassett, G. (1978). “Regression quantiles,” Econometrica 46, 33–50.
[14] Koenker, R. and Park, B.J. (1996). “An interior point algorithm for nonlinear quantile regression,”
Journal of Econometrics 71, 265–283.
[15] Kottas, A. and Gelfand, A.E. (2001). “Bayesian semiparametric median regression modeling,”
Journal of American Statistical Association 96, 1458–1468.

11

[16] Kotz, S., Kozubowski, T.J., and Podgórski, K. (2001). The Laplace Distribution and Generalizations: A Revisit with Applications to Communications, Economics, Boston: Engineering, andFinance. Birkhäuser.
[17] Liu, J.S., Wong, W.H., and Kong, A. (1994). “Covariance structure of the Gibbs sampler withapplications to the comparisons of estimators and augmentation schemes,” Biometrika 81, 27–40.
[18] Mroz, T. (1987). “The sensitivity of an empirical model of married women’s hours of work toeconomic and statistical assumptions,” Econometrica 55, 765–799.
[19] Park, T. and Casella, G. (2008). “The Bayesian Lasso,” Journal of the American Statistical Association 103, 681–686.
[20] Portnoy, S. and Koenker, R. (1997). “The Gaussian hare and the Laplacian tortoise: computabilityof squared-error versus absolute-error estimators (with discussion),” Statistical Science 12, 279–300.
[21] Powell, J. (1986). “Censored regression quantiles,” Journal of Econometrics 32, 143–155.
[22] Tibshirani, R. (1996). “Regression shrinkage and selection via the lasso,” Journal of the Royal
Statistical Society B58, 267–288.
[23] Tsioans, E. (2003). “Bayesian quantile inference,” Journal of Statistical Computation and Simulation 73, 659–674.
[24] Walker, S.G. and Mallick, B.K. (1999). “A Bayesian semiparametric accelerated failure timemodel,” Biometrics 55, 477–483.
[25] Wang, P., Cockburn, I., and Puterman, M. K. (1998). “Analysis of patent data – A mixed-Poissonregression-model approach,“ Journal of Business and Economic Statistics 16, 2741.
[26] Yu, K., Lu, Z., and Stander, J. (2003). “Quantile regression: applications and current researchareas,” The Statistician 52, 331–350.
[27] Yu, K. and Moyeed, R.A. (2001). “Bayesian quantile regression,” Statistics & Probability Letters
54, 437–447.
[28] Yu, K. and Stander, J. (2005). “Bayesian analysis of a Tobit quantile regression model,” Journalof Econometrics 137, 260–276.
[29] Yu, K. and Zhang, J. (2005). “A three-parameter asymmetric Laplace distribution and its extension,” Communications in Statistics – Theory and Methods 34, 1867–1879.

12

Table 1: Simulated data: Ineﬃciency factors are shown.

Gibbs sampler

Random walk sampler

p

0.05

0.5

0.95

0.05

0.5

0.95

β1

11.902

3.244

14.14159.48616.34845.394β2

12.1503.904

16.66356.312

14.83251.253β320.085

3.540

11.26271.70824.90420.252β413.7435.564

17.46181.10424.57859.490β510.572

3.615

11.48639.88021.22169.212

Table 2: Labor supply data: Posterior means and standard deviations are shown.

Normal priorExplanatory

p = 0.05

p = 0.5

p = 0.95

Variable

Mean

SD

Mean

SD

Mean

SD

constant-5.8333.083

11.8834.120

22.8033.558

nwi f einc

-0.132

0.043-0.0990.045-0.059

0.036

educ

0.858

0.139

0.859

0.209

0.331

0.178

exper

0.948

0.155

1.413

0.181

0.798

0.181

exper2-0.0180.004-0.0180.006-0.0120.006

age-0.3540.057-0.6070.072-0.2690.065

kidslt6-7.2371.336-9.6301.142-4.855

1.145

kidsge6-0.1310.278-0.437

0.394

0.223

0.400

Double exponential priorExplanatory

p = 0.05

p = 0.5

p = 0.95

Variable

Mean

SD

Mean

SD

Mean

SD

constant-5.2903.201

11.2744.551

23.8983.653

nwi f einc

-0.1290.041-0.1000.045-0.057

0.037

educ

0.839

0.150

0.883

0.216

0.307

0.177

exper

0.929

0.142

1.410

0.182

0.769

0.188

exper2-0.0170.004-0.0180.006-0.0110.006

age-0.3610.058-0.6000.074-0.2820.067

kidslt6-7.1791.243-9.611

1.160-4.820

1.236

kidsge6-0.1590.256-0.402

0.404

0.160

0.380

13

1.0

constant

1.0

0.5

0.5

0.0

0.0

−0.5

−0.5

0
log(RD)2
1.0

10

20

30

0
10
log(RD/Sale)
1.0

0.5

0.5

0.0

0.0

−0.5

−0.5

0

10

20

log(RD)

30

0

10

20

30

20

30

Figure 1: Patent data: Sample autocorrelation functions are shown for p = 0.5.

14

50

constant

0.1

nwifeinc

0.0
25
−0.1
−0.2

0
0.0
0.2
educ

0.4

0.6

0.8

1.0

0.0
0.2
exper

0.4

0.6

0.8

1.0

0.0
age

0.2

0.4

0.6

0.8

1.0

0.0
0.2
kidsge6

0.4

0.6

0.8

1.0

0.0

0.4

0.6

0.8

1.0

1.5
1.5
1.0

0.5

0.5
0.0
0.2
exper 2

0.4

0.6

0.8

1.0

−0.25

0.00

−0.50
−0.02
−0.75
0.0
0.2
kidslt6

0.4

0.6

0.8

1.0
1

0
−5

−1

−10
0.0

0.2

0.4

0.6

0.8

1.0

0.2

Figure 2: Labor supply data: Posterior means (solid line) and 95% credible intervals (dotted lines) are plotted for the normal
prior.

15