1

i am working on a two stage threshold quantile regression model in R, and my aim is to estimate the threshold of the reduced-form equation (call it rhohat), and the threshold of the structural equation (call it qhat), in two stages. On the first stage, i estimate rhohat by quantile regression and obtain the fitted values. I use these fitted values to estimate qhat on the second stage. The code is as follows (thanks to Prof. Bruce Hansen, whose code i modified): 

#************************************************************#
#Quantile Regression.
qr.regress <- function(y,x){
beta <- c(rq(y~x,tau)$coefficients[1],rq(y~x,tau)$coefficients[2])
beta
}

#Threshold Estimation with one independent variable + constant.
joint_thresh <- function(y,x,q){
n=nrow(y)
k=ncol(x)
e=y-x%*%rq(y~x,tau)$coefficients[2]-rq(y~x,tau)$coefficients[1]
s0 <- det(t(e)%*%e)    
n1 <- round(.05*n)+k
n2 <- round(.95*n)-k
qs <- sort(q)
qs <- qs[n1:n2]
qs <- as.matrix(unique(qs))
qn <- nrow(qs)
sn <- matrix(0,qn,1)
for (r in 1:qn){
   d <- (q<=qs[r])
  xx <- (x)*(d%*%matrix(1,1,k))
  xx <- xx-x%*%rq(xx~x,tau)$coefficients[2]-rq(xx~x,tau)$coefficients[1]
  ex <- e-xx%*%rq(e~xx,tau)$coefficients[2]-rq(e~xx,tau)$coefficients[1]
  exw <- ex*(tau-(ex<0))
  sn[r] <- sum(exw)   
}
r <- which.min(sn)
smin <- sn[r] 
qhat <- qs[r]
d <- (q<=qhat)
x1 <- x*(d%*%matrix(1,1,k))
x2 <- x*((1-d)%*%matrix(1,1,k))
beta1 <- rq(y~x1,tau)$coefficients[2]
beta2 <- rq(y~x2,tau)$coefficients[2]
yhat <- x1%*%beta1+x2%*%beta2
list(yhat=yhat,qhat=qhat)
}

#Threshold Estimation with two independent variables + constant.
joint_thresh_2 <- function(y,x,q){
n <- nrow(y)
k <- ncol(x)
e=y-x[,1]%*%t(rq(y~x-1,tau)$coefficients[3])-x[,2]%*%t(rq(y~x-1,tau)$coefficients[2])-x[,3]%*%t(rq(y~x-1,tau)$coefficients[3])
s0 <- det(t(e)%*%e)    
n1 <- round(.05*n)+k
n2 <- round(.95*n)-k
qs <- sort(q)
qs <- qs[n1:n2]
qs <- as.matrix(unique(qs))
qn <- nrow(qs)
sn <- matrix(0,qn,1)
for (r in 1:qn){
  d <- (q<=qs[r])
  xx <- (x)*(d%*%matrix(1,1,k))
  xx <- xx[,1]-x%*%rq(x[,1]~xx-1,tau)$coefficients
  ex <- e-xx%*%qr.regress(e,xx)[2]-xx%*%qr.regress(e,xx)[1]
  exw <- ex*(tau-(ex<0))
  sn[r] <- sum(exw)  
}
r <- which.min(sn)
smin <- sn[r]
qhat <- qs[r]
d <- (q<=qhat)
x1 <- x*(d%*%matrix(1,1,k))
x2 <- x*((1-d)%*%matrix(1,1,k))
beta1 <- rq(y~x1-1,tau)$coefficients[1]
beta2 <- rq(y~x2-1,tau)$coefficients[3]
c1 <- rq(y~x1-1,tau)$coefficients[2]
c2 <- rq(y~x2-1,tau)$coefficients[2]
yhat <- x1[,1]%*%t(beta1)+x2[,3]%*%t(beta2)+c1+c2
list(yhat=yhat,qhat=qhat)
}

#Threshold Reduced-form eqn. 
tau=0.50

stqr_thresh_loop <- function(n,reps){
qhat=as.vector(reps)
rhohat=as.vector(reps)
kx <- 1 
sig <- matrix(c(1,0.5,0.5,1),2,2)  
x<- matrix(rnorm(n*kx),n,kx)
q <- matrix(rnorm(n),n,1) 
z2 <- cbind(matrix(1,n,1),q)
for(i in 1:reps){
e <- matrix(rnorm(n*2,quantile(rnorm(n),tau),1),n,2)%*%chol(sig)
z1=0.5+0.5*x*(q<=0)+1*x*(q>0)+e[,2] 
y=0.5+1*z1*(q<=1)+1.5*z1*(q>1)+1*z2[,2]+e[,1]
    out1 <-  joint_thresh(y=z1,x=x,q=q)
    z1hat<- out1$yhat
    rhohat[i] <- out1$qhat
zhat <- cbind(z1hat,z2)

out2 <- joint_thresh_2(y=y,x=zhat,q=q)
qhat[i] <- out2$qhat
        } #Close for loop.
list(rhohat=rhohat,qhat=qhat)
}

