What type of data should I generate to observe/amplify a crossing problem in quantile regression? 1. Background
Crossing problem in quantile regression can be observed when we want to estimate several conditional quantiles (e.g. τ = 0.1, 0.2, . . . , 0.9), as two or more estimated
conditional quantile functions could cross or overlap. However, a higher percentile must have more points below it than a lower percentile. That is, we should not have more points below the 75th percentile than we do at the 80th percentile. This would violate the basic definition of percentiles, making the quantile regression invalid.
(Takeuchi et al., 2005 ; Saleh, 2021)
2. Question
This being said, I would like to generate synthetic data to observe and maximise this issue. To perform the quantile regression, I am using the quantreg framework in R.
(https://cran.r-project.org/web/packages/quantreg/quantreg.pdf)
Since starting my research, I have tried to generate various datasets that present mean or SD non-linearly changes and have yet to be successful.
Does anyone know how to do this, or any direction to look for?
 A: If we are to avoid trivialities related to outlying values, evidently crossing is an indication the regression model isn't correct.  So, all you need is suitably poor goodness of fit along with a little leverage and not too many data points.
For example, generate conditional distributions, some of whose quantiles vary nonlinearly with the explanatory variables, and fit linear models to the quantiles.  Something like $Y\mid x \sim\mathcal{N}(0,\exp(x))$ ought to do the trick.
Here is an example, zoomed at the right to show the details of the crossing regression lines.  It was created by setting df = 1 at the beginning of the R code below, which generates 601 data points.

The colors of the fitted lines start dark at low quantiles and graduate to lighter colors at higher quantiles.  They are in the correct order for $x gt 2$ or so, but cross at lower values of $x.$
The crossings are severe: ultimately, near $x=0,$ the fitted $5^\text{th}$ percentile exceeds the fitted $95^\text{th}$ percentile!
Even more complex regression models are subject to this problem.  Here are the same data fitted with a natural spline with three degrees of freedom (using the code below unchanged).

The nonlinearity is achieved through the strong variation of the conditional response variance, as suggested above.  The points at the right are given greater leverage by sampling them disproportionately more often.  This causes the trends at the right to overwhelm the fewer amounts of data at the left.
These patterns are consistent with other datasets (comment out the set.seed call to see).  In most cases the software does not issue a warning that the fit might not be unique.
library(quantreg)
#
# Generate a dataset.
#
set.seed(17)
n <- 601                  # Size of dataset
lambda <- 5               # Upper limit of `x`
N <- 31                   # Number of distinct possible values of `x`, equally spaced
probs <- seq_len(N) ^ 2.5 # Relative sampling probabilities for `x`
df <- 3                   # Spline degrees of freedom (1 = linear)
X <- data.frame(x = sort(sample(seq(0, lambda, length.out = N), 
                                n, replace = TRUE, prob = probs)))
X$y <- with(X, rnorm(n, sd = exp(x)))
#
# Fit a range of quantiles, not too extreme.
#
tau <- c(0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95) # quantiles
fit <- rq(y ~ splines:::ns(x, df), tau = tau, X)
#
# Compute the predicted values.
#
x.0 <-seq(0, lambda, length.out = 501)
y.hat <- predict(fit, newdata = data.frame(x = x.0))
#
# Plot the data and the results.
#
hue <- function(x) log10(1/2 - abs(x - 1/2))
hmax <- max(abs(hue(tau)))
par(mfrow = c(1, 2))
for (mu in c(1, 1/10))
with(X, {
  plot(range(x), range(c(y, y.hat)) * mu, type = "n", 
       main = bquote(Quantiles * phantom(.) * .(paste(tau, collapse = ", "))),
       cex.main = 1,
       xlab = "x", ylab = "y")
  invisible(sapply(seq_len(ncol(y.hat)), function(i) 
    lines(x.0, y.hat[,i], 
          col = hsv(-0.8 / hmax * hue(tau[i]), .9, (1 + tau[i])/2), lwd = 2)))
  points(x, y, col = "#00000060")
})
par(mfrow = c(1, 1))

