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