function beta = ipqr (q, p, y, x, func, toler, beta)
% Use interior point algorithm for quantile regression.
% (See Koenker and Park (1996), J. Econometrics)
%
% Inputs:
% -- q is the quantile desired (number in (0,1)).
% -- p is the size of the parameter vector. This is actually
% redundant since it's just the length of beta, which is also
% input, but I didn't change the code to eliminate this
% redundancy. Note that p is not necessarily the number of columns
% of x, though this is usually the case.
% -- y is the response vector (n x 1).
% -- x is the predictor matrix (n x m).
% -- func is a character array containing the name of a function
% that evaluates f(x, beta). This MATLAB function takes (x, beta)
% as inputs where x is the n x m predictor matrix and beta is some
% p x 1 vector of parameter values. It returns an n x 1 vector.
% In the case of linear QR, m=p and the function (whatever it's
% called) returns x*beta.
% There must also be a second MATLAB function, with the same name
% as func except with the letter d appended to the front, that
% also takes (x,beta) as input but returns the $n x p$ matrix
% of partial derivatives of f(x,beta) evaluated at beta. In the
% linear QR example, m=p and this second function simply returns x.
% -- toler is a small number giving the minimum change in the value
% of the surrogate function before convergence is declared.
% -- beta is the starting value of beta (p x 1 vector).
flops(0); % reset flops count.
dfunc = ['d' func]; % name of differential function.
eta=.97; % Value recommended by K&P
omq=1-q;
flops2=0;
iteration = 0;
n=length(y);
d=zeros(n,1);
change = realmax;
residuals = y-feval(func,x,beta);
lastobj=q*sum(residuals)-sum(residuals(residuals<0));
while abs(change) > toler
iteration = iteration + 1;
J=feval(dfunc,x,beta); % J is first differential matrix.
dhat=d-J*(inv(J'*J)*(J'*d));
m=toler+max([-dhat./(1-q);dhat./q]);
d = dhat ./ m;
% Now we perform the inner iterations. K&P recommend two of them per
% outer iteration.
for k=1:2
dmatrix = min([q-d 1-q+d]')';
d2=dmatrix.^2;
direc = inv(J'*(d2(:,ones(p,1)).*J))*(J'*(d2.*residuals));
s = d2.*(residuals - J*direc);
alpha = max(max([s./(q-d) -s./(1-q+d)]'));
d=d+eta/alpha*s;
end
% Now we do the line search recommended by K&P. Some of the parameters
% passed to the fmin function are ad hoc (K&P give no specific
% recommendations).
step=fmin('ipqr_objfunc',-1,1,[],func,q,y,x,beta,direc) * direc;
beta = beta + step;
residuals = y - feval(func, x, beta);
obj=q*sum(residuals)-sum(residuals(residuals<0));
change=lastobj-obj;
lastobj=obj;
end
% At this point, beta is the estimate of regression coefficients;
% lastobj is the value of the objective function at beta;
% flops is the number of floating-point operations needed
% (according to the way MATLAB counts them, at least).