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