In this paper, we present an algorithm for Censored Quantile Regression (CQR) estimation problems. Our method permits CQR estimation problems to be solved more efficiently and reliably than was hitherto possible. It guarantees to find a high quality estimator in O(kÃ—nÂ²) operations with k regressors and n observations, which is much less than the existing algorithms for CQR problems.Cencored Quantile Regression

Dianqin Wang

Thanasis Stengos

English

Research Papers in Economics

An algorithm for censored quantile regressions 
Thanasis Stengos Dianqin Wang
University of Guelph University of Guelph
Abstract
In this paper, we present an algorithm for Censored Quantile Regression (CQR) estimation
problems. Our method permits CQR estimation problems to be solved more efficiently and
reliably than was hitherto possible. It guarantees to find a high quality estimator in O(k×n²)
operations with k regressors and n observations, which is much less than the existing
algorithms for CQR problems.
The first author would like to acknowledge financial support from the SSHRC of Canada.
Citation: Stengos, Thanasis and Dianqin Wang, (2007) "An algorithm for censored quantile regressions." Economics
Bulletin, Vol. 3, No. 1 pp. 1-9
Submitted: October 5, 2006.  Accepted: January 9, 2007.
URL: http://economicsbulletin.vanderbilt.edu/2007/volume3/EB-06C20071A.pdf1 Introduction
It is well known that censoring poses serious problems for regression models, see Goldberger
(1983) and Honoré (1992), as the resulting estimators in general may not be consistent.
However, in the context of quantile regression, censoring can be handled under very weak
distributional assumptions. For the Censored Quantile Regression (CQR) model, Powell
(1986) obtained estimator consistency and asymptotic normality. However, although the
CQR model has many appealing theoretical features, its application has been hampered by
the computational complexity of standard methods in obtaining CQR estimators.
The CQR estimation problem is to minimize a non-convex and piecewise linear distance
function which may have multiple local optima. The optimal estimators occur at non-smooth
points. The feature of the distance function presents a challenge to the standard optimization
techniques which require the objective function to be convex and/or diﬀerentiable. Several
algorithms have been developed by modifying standard optimization algorithms, yet their
performance has not been satisfactory. An optimal estimator for a CQR problem is diﬃcult to
obtain by modifying other existing algorithms, especially for large-scale problems in practice,
see Fitzenberger (1997).
Since the CQR model was introduced by Powell (1986), several algorithms have been
presented in the literature to deal with this problem. A detailed survey of the existing algo-
rithms and their performances may be found in Fitzenberger (1997) and Buchinsky (1994).
To name a few of the algorithms, they include a modiﬁed reduced-gradient algorithm by
Womersley (1986), an adaptation of the standard Barrodale-Roberts algorithm, see Barro-
dale and Roberts, (1973), for Quantile Regressions to the case of Censored Quantile Regres-
sions by Fitzenberger (1997), an interior point approach by Koenker and Park (1996), and
an iterative linear programming algorithm by Buchinsky (1994). Emulation algorithms (EA)
presented in Pinkse (1993) and Fitzenberger (1997) compute the exact CQR estimator by
checking every critical point. Some hybrid algorithms were proposed by Fitzenberger and
Winker (1999). However, these algorithms have run into diﬃculties in solving CQR esti-
mation problems, since they exhibit a high degree of complexity in their implementation.
Typically these algorithms achieve convergence to local minima, whereas obtaining a global
minimum requires a heavy computational load, something that renders their use in solving
real economic problems impractical. Among all these methods, only emulation algorithms
could be used to guarantee convergence.
2The purpose of this paper is to oﬀer an alternative algorithm that simpliﬁes the computa-
tions of the CQR estimation problem. The method presented in this paper uses a systematic
procedure to improve the reliability of the estimator, by using approximately O(k ×n2) op-
erations for a CQR problem with n observations and k regressors. The algorithm is simple
to implement. The next section of the paper presents the CQR estimation problem. We
then proceed to present the algorithm and then oﬀer some simulation results regarding its
relative performance when compared with the EA alternative.
2 The CQR estimation problem
The CQR estimation problem is to minimize a nonconvex and piecewise linear distance
function
min
n X
i=1
(θ ∗ I(di > 0) + (1 − θ) ∗ I(di < 0))|di| (1)
where di = yi − max(xiβ,cyi), and I(x) is an indication function, such that
I(x)=
(
1i f x i s t r u e
0 otherwise.
Also θ ∈ (0,1) i st h eq u a n t i l ea n dyi is censored from below at cyi.T h ev e c t o rβ is k ×1
dimensional. Other types of CQR can be easily transformed into the form of Equation (1).
We deﬁne Φ(β) as
Φ(β)=
n X
i=1
(θ ∗ I(di > 0) + (1 − θ) ∗ I(di < 0))|di|.
The distance function Φ(β) is piecewise linear, non-convex, and it has local minima, and
may not have a unique solution as illustrated in the following example. Consider the case
where θ =0 .5, and
Ã
x1
x2
!
=
Ã
3
4
!
,
Ã
y1
y2
!
=
Ã
1
3
!
,
Ã
cy1
cy2
!
=
Ã
1
2
!
.
The distance function is depicted below
30.8 0.75 0.7 0.65 0.6 0.55 0.5 0.45 0.4 0.35 0.3 0.25 0.2
0.8
0.775
0.75
0.725
0.7
0.675
0.65
0.625
0.6
0.575
0.55
0.525
0.5
ß ß
The distance function Φ(β)
The above ﬁgure shows that the distance function Φ(β) has a local minimum 5
8 at β = 3
4.
It is ﬂat when β ≤ 1
3, and stays at the minimum of zero. From the above ﬁgure, we see
that the CQR distance function is not smooth at the critical points β =
yi
xi, and
cyi
xi , for
i =1 ,···,n. Otherwise, Φ(β) is linear. The minimum of Φ(β) must occur at one of these
points.
3 A new algorithm for CQR estimation
As mentioned above CQR estimation deals with the minimization of a non-convex function.
Non-convex functions have been proven very diﬃcult to minimize or maximize. However,
from the example above we notice that CQR estimation involves a particular type of non-
convexity. Given a direction and a starting point, the objective function for CQR estimation
is a piecewise linear function. It consists of at most 2n +1line-segments, or there are at
most 2n critical points (knots). This unique property guarantees that a local minimum could
occur only in 2n critical points or knots.
Based on the above observations, our method is designed to ﬁnd the minimum point by
improving the value of the objective function at a given critical point with a smaller objective
function value at another critical point. In other words, the method oﬀers a way of moving
from a given critical point to the next critical point with a smaller objective function value.
This process continues until no critical points with smaller objective function value can be
found.
More speciﬁcally, our algorithm starts with an initial critical point, and then we choose
a direction. Next, we check the 2n critical points for a given direction. We choose the
critical point with the smallest objective function value to be the new starting point. Then
4we choose another direction, and we repeat the process, until no new starting point can be
found.
4 Implementation of the algorithm
Since the minimum only occurs at critical points, this motivates us to design an algorithm to
improve the estimated β through critical points, such that Φ(βl+1)is less than Φ(βl). In other
words, given a starting value βl, we try to ﬁnd βl+1,s u c ht h a tΦ(βl+1)is less than Φ(βl).
We proceed as follows. We ﬁrst deﬁne βl+1 = βl + dj ∗ ej for j ∈ (1,k), where ej is a
k−dimensional vector with its jth element being one, the rest zero and dj is a step size.
We next have xiβl+1 = xi (βl + dj ∗ ej). The potential critical points are those such that
xiβl+1 = yi for each i ∈ (1,n). That is xi (βl + dj ∗ ej)=yi, or xi,jdj = yi−xiβl.I fxi,j is not
zero, then we have dj =
yi−xiβl
xi,j . We proceed by computing Φ(βl+1). We deﬁne βl+1 = β
0
l+1 to
be that value that minimizes
©
Φ(βl+1)
ª
for all i ∈ (1,n) and j ∈ (1,k).I fΦ(βl+1)turns out
to be less that Φ(βl), then a new βl+1 is found. We repeat the above process by replacing βl
by βl+1, until a minimum point is reached. The process will end in ﬁnite steps since there
is only a ﬁnite number of critical points. To ﬁnd an improvement of βl involves two loops,
which need about O(k×n2) operations. The above critical point search algorithm is outlined
below1.
1. for a given βl.
2. for i =1 ,···,n.
(a) for j =1 ,···,k.
i. if xi,j is not zero, then
A. dj =
yi−xiβl
xi,j .
B. compute Φ(βl + dj × ej).
C. if Φ(βl + dj × ej) < Φ(βl), then βl+1 = βl + dj × ej.
ii. end of “if” loop.
(b) end of loop for j
3 .e n do fl o o pf o ri.
1The implementation of the above method is written in Gauss code and is available from the authors on
request.
54. if βl+1 is found during the loop for i,t h e n
(a) replace βl by βl+1.
(b) repeat steps 2 to 3.
5. otherwise the program terminates, a minimum is found.
5 Simulation Study
The timing of ﬁnding estimates of CQR problems is reported in Table 1. The method
presented in this paper is coded in Gauss, while the computer used to perform the simulations
is a Pentium 4 with 1.5 GHZ CPU. We also give the estimated times from using the EA of
Pinkse (1993) and Fitzenberger (1997), which is the only one among the standard methods
to guarantee convergence. We know that the EA methods checks
Ã
n
k
!
critical points.
The number of critical points checked by our method are reported in Table 2. Based on the
critical points checked, we can estimate the time that it would take the EA method to solve
the same CQR problems. To compare our method with the EA method, we compute the
time ratio of the EA and our method. We randomly generated 2000×20 uniform numbers in
the interval (−2.5,2.5). The CQR problems are solved with (n =5 0 0 ,1000,1500,2000) and
(k =2 ,5,10,15,20). The variable y is censored below a number which is randomly generated
between (−1,1). The dependent variable y is generated as follows:
y =m a x ( x1 + ···+ xk + ε,cy).
Here, the distribution of ε is the Student’s t-distribution with 5 degrees of freedom. The
results are reported in Table 2. The total time unit is seconds.
From Table 2, we see that the timing ratio increases dramatically as k grows from 2 to
20 for given n. It is a very diﬃcult task to ﬁnd the optimal points by using the EA method
for k over 10 using standard computers. On the other hand our proposed method can solve
large scale CQR problems that appear in practice using only O(k ×n2) operations, which is
a fraction of what is required by other existing methods.
6Table 1: The timing results for the method presented in this paper.
Our Method
n k Time (Seconds)
500 2 4.98E+01
500 5 7.13E+02
500 10 1.67E+04
500 15 1.43E+04
500 20 3.18E+04
1000 2 5.19E+01
1000 5 1.68E+04
1000 10 2.55E+04
1000 15 1.38E+07
1000 20 3.00E+08
1500 2 2.63E+03
1500 5 2.05E+04
1500 10 3.46E+08
1500 15 7.81E+08
1500 20 1.13E+09
2000 2 3.28E+03
2000 5 1.64E+04
2000 10 6.66E+08
2000 15 2.67E+09
2000 20 1.83E+09
7Table 2: The comparison of timings between our method and EA method. NCPC denotes
the number of critical points checked.
Our Method EA Method Time ratio of
n k NCPC NCPC EA over ours
500 2 105000 124750 1.19E+00
500 5 125000 2.55245E+11 2.04E+06
500 10 780000 2.45811E+20 3.15E+14
500 15 3427500 1.88779E+28 5.51E+21
500 20 3120000 2.6672E+35 8.55E+28
1000 2 58000 499500 8.61E+00
1000 5 375000 8.25029E+12 2.20E+07
1000 10 2330000 2.6341E+23 1.13E+17
1000 15 5715000 6.88141E+32 1.20E+26
1000 20 19860000 3.39483E+41 1.71E+34
1500 2 78000 1124250 1.44E+01
1500 5 412500 6.28604E+13 1.52E+08
1500 10 5730000 1.54203E+25 2.69E+18
1500 15 9675000 3.12155E+35 3.23E+28
1500 20 12420000 1.20351E+45 9.69E+37
2000 2 112000 1999000 1.78E+01
2000 5 800000 2.65336E+14 3.32E+08
2000 10 6280000 2.75899E+26 4.39E+19
2000 15 10590000 2.37736E+37 2.24E+30
2000 20 12560000 3.91816E+47 3.12E+40
8References
[1] Barrodale, I. and F. Roberts, (1973), “An Improved Algorithm for Discrete l1 Linear
Approximation,” SIAM Journal of Numerical Analysis 10, 839-848.
[2] Buchinsky, M., (1994), “Changes in the U.S. Wage Structure 1963-1987: Application of
Quantile Regression,” Econometrica 62, 405-458.
[3] Fitzenberger, B., (1997), “A Guide to Censored Quantile Regressions,” in Handbook
of Statistics, Robust Inference by G.S. Maddala and C.R. Rao eds., North-Holland,
Amsterdam.
[4] Fitzenberger, B. and P. Winker, (1999), “Improving the Computation of Censored
Quantile Regressions,” Universität Mannheim, Institute für VWL und Statistik, Band
Nr. 568-99.
[5] Goldberger, A.R., (1983), “Abnormal Selection Bias,” in Studies in Econometrics, Time
Series, and Multivariate Statistics by S. Karlin, ed., Academic Press, New York.
[6] Honoré, B. E., (1992), “Trimmed LAD and Least Squares Estimation of Truncated and
Censored Regression Models with Fixed Eﬀects,” Econometrica 60, 553-565.
[7] Koenker, R. and B.J. Park, (1996), “An Interior Point Algorithm for Nonlinear Quantile
Regression,” Journal of Econometrics 71, 265-283.
[8] Pinkse, C.A.P., (1993), “On the Computation of Semiparametric Estimates in Limited
Dependent Variable Models,” Journal of Econometrics 58, 185-205.
[9] Powell, J., (1986), “Censored Regression Quantiles,” Journal of Econometrics 32, 143-
155.
[10] Womersley, R.S., (1986), “Censored Discrete Linear l1 Approximation,” SIAM Journal
of Scientiﬁc and Statistical Computing 7, 105-122.
9

An algorithm for censored quantile regressions

http://www.accessecon.com/pubs/EB/2007/Volume3/EB-06C20071A.pdf

An algorithm for censored quantile regressions

Abstract

Similar works

Full text

Available Versions

Research Papers in Economics