Quantile Regression
USAGE: [p,stats]=quantreg(x,y,tau[,order,nboot]);
INPUTS:
x,y: data that is fitted. (x and y should be columns)
Note: that if x is a matrix with several columns then multiple
linear regression is used and the "order" argument is not used.
tau: quantile used in regression.
order: polynomial order. (default=1)
nboot: number of bootstrap surrogates used in statistical inference.(default=200)
stats is a structure with the following fields:
.pse: standard error on p. (not independent)
.pboot: the bootstrapped polynomial coefficients.
.yfitci: 95% confidence interval on polyval(p,x)
Note: uses bootstrap on residuals for statistical inference. (see help bootstrp)
check also: http://www.econ.uiuc.edu/~roger/research/intro/rq.pdf
EXAMPLE:
x=(1:1000)';
y=randn(size(x)).*(1+x/300)+(x/300).^2;
[p,stats]=quantreg(x,y,.9,2);
plot(x,y,x,polyval(p,x),x,stats.yfitci,'k:')
legend('data','2nd order 90th percentile fit','95% confidence interval','location','best')
For references on the method check e.g. and refs therein:
http://www.econ.uiuc.edu/~roger/research/rq/QRJEP.pdf
Copyright (C) 2008, Aslak Grinsted
Aslak Grinsted (2020). quantreg(x,y,tau,order,Nboot) (https://www.mathworks.com/matlabcentral/fileexchange/32115quantregxytauordernboot), MATLAB Central File Exchange. Retrieved .
1.3.0.0  implemented suggested change from Simeon Yurek in a FEX comment 

1.2.0.0  Fixed another small bug. 

1.1.0.0  Fixed a few issues with input parameter parsing. 
Create scripts with code, output, and formatted text in a single executable document.
Hi Aslak,
Thanks for the code.
However, I got the message like below. Do you know the reason?
"Exiting: Maximum number of function evaluations has been exceeded
 increase MaxFunEvals option.
Current function value: 334661.090416 "
I also found that Koenker  referenced here  translated his R code to matlab at one point. the link to that code is on this page: http://www.econ.uiuc.edu/~roger/research/rq/rq.html
The code is well written, but you should not optimize using fminsearch. Quantile regression can be framed as a linear programming problem, guaranteeing finding the optimum solution and doing so in little time. I recommend the R package quantreg or  if you need matlab  this stackoverflow post. https://stackoverflow.com/questions/21657158/quantileregressionwithlinproginmatlab
@Antony: I have the same question with you for a big dataset. But the results seem to be right.
Hi Aslak, I am inexperienced using quantile regression and was wondering if you or anyone could break down what the different numbers in the p statistic mean.
hi,Aslak.I have some questions to ask you.Now I want to use this code to estimate the parameters of model.For example.x1=a*x2+b*x3+c*x4+d*x2.^2+e*x3.^2+f*x4.^2...
Correction: Not coefficient estimates. I mean standard errors.
Hi Aslak, thanks for your contribution. I'm wondering why when I switch the order of [x1,x2] to [x2,x1], the coefficients estimates lose the rotation invariance property?
Hi Aslak, first thank you very much for sharing the code.
A quick question is about the way the standard errors of the parameter estimates are bootstrapped in your code. For many crosssectional data, it is not unreasonable to assume independence btw data points where it is fine to use the bootstrap method in your code. However, if one is dealing with time series data, there is usually serial correlation btw data points, a more suitable bootstrap method should be used, e.g. the stationary bootstrap of Politis&Romano(1994). I would suggest that the independence assumption is explicit in the code description and perhaps have another code specifically for correlated data points. Best regards, Zhongfang
@lisa: they are the polynomial coefficients. They can be interpreted similar to the ouput from polyval
Hi, very helpful code, thank you.
A quick question: p outputs several values depending on the order, what do they represent?
thanks!
@SimeonYurek: thank you for the suggestion. It has been implemented.
hi. just a quick question. How do we calculate the goodness of fit for quantile regression? thanks
Very nice code for Koenker and Hallock (2001). Thanks for posting.
One question: in your statement of the function rho (line 85), when r >= 0 (all positive residuals above x*p), does the function reduce to abs(r), and thus is not weighted by tau? It's true that [r  0.*r/tau] = r.*tau, but is the tau lost (being multiplied by zero)? Could you have stated instead:
rho=@(r)sum(abs(r).*abs(tau(r<0)))
and this way weight both over and under residuals? The stats for the bootstrap are slightly less robust but not by much. Please let me know if I'm off.
Antony: help fminsearch. 'MaxFunEvals' and 'MaxIter' can be defined as options.
Hi, first of all, thank you very much for the code.
When I run the code I get the following message:
"Exiting: Maximum number of function evaluations has been exceeded
 increase MaxFunEvals option.
Current function value: 3183.464509 "
I get this message a number of times with just the current function value changing.
While it is not an error message and I still obtain the estimations I need, I was wondering if this had any influence on the validity of my results.
Thank you in advance for anyone who can help me out here.
Excellent. Well written help and code! Runs as advertised!
If X and y have some missing observations, then
[p,stats]=quantreg(X,y,.9); does not work. I got error message.
Could you please suggest me how to resolve this problem??
@mohammad: The code does do MLR. Here's an example with multiple predictors:
X=randn(100,4);
y=X*[4 2 4 1]'+randn(100,1)*.2;
[p,stats]=quantreg(X,y,.9);
In fact I wanted to run this Mfile for multi column X, but I got error, so it is not implement for MLR? if it is , could you please put an example? because when I do not mention Order , I get error when i mention it is not MLR
Thanks for the comment, and suggestions for improvement. I've just uploaded a new version (it should be online shortly).
One unexpected thing with this code. Suppose I have a Y vector and want to regress it on one explanatory variable, but also include a constant in my regression. Then my X matrix has 2 columns, the first column is just ones, the second is the explanatory variable. It should still be possible to plot this along with the results from the regression.
Currently, however, if you run:
quantreg([ones(length(x),1) x], y,.5)
You get an error because it tries to plot it, but order has been set (line 44) to [], so things get messed up.
Additionally, the plots that do get produced look odd because the default is to draw lines between all of the points, which usually isn't what you want. For example, this doesn't look like what it should:
x=randn(1000,1);
y=1+5*x+randn(1000,1);
quantreg(x,y,0.5)
Finally, there's a problem with the error checking of inputs (lines 4147) because it's all one big ifelse statement. If you only put in three inputs, then line 44 runs and order gets set to [], but then the program exits the ifelse statement and so Nboot never gets set correctly. You should split these up into separate if statements because if there are only 3 inputs then you need to set both order and Nboot.