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