C This is a very primative version of the barrodale-roberts rq
C algorithm intended only to be used for resampling procedures.
C bug reports should be sent to roger@ysidro.econ.uiuc.edu or to
C
C Roger Koenker
C Department of Economics
C University of Illinois
C Champaign, Illinois, 61820
C
C
SUBROUTINE RQ1(M,N,M5,N2,A,B,TAU,TOLER,IFT,X,E,S,WA,WB,
*NSOL,SOL,DSOL,LSOL)
C
C Number of Observations
C Number of Parameters
C M+5
C N+2
C A is the X matrix
C B is the y vector
C tau, the desired quantile
C toler, smallest detectable |x-y|/x machine precision to the 2/3
C ift exit code:
C 0-ok
C else dimensions inconsistent
C X the parameter estimate betahat
C E is the residual vector
C S is an integer work array
C WA is a work array
C WB is another work array
C NSOL is an estimated (row) dimension of the solution array say 3*M
C SOL is the primal solution array
C DSOL is the dual solution arry
C LSOL is the actual dimension of the solution arrays
C Utilization: If you just want a solution at a single quantile you
C needn't bother with SOL, NSOL, etc, if you want all the solutions
C then set TAU to something <0 and SOL and DSOL will return all the
C estimated quantile solutions.
C The algorithm is a slightly modified version of Algorithm AS 229
C described in Koenker and D'Orey, "Computing Regression Quantiles,
C Applied Statistics, pp. 383-393.
C
INTEGER I,J,K,KL,KOUNT,KR,L,LSOL,M,M1,M2,M3,M4,M5,IFT
INTEGER N,N1,N2,NSOL,OUT,S(M)
LOGICAL STAGE,TEST,INIT,IEND
DOUBLE PRECISION A1,AUX,B1,BIG,D,DIF,PIVOT,SMAX,TAU,T0,T1,TNT
DOUBLE PRECISION MIN,MAX,TOLER,ZERO,HALF,ONE,TWO
DOUBLE PRECISION B(M),SOL(N2,NSOL),A(M,N),X(N),WA(M5,N2),WB(M)
DOUBLE PRECISION SUM,E(M),DSOL(M,NSOL)
DATA BIG/1.D37/
DATA ZERO/0.D00/
DATA HALF/0.5D0/
DATA ONE/1.D0/
DATA TWO/2.D0/
C
C CHECK DIMENSION PARAMETERS
C
IFT=0
IF(M5 .NE. M+5)IFT = 3
IF(N2 .NE. N+2)IFT = 4
IF(M.LE.ZERO.OR.N.LE.ZERO)IFT = 5
IF(IFT .GT. TWO)RETURN
C
C INITIALIZATION
C
M1 = M+1
N1 = N+1
M2 = M+2
M3 = M+3
M4 = M+4
DO 2 I=1,M
WB(I)=B(I)
DO 1 J=1,N
WA(I,J)=A(I,J)
1 CONTINUE
2 CONTINUE
WA(M2,N1)=ZERO
DIF = ZERO
IEND = .TRUE.
IF(TAU .GE. ZERO .AND. TAU .LE. ONE)GOTO 3
T0 = ONE/FLOAT(M)-TOLER
T1 = ONE - T0
TAU = T0
IEND = .FALSE.
3 CONTINUE
INIT = .FALSE.
LSOL = 1
KOUNT = 0
DO 9 K=1,N
WA(M5,K) = ZERO
DO 8 I=1,M
WA(M5,K) = WA(M5,K) + WA(I,K)
8 CONTINUE
WA(M5,K) = WA(M5,K)/FLOAT(M)
9 CONTINUE
DO 10 J=1,N
WA(M4,J) = J
X(J) = ZERO
10 CONTINUE
DO 40 I=1,M
WA(I,N2) = N+I
WA(I,N1) = WB(I)
IF(WB(I).GE.ZERO)GOTO 30
DO 20 J=1,N2
WA(I,J) = -WA(I,J)
20 CONTINUE
30 E(I) = ZERO
40 CONTINUE
DO 42 J=1,N
WA(M2,J) = ZERO
WA(M3,J) = ZERO
DO 41 I=1,M
AUX = SIGN(ONE,WA(M4,J)) * WA(I,J)
WA(M2,J) = WA(M2,J) + AUX * (ONE - SIGN(ONE,WA(I,N2)))
WA(M3,J) = WA(M3,J) + AUX * SIGN(ONE,WA(I,N2))
41 CONTINUE
WA(M3,J) = TWO * WA(M3,J)
42 CONTINUE
GOTO 48
43 CONTINUE
LSOL = LSOL + 1
DO 44 I=1,M
S(I) = ZERO
44 CONTINUE
DO 45 J=1,N
X(J) = ZERO
45 CONTINUE
C
C COMPUTE NEXT TAU
C
SMAX = TWO
DO 47 J=1,N
B1 = WA(M3,J)
A1 = (-TWO - WA(M2,J))/B1
B1 = -WA(M2,J)/B1
IF(A1 .LT. TAU)GOTO 46
IF(A1 .GE. SMAX) GOTO 46
SMAX = A1
DIF = (B1 - A1 )/TWO
46 IF(B1 .LE. TAU) GOTO 47
IF(B1 .GE. SMAX)GOTO 47
SMAX = B1
DIF = (B1 - A1)/TWO
47 CONTINUE
TNT = SMAX + TOLER * (ONE + ABS(DIF))
IF(TNT .GE. T1 + TOLER)IEND = .TRUE.
TAU = TNT
IF(IEND)TAU = T1
48 CONTINUE
C
C COMPUTE NEW MARGINAL COSTS
C
DO 49 J=1,N
WA(M1,J) = WA(M2,J) + WA(M3,J) * TAU
49 CONTINUE
IF(INIT) GOTO 265
C
C STAGE 1
C
C DETERMINE THE VECTOR TO ENTER THE BASIS
C
STAGE=.TRUE.
KR=1
KL=1
70 MAX=-ONE
DO 80 J=KR,N
IF(ABS(WA(M4,J)).GT.N)GOTO 80
D=ABS(WA(M1,J))
IF(D.LE.MAX)GOTO 80
MAX=D
IN=J
80 CONTINUE
IF(WA(M1,IN).GE.ZERO)GOTO 100
DO 90 I=1,M4
WA(I,IN)=-WA(I,IN)
90 CONTINUE
C
C DETERMINE THE VECTOR TO LEAVE THE BASIS
C
100 K=0
DO 110 I=KL,M
D=WA(I,IN)
IF(D.LE.TOLER)GOTO 110
K=K+1
WB(K)=WA(I,N1)/D
S(K)=I
TEST=.TRUE.
110 CONTINUE
120 IF(K.GT.0)GOTO 130
TEST=.FALSE.
GOTO 150
130 MIN=BIG
DO 140 I=1,K
IF(WB(I).GE.MIN)GOTO 140
J=I
MIN=WB(I)
OUT=S(I)
140 CONTINUE
WB(J)=WB(K)
S(J)=S(K)
K=K-1
C
C CHECK FOR LINEAR DEPENDENCE IN STAGE 1
C
150 IF(TEST.OR..NOT.STAGE)GOTO 170
DO 160 I=1,M4
D=WA(I,KR)
WA(I,KR)=WA(I,IN)
WA(I,IN)=D
160 CONTINUE
KR=KR+1
GOTO 260
170 IF(TEST)GOTO 180
WA(M2,N1)=TWO
GOTO 390
180 PIVOT=WA(OUT,IN)
IF(WA(M1,IN)-PIVOT-PIVOT.LE.TOLER)GOTO 200
DO 190 J=KR,N1
D=WA(OUT,J)
WA(M1,J)=WA(M1,J)-D-D
WA(M2,J)=WA(M2,J)-D-D
WA(OUT,J)=-D
190 CONTINUE
WA(OUT,N2)=-WA(OUT,N2)
GOTO 120
C
C PIVOT ON WA(OUT,IN)
C
200 DO 210 J=KR,N1
IF(J.EQ.IN)GOTO 210
WA(OUT,J)=WA(OUT,J)/PIVOT
210 CONTINUE
DO 230 I=1,M3
IF(I.EQ.OUT)GOTO 230
D=WA(I,IN)
DO 220 J=KR,N1
IF(J.EQ.IN)GOTO 220
WA(I,J)=WA(I,J)-D*WA(OUT,J)
220 CONTINUE
230 CONTINUE
DO 240 I=1,M3
IF(I.EQ.OUT)GOTO 240
WA(I,IN)=-WA(I,IN)/PIVOT
240 CONTINUE
WA(OUT,IN)=ONE/PIVOT
D=WA(OUT,N2)
WA(OUT,N2)=WA(M4,IN)
WA(M4,IN)=D
KOUNT=KOUNT+1
IF(.NOT.STAGE)GOTO 270
C
C INTERCHANGE ROWS IN STAGE 1
C
KL=KL+1
DO 250 J=KR,N2
D=WA(OUT,J)
WA(OUT,J)=WA(KOUNT,J)
WA(KOUNT,J)=D
250 CONTINUE
260 IF(KOUNT+KR.NE.N1)GOTO 70
C
C STAGE 2
C
265 STAGE=.FALSE.
C
C DETERMINE THE VECTOR TO ENTER THE BASIS
C
270 MAX=-BIG
DO 290 J=KR,N
D=WA(M1,J)
IF(D.GE.ZERO)GOTO 280
IF(D.GT.(-TWO))GOTO 290
D=-D-TWO
280 IF(D.LE.MAX)GOTO 290
MAX=D
IN=J
290 CONTINUE
IF(MAX.LE.TOLER)GOTO 310
IF(WA(M1,IN).GT.ZERO)GOTO 100
DO 300 I=1,M4
WA(I,IN)=-WA(I,IN)
300 CONTINUE
WA(M1,IN)=WA(M1,IN)-TWO
WA(M2,IN)=WA(M2,IN)-TWO
GOTO 100
C
C COMPUTE QUANTILES
C
310 CONTINUE
DO 320 I=1,KL-1
K=WA(I,N2)*SIGN(ONE,WA(I,N2))
X(K) = WA(I,N1) * SIGN(ONE,WA(I,N2))
320 CONTINUE
SUM=ZERO
DO 330 I=1,N
SUM=SUM + X(I) * WA(M5,I)
SOL(I+2,LSOL)=X(I)
KD=ABS(WA(M4,I))-N
DSOL(KD,LSOL)=ONE+WA(M1,I)/TWO
IF(WA(M4,I).LT.ZERO)DSOL(KD,LSOL)=ONE-DSOL(KD,LSOL)
330 CONTINUE
DO 335 I=KL,M
KD=ABS(WA(I,N2))-N
IF(WA(I,N2).LT.ZERO)DSOL(KD,LSOL)=ZERO
IF(WA(I,N2).GT.ZERO)DSOL(KD,LSOL)=ONE
335 CONTINUE
SOL(1,LSOL) = SMAX
SOL(2,LSOL) = SUM
IF(IEND)GOTO 340
INIT = .TRUE.
GOTO 43
340 CONTINUE
IF(LSOL.LE.2)GOTO 355
SOL(1,1)=ZERO
SOL(1,LSOL)=ONE
DO 350 I=1,M
DSOL(I,1)=ONE
DSOL(I,LSOL)=ZERO
350 CONTINUE
355 CONTINUE
L =KL-1
DO 370 I=1,L
IF(WA(I,N1).GE.ZERO)GOTO 370
DO 360 J=KR,N2
WA(I,J)=-WA(I,J)
360 CONTINUE
370 CONTINUE
WA(M2,N1)=ZERO
IF(KR.NE.1)GOTO 390
DO 380 J=1,N
D=ABS(WA(M1,J))
IF(D.LE.TOLER.OR.TWO-D.LE.TOLER)GOTO 390
380 CONTINUE
WA(M2,N1)=ONE
390 SUM = ZERO
DO 400 I=KL,M
K = WA(I,N2) * SIGN(ONE,WA(I,N2))
D = WA(I,N1) * SIGN(ONE,WA(I,N2))
SUM = SUM + D * SIGN(ONE,D) * (HALF + SIGN(ONE,D)*(TAU-HALF))
K=K-N
E(K)=D
400 CONTINUE
WA(M2,N2)=KOUNT
WA(M1,N2)=N1-KR
WA(M1,N1) = SUM
RETURN
END