Multivariate Quantile Regression: Methods & Applications

• Multivariate quantile regression is a diverse framework for modeling conditional distributions of vector-valued responses, addressing the lack of a natural total ordering.
• Key approaches include directional quantiles, spatial quantiles via optimal transport, and sign-concordance methods, each targeting different inferential goals.
• MQR is applied in fields like biostatistics, finance, and forecasting to improve tail risk assessment and reveal joint dependency structures.

Multivariate quantile regression (MQR) studies how conditional distributional features of a vector-valued response vary with covariates when the target is not scalar and therefore lacks a canonical total ordering. In the recent literature, MQR denotes a family of constructions rather than a single model: some approaches define directional quantile hyperplanes and regions, some represent vector or geometric quantiles through optimal transport or spatial quantiles, some jointly estimate marginal quantiles under a multivariate asymmetric Laplace likelihood, some model joint sign patterns of marginal quantile residuals, and some select a representative point on a joint conditional quantile contour through sequential conditioning (Santos et al., 2019, Bhattacharya et al., 2020, Petrella et al., 2018, Columbu et al., 2021, Galvao et al., 21 Aug 2025). This plurality reflects a central technical fact: a unique definition of multivariate quantiles is lacking, so the field is organized around alternative quantile objects, estimands, and inferential targets rather than a single canonical extension of scalar quantile regression (Columbu et al., 2021).

1. Scalar quantile regression and the multivariate extension problem

For a scalar response, quantile regression at level $\tau \in (0,1)$ is typically written as

$Q_{Y_j}(\tau \mid X_i)=X_i^\top \beta_j(\tau),$

with $\hat\beta_j(\tau)$ obtained by minimizing the check loss

$$\sum_{i=1}^n \rho_\tau(Y_{ij}-X_i^\top\beta), \qquad \rho_\tau(u)=u\bigl(\tau-I[u<0]\bigr).$$

This formulation is standard in several multivariate extensions, where it is applied either margin by margin or as a component of a larger joint model (Columbu et al., 2021).

The multivariate case is difficult because the response $Y$ takes values in $\mathbb R^k$ or, in more recent work, on a manifold, and there is no natural scalar ordering analogous to $\le$ on $\mathbb R$. Several papers state this obstacle explicitly: multivariate quantiles are not straightforward to define, and classical Koenker–Bassett quantile regression extends poorly to multivariate $Y$ because of the lack of total ordering (Santos et al., 2019, Bhattacharya et al., 2020). When the target is non-Euclidean, such as data on spheres or tori, the issue is compounded by the geometry of the response space (Pegoraro et al., 2023).

A useful distinction runs through the literature. Some methods aim to define a genuinely joint multivariate quantile object—such as a region, surface, transport map, or contour—while others retain marginal scalar quantiles and then model the dependence structure induced by those marginal fits. This distinction is especially clear in the contrast between directional or transport-based formulations and the sign-concordance approach, which deliberately avoids introducing an explicit multivariate quantile function (Columbu et al., 2021).

2. Principal formulations of multivariate quantiles

Formulation	Quantile object	Typical estimation strategy
Directional quantiles	Hyperplanes, halfspaces, nested regions	Direction-wise check-loss regression with noncrossing adjustment
Quantile surfaces	Direction-indexed radial surface around a center	Directional neural quantile regression
Geometric/vector quantiles	Spatial quantile or optimal-transport map	Bayesian mixture models or OT-based optimization
Sign-concordance	Joint distribution of residual sign patterns	Separate marginal QR plus multinomial regression
MAL-based joint regression	Componentwise marginal quantiles under a joint likelihood	EM or Bayesian MCMC
Sequential conditional-CDF MQR	Ordered point on a joint $\tau$-contour	Sequential univariate QR with generated regressors
Multivariate M-quantiles	M-quantile profiles with shared random effects	Finite-mixture likelihood and EM

In the directional family, a unit direction $Q_{Y_j}(\tau \mid X_i)=X_i^\top \beta_j(\tau),$0 defines a projection of the response and an associated quantile hyperplane. In the structured-additive Bayesian framework, one writes $Q_{Y_j}(\tau \mid X_i)=X_i^\top \beta_j(\tau),$1 and $Q_{Y_j}(\tau \mid X_i)=X_i^\top \beta_j(\tau),$2, estimates directional hyperplanes by minimizing the check loss, and forms the $Q_{Y_j}(\tau \mid X_i)=X_i^\top \beta_j(\tau),$3-quantile region as

