# This is a special version of K-dO rq modified to compute the LR process # Note that the sol array is now p+3 by J and the third row contains the # value of the objective function at the tau specified in row one. # # This software is in the public domain and may be freely # used and redistributed for non-commercial purposes. No guarantees # are offered or implied. comments, bug reports, etc are welcome # and should be sent to roger@ysidro.econ.uiuc.edu or to # # roger koenker # department of economics # university of illinois # champaign, illinois, 61820 # # subroutine rqbr(m,nn,m5,n3,n4,a,b,t,toler,ift,x,e,s,wa,wb,nsol,ndsol,sol,dsol,lsol,h,qn,cutoff,ci,tnmat,big,lci1) # # m = number of observations # n = number of parameters # m5 = m+5 # n3 = n+3 # n4 = n+4 # a is the x matrix # b is the y vector # t, the desired quantile # toler, smallest detectable |x-y|/x machine precision to the 2/3 # ift exit code: # 0-ok # else dimensions inconsistent see below # x the parameter estimate betahat # e is the residual vector # s is an integer work array # wa is a real work array # wb is another real work array # nsol is an estimated (row) dimension of the primal solution array say 3*m; minimum value = 2 # ndsol is an estimated (row) dimension of the dual solution array say 3*m; minimum value = 2 # sol is the primal solution array # dsol is the dual solution arry # lsol is the actual dimension of the solution arrays # h is the matrix of basic observations indices # qn is the vector of residuals variances from the projection of each column of the x matrix on the remaining columns # cutoff, the critical point for N(0,1) # ci is the matrix of confidence intervals # tnmat is the matrix of the JGPK rank test statistics # big, large positive finite floating-point number # utilization: if you just want a solution at a single quantile you # neednt bother with sol, nsol, etc, if you want all the solutions # then set theta to something <0 and sol and dsol will return all the # estimated quantile solutions. # the algorithm is a slightly modified version of algorithm as 229 # described in koenker and dorey, computing regression quantiles, # applied statistics, pp. 383-393. # integer i,j,k,kl,kount,kr,l,lsol,m,m1,m2,m3,m4,m5,ift integer n,n1,n2,n3,nsol,ndsol,out,s(m),h(nn,nsol) integer nn,n4,idxcf logical stage,test,init,iend,lup logical lci1,lci2,skip double precision a1,aux,b1,big,d,dif,pivot,smax,t,t0,t1,tnt double precision min,max,toler,zero,half,one,two double precision b(m),sol(n3,nsol),a(m,nn),x(nn),wa(m5,n4),wb(m) double precision sum,e(m),dsol(m,ndsol) double precision qn(nn),cutoff,ci(4,nn),tnmat(4,nn),tnew,told,tn parameter( zero = 0.d0) parameter( one = 1.d0) parameter( two = 2.d0) # # check dimension parameters # n=nn ift = 0 wa(m+2,nn+1) = one if (m5!=m+5) ift = 3 if (n3!=n+3) ift = 4 if (n4!=n+4) ift = 5 if (m<=zero||n<=zero) ift = 6 if (ift<=two) { # # initialization # half = one/two iend = .true. lci2 = .false. lup = .true. skip = .false. idxcf = 0 tnew = zero tn = zero m1 = m+1 n1 = n+1 n2 = n+2 m2 = m+2 m3 = m+3 m4 = m+4 do j = 1,n{ x(j) = zero } do i = 1,m e(i) = zero if (tone) { t0 = one/(two*float(m)) t1 = one-t0 t = t0 iend = .false. lci1 = .false. } repeat { #0 do i = 1,m { k = 1 do j = 1,nn if (k<=nn){ if(j==idxcf) skip = .true. else skip = .false. if(!skip){ wa(i,k) = a(i,j) k = k+1 } } wa(i,n4) = n+i wa(i,n2) = b(i) if (idxcf != 0) wa(i,n3) = tnew*a(i,idxcf) else wa(i,n3) = zero wa(i,n1) = wa(i,n2)-wa(i,n3) if (wa(i,n1)(-two)) next 1 d = -d-two } if (d>max) { max = d in = j } } if (max<=toler) break 2 if (wa(m1,in)<=zero) { do i = 1,m4 wa(i,in) = -wa(i,in) wa(m1,in) = wa(m1,in)-two wa(m2,in) = wa(m2,in)-two } repeat { #4 # # determine the vector to leave the basis # k = 0 do i = kl,m { d = wa(i,in) if (d>toler) { k = k+1 wb(k) = wa(i,n1)/d s(k) = i test = .true. } } repeat { #5 if (k<=0) test = .false. else { min = big do i = 1,k if (wb(i)max) { max = d in = j } } if (wa(m1,in)zero) dsol(kd,lsol) = one } if (!lci2){ sol(1,lsol) = smax sol(2,lsol) = sum sum = zero do j=kl,m{ d = wa(j,n1)*sign(one,wa(j,n4)) sum = sum + d*(smax + half*(sign(one,d) - one)) } sol(3,lsol) = sum do i=1,m dsol(i,lsol+1) = dsol(i,lsol) } if (lci2){ # compute next theta a1 = zero do i = 1,m { a1 = a1+a(i,idxcf)*(dsol(i,lsol)+t-one) } tn = a1/sqrt(qn(idxcf)*t*(one-t)) if (abs(tn)tnew) if (tntsmax){ smax = tnt out = i } } } if (lup){ told = tnew tnew = smax + toler ci(3,idxcf) = told - toler tnmat(3,idxcf) = tn if (!(tnew < big-toler)){ ci(3,idxcf) = big ci(4,idxcf) = big tnmat(3,idxcf) = tn tnmat(4,idxcf) = tn lup = .false. go to 70 } } else{ told = tnew tnew = smax - toler ci(2,idxcf) = told + toler tnmat(2,idxcf) = tn if (!(tnew > -big+toler)){ ci(2,idxcf) = -big ci(1,idxcf) = -big tnmat(2,idxcf) = tn tnmat(1,idxcf) = tn lup = .true. go to 60 } } #update the new marginal cost do i = 1,m{ wa(i,n3) = wa(i,n3)/told*tnew wa(i,n1) = wa(i,n2) - wa(i,n3) } do j = kr,n3{ d = wa(out,j) wa(m1,j) = wa(m1,j) -d -d wa(m2,j) = wa(m2,j) -d -d wa(out,j) = -d } wa(out,n4) = -wa(out,n4) init = .true. } else{ if (lup){ ci(4,idxcf) = tnew - toler tnmat(4,idxcf) = tn lup = .false. go to 70 } else{ ci(1,idxcf) = tnew + toler tnmat(1,idxcf) = tn lup = .true. go to 60 } } } if ((iend)&&(!lci2)) go to 40 if (!lci2){ init = .true. lsol = lsol+1 do i = 1,m s(i) = zero do j = 1,n x(j) = zero # # compute next t # smax = two do j = 1,n { b1 = wa(m3,j) a1 = (-two-wa(m2,j))/b1 b1 = -wa(m2,j)/b1 if (a1>=t) if (a1t) if (b1=t1+toler) iend = .true. t = tnt if (iend) t = t1 } } #1 wa(m2,nn+1) = two ift = 2 go to 50 40 if (lsol>2) { sol(1,1) = zero #sol(2,1) = zero sol(3,1) = zero sol(1,lsol) = one #sol(2,lsol) = zero sol(3,lsol) = zero do i = 1,m { dsol(i,1) = one dsol(i,lsol) = zero dsol(i,lsol+1) = zero } } l = kl-1 do i = 1,l if (wa(i,n1)nn){ break 1 } 70 if (lup){ tnew = x(idxcf)+toler told = tnew ci(3,idxcf) = x(idxcf) tnmat(3,idxcf) = zero } else{ tnew = x(idxcf)-toler told = tnew ci(2,idxcf) = x(idxcf) tnmat(2,idxcf) = zero } } } #0 # restore the original value of dual when ci is true do i=1,m dsol(i,lsol) = dsol(i,lsol+1) } return end