https://github.com/cran/quantreg
Tip revision: d37f852036a87c4097be1667f8d54418fc665972 authored by Unknown author on 08 August 1977, 00:00:00 UTC
version 2.0-2
version 2.0-2
Tip revision: d37f852
rq.f
C Output from Public domain Ratfor, version 1.0
C (Version from Brian Ripley's collection SWin, on his WEB site).
C In addition all real declarations changed to DOUBLE PRECISION.
C Kjetil Halvorsen
C
C Below is original text:
C this software has been submitted to statlib and may be freely
C used and redistributed for non-commercial purposes. no guarantees
C are offered or implied. comments, bug reports, etc are welcome
C and should be sent to roger@ysidro.econ.uiuc.edu or to
C
C roger koenker
C department of economics
C university of illinois
C champaign, illinois, 61820
C
C
subroutine rq(m,nn,m5,n2,n4,a,b,t,toler,ift,x,e,s,wa,wb,nsol,ndsol
*,sol,dsol,lsol,h,qn,cutoff,ci,tnmat,big,lci1)
C
C number of observations
C number of parameters
C m+5
C n+2
C n+4
C a is the x matrix
C b is the y vector
C t, the desired quantile
C toler, smallest detectable |x-y|/x machine precision to the 2/3
C ift exit code:
C 0-ok
C else dimensions inconsistent
C x the parameter estimate betahat
C e is the residual vector
C s is an integer work array
C wa is a real work array
C wb is another real work array
C nsol is an estimated (row) dimension of the primal solution array say 3*m;
C minimum value = 2
C ndsol is an estimated (row) dimension of the dual solution array say 3*m;
C minimum value = 2
C sol is the primal solution array
C dsol is the dual solution arry
C lsol is the actual dimension of the solution arrays
C h is the matrix of basic observations indices
C qn is the vector of residuals variances from the projection of each column
C of the x matrix on the remaining columns
C cutoff, the critical point for N(0,1)
C ci is the matrix of confident intervals
C tnmat is the matrix of the JGPK rank test statistics
C big, large positive finite floating-point number
C utilization: if you just want a solution at a single quantile you
C needn't bother with sol, nsol, etc, if you want all the solutions
C then set theta to something <0 and sol and dsol will return all the
C estimated quantile solutions.
C the algorithm is a slightly modified version of algorithm as 229
C described in koenker and d'orey, "computing regression quantiles,
C applied statistics, pp. 383-393.
C
integer i,j,k,kl,kount,kr,l,lsol,m,m1,m2,m3,m4,m5,ift
integer n,n1,n2,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(n2,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
data zero/0.00/
data half/0.50/
data one/1.00/
data two/2.00/
C
C check dimension parameters
C
n=nn
ift = 0
wa(m+2,nn+1) = one
if(m5.ne.m+5)then
ift = 3
endif
if(n2.ne.n+2)then
ift = 4
endif
if(n4.ne.n+4)then
ift = 5
endif
if(m.le.zero.or.n.le.zero)then
ift = 6
endif
if(ift.le.two)then
C
C initialization
C
iend = .true.
lci2 = .false.
lup = .true.
skip = .false.
idxcf = 0
tnew = zero
tn = zero
m1 = m+1
n1 = n+1
n3 = n+3
m2 = m+2
m3 = m+3
m4 = m+4
do23010 j = 1,n
x(j) = zero
23010 continue
23011 continue
do23012 i = 1,m
e(i) = zero
23012 continue
23013 continue
if(t.lt.zero.or.t.gt.one)then
t0 = one/float(m)-toler
t1 = one-t0
t = t0
iend = .false.
lci1 = .false.
endif
23016 continue
C0
do23019 i = 1,m
k = 1
do23021 j = 1,nn
if(k.le.nn)then
if(j.eq.idxcf)then
skip = .true.
else
skip = .false.
endif
if(.not.skip)then
wa(i,k) = a(i,j)
k = k+1
endif
endif
23021 continue
23022 continue
wa(i,n4) = n+i
wa(i,n2) = b(i)
if(idxcf .ne. 0)then
wa(i,n3) = tnew*a(i,idxcf)
else
wa(i,n3) = zero
endif
wa(i,n1) = wa(i,n2)-wa(i,n3)
if(wa(i,n1).lt.zero)then
do23033 j = 1,n4
wa(i,j) = -wa(i,j)
23033 continue
23034 continue
endif
23019 continue
23020 continue
do23035 j = 1,n
wa(m4,j) = j
wa(m2,j) = zero
wa(m3,j) = zero
do23037 i = 1,m
aux = sign(one,wa(m4,j))*wa(i,j)
wa(m2,j) = wa(m2,j)+aux*(one-sign(one,wa(i,n4)))
wa(m3,j) = wa(m3,j)+aux*sign(one,wa(i,n4))
23037 continue
23038 continue
wa(m3,j) = two*wa(m3,j)
23035 continue
23036 continue
dif = zero
init = .false.
if(.not.lci2)then
Ccompute the columns means
do23041 k = 1,n
wa(m5,k) = zero
do23043 i = 1,m
wa(m5,k) = wa(m5,k)+a(i,k)
23043 continue
23044 continue
wa(m5,k) = wa(m5,k)/float(m)
23041 continue
23042 continue
endif
lsol = 1
kount = 0
23045 continue
C1
C
C compute new marginal costs
C
do23048 j = 1,n
wa(m1,j) = wa(m2,j)+wa(m3,j)*t
23048 continue
23049 continue
if(.not.init)then
C
C stage 1
C
C determine the vector to enter the basis
C
stage = .true.
kr = 1
kl = 1
go to 30
endif
23052 continue
C2
C
C stage 2
C
stage = .false.
23055 continue
C3
C
C determine the vector to enter the basis
C
max = -big
do23058 j = kr,n
d = wa(m1,j)
if(d.lt.zero)then
if(d.gt.(-two))then
goto 23058
endif
d = -d-two
endif
if(d.gt.max)then
max = d
in = j
endif
23058 continue
23059 continue
if(max.le.toler)then
goto 23054
endif
if(wa(m1,in).le.zero)then
do23070 i = 1,m4
wa(i,in) = -wa(i,in)
23070 continue
23071 continue
wa(m1,in) = wa(m1,in)-two
wa(m2,in) = wa(m2,in)-two
endif
23072 continue
C4
C
C determine the vector to leave the basis
C
k = 0
do23075 i = kl,m
d = wa(i,in)
if(d.gt.toler)then
k = k+1
wb(k) = wa(i,n1)/d
s(k) = i
test = .true.
endif
23075 continue
23076 continue
23079 continue
C5
if(k.le.0)then
test = .false.
else
min = big
do23084 i = 1,k
if(wb(i).lt.min)then
j = i
min = wb(i)
out = s(i)
endif
23084 continue
23085 continue
wb(j) = wb(k)
s(j) = s(k)
k = k-1
C
C check for linear dependence in stage 1
C
endif
if(.not.test.and.stage)then
goto 23081
endif
if(.not.test)then
goto 23047
endif
pivot = wa(out,in)
if(wa(m1,in)-pivot-pivot.le.toler)then
go to 10
endif
do23094 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
23094 continue
23095 continue
wa(out,n4) = -wa(out,n4)
C5
23080 goto 23079
23081 continue
do23096 i = 1,m4
d = wa(i,kr)
wa(i,kr) = wa(i,in)
wa(i,in) = d
23096 continue
23097 continue
kr = kr+1
go to 20
C
C pivot on wa(out,in)
C
10 do23098 j = kr,n3
if(j.ne.in)then
wa(out,j) = wa(out,j)/pivot
endif
23098 continue
23099 continue
do23102 i = 1,m3
if(i.ne.out)then
d = wa(i,in)
do23106 j = kr,n3
if(j.ne.in)then
wa(i,j) = wa(i,j)-d*wa(out,j)
endif
23106 continue
23107 continue
endif
23102 continue
23103 continue
do23110 i = 1,m3
if(i.ne.out)then
wa(i,in) = -wa(i,in)/pivot
endif
23110 continue
23111 continue
wa(out,in) = one/pivot
d = wa(out,n4)
wa(out,n4) = wa(m4,in)
wa(m4,in) = d
kount = kount+1
if(.not.stage)then
goto 23074
C
C interchange rows in stage 1
C
endif
kl = kl+1
do23116 j = kr,n4
d = wa(out,j)
wa(out,j) = wa(kount,j)
wa(kount,j) = d
23116 continue
23117 continue
20 if(kount+kr.eq.n1)then
goto 23057
endif
30 max = -one
do23120 j = kr,n
if(abs(wa(m4,j)).le.n)then
d = abs(wa(m1,j))
if(d.gt.max)then
max = d
in = j
endif
endif
23120 continue
23121 continue
if(wa(m1,in).lt.zero)then
do23128 i = 1,m4
wa(i,in) = -wa(i,in)
23128 continue
23129 continue
endif
C4
23073 goto 23072
23074 continue
C3
23056 goto 23055
23057 continue
C2
23053 goto 23052
23054 continue
if(kr.eq.1)then
do23132 j = 1,n
d = abs(wa(m1,j))
if(d.le.toler.or.two-d.le.toler)then
ift = 1
wa(m2,nn+1) = zero
go to 80
endif
23132 continue
23133 continue
endif
80 kount = 0
sum = zero
if(.not.lci2)then
do23138 i = 1,kl-1
k = wa(i,n4)*sign(one,wa(i,n4))
x(k) = wa(i,n1)*sign(one,wa(i,n4))
23138 continue
23139 continue
endif
do23140 i = 1,n
kd = abs(wa(m4,i))-n
dsol(kd,lsol) = one+wa(m1,i)/two
if(wa(m4,i).lt.zero)then
dsol(kd,lsol) = one-dsol(kd,lsol)
endif
if(.not.lci2)then
sum = sum+x(i)*wa(m5,i)
sol(i+2,lsol) = x(i)
h(i,lsol) = kd
endif
23140 continue
23141 continue
do23146 i = kl,m
kd = abs(wa(i,n4))-n
if(wa(i,n4).lt.zero)then
dsol(kd,lsol) = zero
endif
if(wa(i,n4).gt.zero)then
dsol(kd,lsol) = one
endif
23146 continue
23147 continue
if(.not.lci2)then
sol(1,lsol) = smax
sol(2,lsol) = sum
do23154 i=1,m
dsol(i,lsol+1) = dsol(i,lsol)
23154 continue
23155 continue
endif
if(lci2)then
C compute next theta
a1 = zero
do23158 i = 1,m
a1 = a1+a(i,idxcf)*(dsol(i,lsol)+t-one)
23158 continue
23159 continue
tn = a1/sqrt(qn(idxcf)*t*(one-t))
if(abs(tn).lt.cutoff)then
if(lup)then
smax = big
else
smax = -big
endif
do23164 i =1,kl-1
k = wa(i,n4)*sign(one,wa(i,n4))
sol(k,1) = wa(i,n2)*sign(one,wa(i,n4))
sol(k,2) = wa(i,n3)*sign(one,wa(i,n4))/tnew
23164 continue
23165 continue
do23166 i = kl,m
a1 = zero
b1 = zero
k = wa(i,n4)*sign(one,wa(i,n4))-n
l = 1
do23168 j = 1,n
if(j.eq.idxcf)then
l = l+1
endif
a1 = a1 + a(k,l)*sol(j,1)
b1 = b1 + a(k,l)*sol(j,2)
l = l+1
23168 continue
23169 continue
tnt = (b(k)-a1)/(a(k,idxcf)-b1)
if(lup)then
if(tnt.gt.tnew)then
if(tnt.lt.smax)then
smax = tnt
out = i
endif
endif
else
if(tnt.lt.tnew)then
if(tnt.gt.smax)then
smax = tnt
out = i
endif
endif
endif
23166 continue
23167 continue
if(lup)then
told = tnew
tnew = smax + toler
ci(3,idxcf) = told - toler
tnmat(3,idxcf) = tn
if(.not.(tnew .lt. big-toler))then
ci(3,idxcf) = big
ci(4,idxcf) = big
tnmat(3,idxcf) = tn
tnmat(4,idxcf) = tn
lup = .false.
go to 70
endif
else
told = tnew
tnew = smax - toler
ci(2,idxcf) = told + toler
tnmat(2,idxcf) = tn
if(.not.(tnew .gt. -big+toler))then
ci(2,idxcf) = -big
ci(1,idxcf) = -big
tnmat(2,idxcf) = tn
tnmat(1,idxcf) = tn
lup = .true.
go to 60
endif
Cupdate the new marginal cost
endif
do23188 i = 1,m
wa(i,n3) = wa(i,n3)/told*tnew
wa(i,n1) = wa(i,n2) - wa(i,n3)
23188 continue
23189 continue
do23190 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
23190 continue
23191 continue
wa(out,n4) = -wa(out,n4)
init = .true.
else
if(lup)then
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
endif
endif
endif
if((iend).and.(.not.lci2))then
go to 40
endif
if(.not.lci2)then
init = .true.
lsol = lsol+1
do23198 i = 1,m
s(i) = zero
23198 continue
23199 continue
do23200 j = 1,n
x(j) = zero
C
C compute next t
C
23200 continue
23201 continue
smax = two
do23202 j = 1,n
b1 = wa(m3,j)
a1 = (-two-wa(m2,j))/b1
b1 = -wa(m2,j)/b1
if(a1.ge.t)then
if(a1.lt.smax)then
smax = a1
dif = (b1-a1)/two
endif
endif
if(b1.gt.t)then
if(b1.lt.smax)then
smax = b1
dif = (b1-a1)/two
endif
endif
23202 continue
23203 continue
tnt = smax+toler*(one+abs(dif))
if(tnt.ge.t1+toler)then
iend = .true.
endif
t = tnt
if(iend)then
t = t1
endif
endif
C1
23046 goto 23045
23047 continue
wa(m2,nn+1) = two
ift = 2
go to 50
40 if(lsol.gt.2)then
sol(1,1) = zero
sol(1,lsol) = one
do23218 i = 1,m
dsol(i,1) = one
dsol(i,lsol) = zero
dsol(i,lsol+1) = zero
23218 continue
23219 continue
endif
l = kl-1
do23220 i = 1,l
if(wa(i,n1).lt.zero)then
do23224 j = kr,n4
wa(i,j) = -wa(i,j)
23224 continue
23225 continue
endif
23220 continue
23221 continue
50 sum = zero
if(.not.lci2)then
do23228 i = kl,m
k = wa(i,n4)*sign(one,wa(i,n4))
d = wa(i,n1)*sign(one,wa(i,n4))
sum = sum+d*sign(one,d)*(half+sign(one,d)*(t-half))
k = k-n
e(k) = d
23228 continue
23229 continue
wa(m2,n2) = kount
wa(m1,n2) = n1-kr
wa(m1,n1) = sum
endif
if(wa(m2,nn+1).eq.two)then
goto 23018
endif
if(.not.lci1)then
goto 23018
endif
if(.not.lci2)then
lci2 = .true.
n = nn-1
n1 = n+1
n2 = n+2
n3 = n+3
n4 = n+4
60 idxcf = idxcf+1
if(idxcf.gt.nn)then
goto 23018
endif
70 if(lup)then
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
endif
endif
C0
C restore the original value of dual when ci is true
23017 goto 23016
23018 continue
do23240 i=1,m
dsol(i,lsol) = dsol(i,lsol+1)
23240 continue
23241 continue
endif
return
end