Solved – The issue of quantile curves crossing each other

Question

While trying to perform some quantile regressions on some data I have encountered a classic problem: the 0.25,0.50, 0.75, 0.95 quantile curves are crossing each other (I am dealing with linear curves).

I am currently working with Stata using the sqreg command.

I have read Howard Bondell's paper "Non-crossing quantile regression curve estimation"(2010) and I know he used R to code a solution.

I have also read papers from Andriyana and from Chernozhukov (on the issue of non-monotonicity issues in quantile regressions).

Just to be clear about what we are talking, here below is an example of quantile curves that DO NOT CROSS (and my issue is what to do when they cross):

I was wondering if any of you could recommend what software solutions are most commonly used in the working practice, to deal with crossing quantile curves (=to avoid their crossing), I know there are various methods, so the question is both: what method you prefer and what software would you to use to solve the issue?

Note: I am not posting a sample of my data because it's irrelevant, I am just interested in the general theory of the problem and how to deal with it specifically with software (ideally STATA or Mathematica, or R or Python or Perl).

Answer 1

In general case, the problem of non-crossing linear quantile regresion has no solutions.

Indeed, either your quantile lines are parallel (and trivial), or they do cross somewhere.

By placing restrictions, you could ensure that the lines do not cross in the training sample, but it will in no way guarantee that crossing will not occur at the next observation you see after training the model. And if the quantile lines intersect within the training sample, it most probably means that your model is specified incorrectly: either mean or standard deviation change nonlinearly, or you apply a wrong cost function when fitting the model.

If you still want to estimate such a model, I would propose fitting it as a neural network.

Your model should take input $X$, multiply it by matrix(!) $\beta$ to get a matrix of forecasts $f=X\beta$ with size $n \times k$, where $n$ is number of observations and $k$ is number of estimated quantiles. I assume that quantile percentages $q$ are in increasing order. You should minimize the function

$$ L = \sum_{i=1}^n \left(\sum_{j=1}^k \max(q_j(y_i-f_{ij}), (q_j-1)(y_i-f_{ij})) + \sum_{j=1}^{k-1} \alpha \max(0, \delta - (f_{i,j+1}-f_{ij})) \right) $$

The first term in the inner sum is just the ordinary quantile regression loss. The second term is the penalty that is applied if two consecutive quantile predictions differ by less then $\delta$.

Minimizing this function by gradient descent will give you your non-crossing quantile lines (if $\alpha$ is large enough). But I still warn you that if "natural" quantile lines intersect, then ther might be problems with the functional form of your model. Maybe you would prefer quantile estimates of Random Forest (like quantregForest in R), which are always consistent.

There is a Python example. As for me, both restricted and unrestricted versions look ugly.

# import everything
import keras
from keras import backend as K
from keras import Sequential
from keras.layers import Dense
import numpy as np
import matplotlib.pyplot as plt

# create the dataset
n = 200
np.random.seed(1)
X = np.sort(np.random.normal(size=n))[:, np.newaxis] + 4
y = 5 + X.ravel()*(1 + np.random.normal(size=n)*0.2)
quantiles = np.array([0.1, 0.25, 0.5, 0.75, 0.9])

# define loss function
def quantile_ensemble_loss(q, y, f, alpha=100, margin=0.1):
    error = (y - f)
    quantile_loss = K.mean(K.maximum(q*error, (q-1)*error))
    diff = f[:, 1:] - f[:, :-1]
    penalty = K.mean(K.maximum(0.0, margin - diff)) * alpha
    return quantile_loss + penalty

# fit two models
for i, alpha, name in [(1, 0, 'w/o penalty'), (2, 10, 'with_penalty')]:
    model = Sequential()
    model.add(Dense(len(quantiles), input_dim=X.shape[1]))
    model.compile(loss=lambda y,f: quantile_ensemble_loss(quantiles,y,f,alpha), optimizer=keras.optimizers.RMSprop(lr=0.003))
    model.fit(X, y, epochs=3000, verbose=0, batch_size=100);
    plt.subplot(1,2,i)
    plt.scatter(X.ravel(), y, s=1)
    plt.plot(X.ravel(), model.predict(X))
    plt.title(name)

plt.show()