  In spite of the recent surge of interest in quantile regression, joint
estimation of linear quantile planes remains a great challenge in statistics
and econometrics. We propose a novel parametrization that characterizes any
collection of non-crossing quantile planes over arbitrarily shaped convex
predictor domains in any dimension by means of unconstrained scalar, vector and
function valued parameters. Statistical models based on this parametrization
inherit a fast computation of the likelihood function, enabling penalized
likelihood or Bayesian approaches to model fitting. We introduce a complete
Bayesian methodology by using Gaussian process prior distributions on the
function valued parameters and develop a robust and efficient Markov chain
Monte Carlo parameter estimation. The resulting method is shown to offer
posterior consistency under mild tail and regularity conditions. We present
several illustrative examples where the new method is compared against existing
approaches and is found to offer better accuracy, coverage and model fit.
