1      SUBROUTINE RQ0(M,N,M5,N2,A,B,T,TOLER,IFT,X,E,S,WA,WB)
2C
3C     Modified to remove SOL and related vars -- only good for single tau
4C     M Number of Observations
5C     N Number of Parameters
6C     M5 = M+5  row dimension for WA
7C     N2 = N+2  col dimension for WA
8C     A is the X matrix
9C     B is the y vector
10C     T, the desired quantile
11C     TOLER, smallest detectable |x-y|/x machine precision to the 2/3
12C     IFT exit code:
13C		0-ok
14C		else dimensions inconsistent or T not in (0,1)
15C     X the parameter estimate betahat
16C     E is the residual vector
17C     S is an integer work array (M)
18C     WA is a real work array (M5,N2)
19C     WB is another real work array (M)
20C     Utilization:  If you just want a solution at a single quantile you
21C     The algorithm  is a slightly modified version of Algorithm AS 229
22C     described in Koenker and D'Orey, "Computing Regression Quantiles,
23C     Applied Statistics, pp. 383-393.
24C
25      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
26      INTEGER I,J,K,KL,KOUNT,KR,M,M1,M2,M3,M4,M5,IFT
27      INTEGER N,N1,N2,OUT,S(M)
28      LOGICAL STAGE,TEST,INIT,IEND
29      DOUBLE PRECISION MIN,MAX
30      DOUBLE PRECISION B(M),A(M,N),X(N),WA(M5,N2),WB(M),E(M)
31      DATA BIG/1.D37/
32      DATA ZERO/0.00D0/
33      DATA HALF/0.50D0/
34      DATA ONE/1.00D0/
35      DATA TWO/2.00D0/
36C
37C  CHECK DIMENSION PARAMETERS
38C
39      IFT=0
40      IF(N2 .NE. N+2)IFT = 4
41      IF(M.LE.ZERO.OR.N.LE.ZERO)IFT = 5
42      IF(IFT .GT. TWO)RETURN
43C
44C  INITIALIZATION
45C
46      M1 = M+1
47      N1 = N+1
48      M2 = M+2
49      M3 = M+3
50      M4 = M+4
51      M5 = M+5
52      DO 2 I=1,M
53      WB(I)=B(I)
54      DO 1 J=1,N
55      WA(I,J)=A(I,J)
56  1   CONTINUE
57  2   CONTINUE
58      WA(M2,N1)=ZERO
59      DIF = ZERO
60      IEND = .TRUE.
61      IF(T .GE. ZERO .AND. T .LE. ONE)GOTO 3
62      IFT = 6
63      RETURN
64  3   CONTINUE
65      INIT = .FALSE.
66      KOUNT = 0
67      DO 9 K=1,N
68      WA(M5,K) = ZERO
69      DO 8 I=1,M
70      WA(M5,K) = WA(M5,K) + WA(I,K)
71  8   CONTINUE
72      WA(M5,K) = WA(M5,K)/FLOAT(M)
73  9   CONTINUE
74      DO 10 J=1,N
75      WA(M4,J) = J
76      X(J) = ZERO
77 10   CONTINUE
78      DO 40 I=1,M
79      WA(I,N2) = N+I
80      WA(I,N1) = WB(I)
81      IF(WB(I).GE.ZERO)GOTO 30
82      DO 20 J=1,N2
83      WA(I,J) = -WA(I,J)
84 20   CONTINUE
85 30   E(I) = ZERO
86 40   CONTINUE
87      DO 42 J=1,N
88      WA(M2,J) = ZERO
89      WA(M3,J) = ZERO
90      DO 41 I=1,M
91      AUX = SIGN(ONE,WA(M4,J)) * WA(I,J)
92      WA(M2,J) = WA(M2,J) + AUX * (ONE - SIGN(ONE,WA(I,N2)))
93      WA(M3,J) = WA(M3,J) + AUX * SIGN(ONE,WA(I,N2))
94 41   CONTINUE
95      WA(M3,J) = TWO * WA(M3,J)
96 42   CONTINUE
97      GOTO 48
98 43   CONTINUE
99      DO 44 I=1,M
100      S(I) = ZERO
101 44   CONTINUE
102      DO 45 J=1,N
103      X(J) = ZERO
104 45   CONTINUE
105C
106C  COMPUTE NEXT T
107C
108      SMAX = TWO
109      DO 47 J=1,N
110      B1 =  WA(M3,J)
111      A1 = (-TWO - WA(M2,J))/B1
112      B1 = -WA(M2,J)/B1
113      IF(A1 .LT. T)GOTO 46
114      IF(A1 .GE. SMAX) GOTO 46
115      SMAX = A1
116      DIF = (B1 - A1 )/TWO
117 46   IF(B1 .LE. T) GOTO 47
118      IF(B1 .GE. SMAX)GOTO 47
119      SMAX = B1
120      DIF = (B1 - A1)/TWO
121 47   CONTINUE
122      TNT = SMAX + TOLER * (ONE + ABS(DIF))
123      IF(TNT .GE. T1 + TOLER)IEND = .TRUE.
124      T = TNT
125      IF(IEND)T = T1
126 48   CONTINUE
127C
128C COMPUTE NEW MARGINAL COSTS
129C
130      DO 49 J=1,N
131      WA(M1,J) = WA(M2,J) + WA(M3,J) * T
132 49   CONTINUE
133      IF(INIT) GOTO 265
134C
135C STAGE 1
136C
137C DETERMINE THE VECTOR TO ENTER THE BASIS
138C
139      STAGE=.TRUE.
140      KR=1
141      KL=1
142 70   MAX=-ONE
143      DO 80 J=KR,N
144      IF(ABS(WA(M4,J)).GT.N)GOTO 80
145      D=ABS(WA(M1,J))
146      IF(D.LE.MAX)GOTO 80
147      MAX=D
148      IN=J
149 80   CONTINUE
150      IF(WA(M1,IN).GE.ZERO)GOTO 100
151      DO 90 I=1,M4
152      WA(I,IN)=-WA(I,IN)
153 90   CONTINUE
154C
155C DETERMINE THE VECTOR TO LEAVE THE BASIS
156C
157 100  K=0
158      DO 110 I=KL,M
159      D=WA(I,IN)
160      IF(D.LE.TOLER)GOTO 110
161      K=K+1
162      WB(K)=WA(I,N1)/D
163      S(K)=I
164      TEST=.TRUE.
165 110  CONTINUE
166 120  IF(K.GT.0)GOTO 130
167      TEST=.FALSE.
168      GOTO 150
169 130  MIN=BIG
170      DO 140 I=1,K
171      IF(WB(I).GE.MIN)GOTO 140
172      J=I
173      MIN=WB(I)
174      OUT=S(I)
175 140  CONTINUE
176      WB(J)=WB(K)
177      S(J)=S(K)
178      K=K-1
179C
180C CHECK FOR LINEAR DEPENDENCE IN STAGE 1
181C
182 150  IF(TEST.OR..NOT.STAGE)GOTO 170
183      DO 160 I=1,M4
184      D=WA(I,KR)
185      WA(I,KR)=WA(I,IN)
186      WA(I,IN)=D
187 160  CONTINUE
188      KR=KR+1
189      GOTO 260
190 170  IF(TEST)GOTO 180
191      WA(M2,N1)=TWO
192      GOTO 390
193 180  PIVOT=WA(OUT,IN)
194      IF(WA(M1,IN)-PIVOT-PIVOT.LE.TOLER)GOTO 200
195      DO 190 J=KR,N1
196      D=WA(OUT,J)
197      WA(M1,J)=WA(M1,J)-D-D
198      WA(M2,J)=WA(M2,J)-D-D
199      WA(OUT,J)=-D
200 190  CONTINUE
201      WA(OUT,N2)=-WA(OUT,N2)
202      GOTO 120
203C
204C PIVOT ON WA(OUT,IN)
205C
206 200  DO 210 J=KR,N1
207      IF(J.EQ.IN)GOTO 210
208      WA(OUT,J)=WA(OUT,J)/PIVOT
209 210  CONTINUE
210      DO 230 I=1,M3
211      IF(I.EQ.OUT)GOTO 230
212      D=WA(I,IN)
213      DO 220 J=KR,N1
214      IF(J.EQ.IN)GOTO 220
215      WA(I,J)=WA(I,J)-D*WA(OUT,J)
216 220  CONTINUE
217 230  CONTINUE
218      DO 240 I=1,M3
219
220      IF(I.EQ.OUT)GOTO 240
221      WA(I,IN)=-WA(I,IN)/PIVOT
222 240  CONTINUE
223      WA(OUT,IN)=ONE/PIVOT
224      D=WA(OUT,N2)
225      WA(OUT,N2)=WA(M4,IN)
226      WA(M4,IN)=D
227      KOUNT=KOUNT+1
228      IF(.NOT.STAGE)GOTO 270
229C
230C INTERCHANGE ROWS IN STAGE 1
231C
232      KL=KL+1
233      DO 250 J=KR,N2
234      D=WA(OUT,J)
235      WA(OUT,J)=WA(KOUNT,J)
236      WA(KOUNT,J)=D
237 250  CONTINUE
238 260  IF(KOUNT+KR.NE.N1)GOTO 70
239C
240C STAGE 2
241C
242 265  STAGE=.FALSE.
243C
244C DETERMINE THE VECTOR TO ENTER THE BASIS
245C
246 270  MAX=-BIG
247      DO 290 J=KR,N
248      D=WA(M1,J)
249      IF(D.GE.ZERO)GOTO 280
250      IF(D.GT.(-TWO))GOTO 290
251      D=-D-TWO
252 280  IF(D.LE.MAX)GOTO 290
253      MAX=D
254      IN=J
255 290  CONTINUE
256      IF(MAX.LE.TOLER)GOTO 310
257      IF(WA(M1,IN).GT.ZERO)GOTO 100
258      DO 300 I=1,M4
259      WA(I,IN)=-WA(I,IN)
260 300  CONTINUE
261      WA(M1,IN)=WA(M1,IN)-TWO
262      WA(M2,IN)=WA(M2,IN)-TWO
263      GOTO 100
264C
265C COMPUTE QUANTILES
266C
267 310  CONTINUE
268      DO 320 I=1,KL-1
269      K=WA(I,N2)*SIGN(ONE,WA(I,N2))
270      X(K) = WA(I,N1) * SIGN(ONE,WA(I,N2))
271 320  CONTINUE
272 390  SUM = ZERO
273      DO 400 I=KL,M
274      K = WA(I,N2) * SIGN(ONE,WA(I,N2))
275      D = WA(I,N1) * SIGN(ONE,WA(I,N2))
276      SUM = SUM + D * SIGN(ONE,D) * (HALF + SIGN(ONE,D)*(T-HALF))
277      K=K-N
278      E(K)=D
279 400  CONTINUE
280      RETURN
281      END
282