#************************************************************#

You can easily run it yourselves. The problem is that when i write,

stqr_thresh_loop(n=200,reps=500)

the code crashes and never gives me any result. What am i doing wrong? Thank you very much!!

Pafjo
  • 4,979
  • 3
  • 23
  • 31
user1531131
  • 11
  • 1
  • 2
    How does it crash, any error messages? – Paul Hiemstra Jul 17 '12 at 08:51
  • Unfortunately, there is no message at all! Nothing! Just "the program does not respond". It is bizzare! Thanks!! – user1531131 Jul 17 '12 at 08:55
  • 1
    It probably means that your program is stuck in an infinite loop, or the number of iterations just becomes very large with `n = 200` and `reps = 500`. – Paul Hiemstra Jul 17 '12 at 08:56
  • 1
    Here are some debugging tools to get you started: http://stackoverflow.com/questions/4442518/general-suggestions-for-debugging-r – Ari B. Friedman Jul 17 '12 at 08:58
  • If it has stuck on an infinite loop, how to correct it? Thanks again!! – user1531131 Jul 17 '12 at 09:04
  • 1
    I'd first take a look at the values of the parameters over which you loop (e.g. `qn`) to see if any have absurdly large values. You could also include a bunch of print statements to narrow down which part of the program (more specifically which loop) causes the problem. – Paul Hiemstra Jul 17 '12 at 09:18
  • Where does the 'rq' function come from? – Spacedman Jul 17 '12 at 11:22
  • This is from the "quantreg" library of Roger Koenker! – user1531131 Jul 17 '12 at 11:31
  • Would be an idea to stick require(quantreg) somewhere then! – Spacedman Jul 17 '12 at 14:27
  • All right: you need the following packages: "quantreg", "SparseM" and "scatterplot3d". Or only "quantreg" is enough. Did you run it yourself? – user1531131 Jul 17 '12 at 15:29
  • I can now! It might just be taking a very long time. Some tests with small n and reps indicate the time taken is linear with n and reps, which makes me believe n=200, reps=500 will take 100x as long as n=20,reps=50, which takes 12s on my system. Have you tried waiting 20 minutes? – Spacedman Jul 17 '12 at 15:40
  • First thanks a lot for the effort! I do not mind waiting. I have done so in the past, running longer codes. The problem is that my R crashes. Vanishes!! It says "this program does not respond" and i need to do control-alt-delete to exit. This is my problem. Can you help me please? – user1531131 Jul 17 '12 at 15:44
  • So the real question is 'how do I stop a running R script?'? What operating system? What R GUI? Sounds like Windows. Are you hitting some "Stop" button? Its possible that inside C code R will be unresponsive because its not handling signals. – Spacedman Jul 17 '12 at 15:46
  • Yes the R code is unresponsive! And i do not think i am stuck on some infinte loop since i checked that and it is working fine. I do not hit stop buttons. The code crashes itself. – user1531131 Jul 17 '12 at 16:02
  • How do you know its unresponsive if you aren't hitting buttons or doing anything? What's it not responding to? – Spacedman Jul 17 '12 at 18:11
  • Sorry for the delay as i am living it Taiwan and there is time lag. The code hangs due to some problem with the rq routine. Maybe sth to do with the data. – user1531131 Jul 18 '12 at 06:19

1 Answers1

2

Try some timing with smaller values:

> system.time({z = stqr_thresh_loop(n=10,reps=5)})
   user  system elapsed 
  0.612   0.000   0.615 
> system.time({z = stqr_thresh_loop(n=20,reps=5)})
   user  system elapsed 
  1.248   0.000   1.249 
> system.time({z = stqr_thresh_loop(n=40,reps=5)})
   user  system elapsed 
  2.740   0.000   2.743 
> system.time({z = stqr_thresh_loop(n=10,reps=10)})
   user  system elapsed 
  1.228   0.000   1.233 
> system.time({z = stqr_thresh_loop(n=10,reps=20)})
   user  system elapsed 
  2.465   0.000   2.472 
> system.time({z = stqr_thresh_loop(n=10,reps=40)})
   user  system elapsed 
  4.968   0.000   4.969

This looks nicely linear in elapsed time with n and reps. So the time for n=200,reps=500 will be at least 100x that for n=20,reps=50, which is (on my system) going to be about 20 minutes. Have you waited that long?

Try printing the value of 'i' in your loop from 1:reps to get a handle on progress.

Spacedman
  • 92,590
  • 12
  • 140
  • 224