$Q_{Y_j}(\tau \mid X_i)=X_i^\top \beta_j(\tau),$4

Under mild regularity, Hallin et al. (2010) show that this region coincides with the Tukey-depth region, so directional quantiles and depth-based geometry are linked directly (Santos et al., 2019).

Quantile surfaces use a related but distinct geometric parameterization. Fixing a central tendency point $Q_{Y_j}(\tau \mid X_i)=X_i^\top \beta_j(\tau),$5, one projects $Q_{Y_j}(\tau \mid X_i)=X_i^\top \beta_j(\tau),$6 onto each direction $Q_{Y_j}(\tau \mid X_i)=X_i^\top \beta_j(\tau),$7, estimates the scalar directional quantile $Q_{Y_j}(\tau \mid X_i)=X_i^\top \beta_j(\tau),$8, and interprets $Q_{Y_j}(\tau \mid X_i)=X_i^\top \beta_j(\tau),$9 as a ray length. The set of endpoints

$\hat\beta_j(\tau)$0

defines the $\hat\beta_j(\tau)$1-quantile surface, which is also characterized as the boundary of the smallest star-shaped set with probability mass $\hat\beta_j(\tau)$2 (Bieshaar et al., 2020).

A different tradition defines multivariate quantiles through spatial or geometric criteria. In the dependent Dirichlet process model, the geometric quantile at level $\hat\beta_j(\tau)$3 is

$\hat\beta_j(\tau)$4

and is recovered from a flexible Bayesian model for the full conditional density $\hat\beta_j(\tau)$5 (Bhattacharya et al., 2020). On manifolds, vector quantile regression is formulated through Riemannian optimal transport with cost $\hat\beta_j(\tau)$6 and conditional vector quantile map

$\hat\beta_j(\tau)$7

where $\hat\beta_j(\tau)$8 is a $\hat\beta_j(\tau)$9-concave potential (Pegoraro et al., 2023).

Other formulations avoid a joint geometric quantile altogether. The sign-concordance method first fits separate univariate quantile regressions, then models the joint distribution of the residual signs

$$\sum_{i=1}^n \rho_\tau(Y_{ij}-X_i^\top\beta), \qquad \rho_\tau(u)=u\bigl(\tau-I[u<0]\bigr).$$0

through multinomial logistic regression. Each sign pattern corresponds to an orthant in the joint outcome space relative to the marginal quantile hyperplanes, so the method targets tail concordance and discordance rather than an explicit vector-valued quantile function (Columbu et al., 2021).

A further formulation defines a multivariate quantile directly through the joint conditional CDF. For $$\sum_{i=1}^n \rho_\tau(Y_{ij}-X_i^\top\beta), \qquad \rho_\tau(u)=u\bigl(\tau-I[u<0]\bigr).$$1, a level-$$\sum_{i=1}^n \rho_\tau(Y_{ij}-X_i^\top\beta), \qquad \rho_\tau(u)=u\bigl(\tau-I[u<0]\bigr).$$2 quantile graph satisfies

$$\sum_{i=1}^n \rho_\tau(Y_{ij}-X_i^\top\beta), \qquad \rho_\tau(u)=u\bigl(\tau-I[u<0]\bigr).$$3

Because this set is typically not a singleton, one paper constructs a representative point by choosing an ordering and factorizing the joint probability into sequential conditional probabilities, with $$\sum_{i=1}^n \rho_\tau(Y_{ij}-X_i^\top\beta), \qquad \rho_\tau(u)=u\bigl(\tau-I[u<0]\bigr).$$4 (Galvao et al., 21 Aug 2025).

