R阈值分位数回归代码R

Question

我正在研究R中的两阶段阈值分位数回归模型，我的目标是估计简化形式方程（称之为rhohat）的阈值，以及结构的阈值方程式（称之为qhat），分两个阶段。在第一阶段，我通过分位数回归估计rhohat并获得拟合值。我用这些拟合的值来估计第二阶段的qhat。代码如下（感谢布鲁斯·汉森教授，他的代码，我修改）：

#************************************************************# 
#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) 
} 

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

您可以轻松地运行它自己。的问题是，当我写，

stqr_thresh_loop（N = 200，代表= 500）

代码崩溃和从未提出的任何结果。 我在做什么错？ 非常感谢！

Answer 1

尝试一些定时具有较小值：

> 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 

这看起来很好地与n和代表经过的时间成线​​性关系。所以n = 200，reps = 500的时间至少为100x，对于n = 20，reps = 50，这是（在我的系统上）大概是20分钟。你等了那么久吗？

尝试从1：代表在循环中打印“我”的值，以获得进度的句柄。