Finally, multivariate M-quantile regression is a related but distinct development. It defines $$\sum_{i=1}^n \rho_\tau(Y_{ij}-X_i^\top\beta), \qquad \rho_\tau(u)=u\bigl(\tau-I[u<0]\bigr).$$5 through an influence function and associated loss, combines the robustness of quantile and the efficiency of expectile regression, and uses shared discrete random effects to induce dependence across outcomes and time in multivariate longitudinal data (Alfo' et al., 2016).

3. Estimation paradigms and computational architectures

The sign-concordance approach is intentionally modular. Step 1 estimates each marginal quantile regression separately. Step 2 fits a multinomial logistic model for the $$\sum_{i=1}^n \rho_\tau(Y_{ij}-X_i^\top\beta), \qquad \rho_\tau(u)=u\bigl(\tau-I[u<0]\bigr).$$6 possible sign patterns, with one pattern chosen as reference. In the bivariate case, this yields log-odds such as

$$\sum_{i=1}^n \rho_\tau(Y_{ij}-X_i^\top\beta), \qquad \rho_\tau(u)=u\bigl(\tau-I[u<0]\bigr).$$7

and the conditional correlation

$$\sum_{i=1}^n \rho_\tau(Y_{ij}-X_i^\top\beta), \qquad \rho_\tau(u)=u\bigl(\tau-I[u<0]\bigr).$$8

can be written in closed form from the multinomial probabilities (Columbu et al., 2021).

Bayesian directional models replace this two-step decomposition with a full hierarchical specification. In the structured-additive multiple-output framework, each directional regression includes basis expansions such as P-splines, spatial Markov random fields, or random effects, and the asymmetric-Laplace likelihood is represented through a normal-exponential mixture. MCMC estimation is then carried out by a block-Gibbs sampler, and a Gaussian-process adjustment across quantile levels enforces noncrossing directional quantiles and nested regions (Santos et al., 2019).

Bayesian nonparametric formulations move one level higher by modeling the entire conditional density. The dependent Dirichlet process approach places a DDP prior on a covariate-indexed mixing measure $$\sum_{i=1}^n \rho_\tau(Y_{ij}-X_i^\top\beta), \qquad \rho_\tau(u)=u\bigl(\tau-I[u<0]\bigr).$$9, convolves it with a countable mixture of multivariate normal kernels, uses a truncated stick-breaking representation for computation, and estimates parameters by a blocked Gibbs sampler. Quantiles are then extracted numerically from the fitted mixture at each posterior draw (Bhattacharya et al., 2020).

Likelihood-based joint estimation under the multivariate asymmetric Laplace (MAL) distribution is another major paradigm. In the maximum-likelihood formulation, the MAL is reparameterized so that the $Y$0th marginal has its $Y$1-quantile at $Y$2, and a location-scale mixture representation yields an EM algorithm with latent exponential variables. A penalized version of the EM adds an $Y$3 penalty for variable selection (Petrella et al., 2018). In the Bayesian time-varying volatility extension, the MAL mixture representation is combined with a Cholesky-type decomposition of the scale matrix and either stochastic volatility or GARCH dynamics, and estimation proceeds by MCMC (Iacopini et al., 2022).

High-dimensional multivariate response settings motivate low-rank and multitask estimators. In factorisable multitask quantile regression, one models each coefficient matrix $Y$4 at quantile level $Y$5 as low rank, estimates it by a nuclear-norm penalized convex program, smooths the nonsmooth empirical risk, and solves the resulting problem with FISTA and singular-value thresholding (Chao et al., 2015).

Neural parameterizations appear in both Euclidean and manifold-valued settings. Quantile surfaces are estimated in two stages: a deterministic point forecast produces the center, and a quantile surface regression neural network (QSNN) predicts directional radial quantiles over a grid of directions and quantile levels (Bieshaar et al., 2020). On manifolds, the conditional transport potential is parameterized by partial-input $Y$6-concave neural networks, soft minimum or entropic regularization makes the optimization differentiable, and stochastic Adam is used for large-scale training (Pegoraro et al., 2023).

Sequential conditional-CDF MQR occupies a middle ground between joint and marginal modeling. One estimates a sequence of univariate quantile regressions, but later stages use generated regressors built from indicators such as $Y$7. The paper develops both the unsmoothed estimator and a smoothed version to handle the non-differentiability of the generated indicators (Galvao et al., 21 Aug 2025). In longitudinal multivariate M-quantile regression, by contrast, the main computational device is an extended EM algorithm over finite-mixture random effects, with weighted Iteratively Reweighted Least Squares in the M-step and optional model selection by BIC or AIC (Alfo' et al., 2016).

4. Statistical properties, identifiability, and evaluation

Several recurring inferential issues structure the MQR literature. One is identifiability. In the sign-concordance model, identifiability requires selecting a reference sign pattern, so all $Y$8 are interpreted relative to that baseline (Columbu et al., 2021). In directional frameworks, noncrossing across quantile levels is central because crossings produce non-nested regions; the Gaussian-process post-processing step is designed precisely to guarantee ordered directional quantiles and nested halfspaces (Santos et al., 2019).

A second issue is the propagation of uncertainty through multi-stage procedures. The sign-concordance paper emphasizes that variability from the first-step marginal quantile fits must be propagated into inference for the second-step multinomial coefficients, and recommends a nonparametric bootstrap for valid standard errors and confidence intervals (Columbu et al., 2021). The sequential conditional-CDF approach derives a Bahadur representation for each stage, proves joint weak convergence for the estimated quantile graph, and shows root-$Y$9 asymptotic normality for the smoothed two-step estimator, with variance inflation from generated regressors accounted for in $\mathbb R^k$0 (Galvao et al., 21 Aug 2025).

Other formulations obtain consistency by modeling the full conditional law. The DDP-based geometric quantile method establishes posterior consistency for the conditional density and, by an arg-max-theorem argument, posterior consistency for the geometric quantile surface (Bhattacharya et al., 2020). The low-rank multitask estimator derives a non-asymptotic error bound in $\mathbb R^k$1 whose rate depends on the latent rank $\mathbb R^k$2, the covariate dimension $\mathbb R^k$3, the response dimension $\mathbb R^k$4, and the sample size $\mathbb R^k$5, while also controlling optimization error from smoothing and approximate numerical solution (Chao et al., 2015). In multivariate longitudinal M-quantile regression, standard errors can be obtained through a sandwich estimator based on Oakes’ identity and the individual scores (Alfo' et al., 2016).

Evaluation targets differ sharply across formulations because the estimands differ. The sign-concordance framework focuses on conditional tail association, especially the conditional correlation $\mathbb R^k$6 of exceedance indicators (Columbu et al., 2021). Quantile-surface models introduce directional CRPS, empirical coverage $\mathbb R^k$7 for calibration, and volume-based sharpness measures (Bieshaar et al., 2020). Manifold vector quantile regression evaluates contour coverage, KDE-$\mathbb R^k$8, ESS%, and the Jacobian-based conditional likelihood induced by the inverse transport map (Pegoraro et al., 2023). In adjacent work on simultaneous inference across multiple quantile levels, a multivariate extension of the rank-score test embedded in a closed-testing procedure controls the familywise error rate and attains higher power than Bonferroni in simulation (Santis et al., 11 Nov 2025).

5. Empirical domains and substantive findings

Applications in biostatistics and social science emphasize dependence structures that separate univariate quantile fits would miss. In lung function data with bivariate outcomes $\mathbb R^k$9, the sign-concordance method found that the estimated conditional correlation $\le$0 is high at $\le$1 across all covariate values, lower at $\le$2 in subjects with comorbidities, and declining at $\le$3 for young or very tall individuals; in the National Merit Twin Study, monozygotic twins had significantly higher log-odds of positive concordance and uniformly higher $\le$4 than dizygotic twins, especially at upper quantiles (Columbu et al., 2021).

In structured-additive directional models, the German inequality application used a bivariate response of standardized health score and log-income with 99 directions, and the Brazilian exam application used a trivariate response with 512 directions. The reported results show that Gaussian-process adjustment removed crossings in sparse-data regions and that the resulting quantile contours revealed interactions of age, education, family status, sex, region, school type, parental background, and gender that were not visible from mean-based summaries (Santos et al., 2019). In the nonparametric Bayesian geometric quantile model, simulation studies with heavy-tailed and skewed bivariate errors gave smaller MSE for the estimated spatial median than two competitors, and the blood-pressure application yielded smoothly increasing posterior spatial-median curves with age and coordinate-wise 95% credible bands described as robust to outliers (Bhattacharya et al., 2020).

Finance and macroeconomic applications exploit joint tails directly. In joint MAL-based estimation for Italian firms, simulation evidence showed RMSE reductions of up to 20–30% relative to separate univariate quantile regressions, especially at more extreme $\le$5, and the financial-distress application identified profit, financial expense, retained earnings, fixed assets, and cash flow as the main covariates associated with leverage and EBITDA quantiles (Petrella et al., 2018). In the energy-commodities study, MQR-SV and MQR-GARCH significantly outperformed the homoskedastic QVAR benchmark at tail quantiles $\le$6 for one- and five-day horizons, while quantile-score-based model combination yielded further gains when no single specification dominated (Iacopini et al., 2022). In the Argentina exchange-rate pass-through application, bivariate MQR on output and inflation traced quantile graphs over a fine $\le$7-grid and showed that larger devaluations shift the output-inflation curve to the right, with tail quantiles responding more sharply than central ones (Galvao et al., 21 Aug 2025).

Forecasting applications demonstrate the value of explicitly multivariate predictive distributions. In renewable-energy forecasting for two wind parks, QSNN used a dense grid of quantile levels and improved average directional CRPS by 80–90% over a homoscedastic Gaussian baseline at all distances; in cyclist trajectory forecasting, directional CRPS improved by up to 99% for short horizons and about 83% at 2.5 seconds, while Q–Q plots remained near the diagonal (Bieshaar et al., 2020).

Longitudinal studies use a different notion of quantile-like heterogeneity. In the Millennium Cohort Study, multivariate M-quantile regression for internalizing and externalizing SDQ scores found stronger positive effects of socio-economic disadvantage, maternal depression, and adverse life events at upper M-quantiles than at the median, larger between-subject than within-subject effects, and markedly non-Gaussian random-effect distributions, which supported the finite-mixture specification (Alfo' et al., 2016).

6. Limitations, misconceptions, and active directions

The most persistent misconception is that MQR is a single settled methodology. The literature instead presents multiple non-equivalent answers to the question of what a “multivariate quantile” should be: directional regions, geometric/spatial quantiles, star-shaped surfaces, sign-pattern probabilities, MAL-based joint marginal quantiles, sequential conditional-CDF contours, and manifold transport maps all target different objects (Santos et al., 2019, Bhattacharya et al., 2020, Bieshaar et al., 2020, Columbu et al., 2021, Pegoraro et al., 2023, Galvao et al., 21 Aug 2025). This suggests that methodological choice in MQR is primarily driven by the intended estimand and the geometry of the response space.

Scalability is a recurrent limitation. In sign-concordance models, the multinomial step has $\le$8 parameters, so dimensionality explodes with the number of outcomes and, in practice, $\le$9 or $\mathbb R$0 is most common (Columbu et al., 2021). Quantile surfaces were designed mainly for $\mathbb R$1 and require sampling many directions on $\mathbb R$2; the paper also states that QSNN currently models unimodal, star-shaped distributions only, with multimodal generalizations requiring mixtures of surfaces or non-starshaped constructions (Bieshaar et al., 2020). On manifolds, the theory assumes a complete Riemannian manifold with known $\mathbb R$3 and $\mathbb R$4 maps and conditional laws absolutely continuous with respect to the manifold volume measure (Pegoraro et al., 2023).

Model specification is another major issue. The sign-concordance paper notes that misspecification in the multinomial step can bias estimates of tail association (Columbu et al., 2021). MAL-based methods rely on a working likelihood, and one paper states explicitly that if the true error law deviates substantially, estimates may be inefficient; the same paper also notes that numerical solution for the scale-skew matrix can be costly for large $\mathbb R$5 and that EM can converge slowly in high dimensions (Petrella et al., 2018). In the time-varying volatility setting, no single specification uniformly outperformed the others across quantiles, time, or variables, motivating quantile-score-based model combination as a hedge against specification risk (Iacopini et al., 2022).

The choice of ordering in sequential conditional-CDF MQR is also consequential. The method fixes a permutation and constructs a representative point on the joint contour by sequential conditioning, rather than recovering the whole contour directly (Galvao et al., 21 Aug 2025). This suggests that the reported quantile graph is tied to the chosen conditioning order. In high-dimensional multitask settings, the low-rank approach raises further open problems already identified in the paper itself, including dependence beyond i.i.d. design, nonparametric factor mappings, extreme quantiles with small $\mathbb R$6, and faster solvers for very large $\mathbb R$7 and $\mathbb R$8 (Chao et al., 2015).

A final source of confusion is terminological. In some time-series literature, “MQR” refers to joint estimation of multiple quantiles for a scalar response, regularized across quantile levels, rather than regression with a multivariate target. The interquantile Lipschitz-regularized wind-power model is a prominent example of this distinct usage (Ruas et al., 2020). This suggests that, in current usage, the acronym alone is insufficient: one must distinguish multivariate-target quantile regression from multiple-quantile regression for scalar targets.