# mulligan rubinstein bounds / simSQB.R

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72``` ```## A brief script illustrating the use of selectionQuantiles.R ## Simulates a sample selection problem and then estimates the set of ## best linear approximations to the conditional quantile function. rm(list=ls()) ## Settings N <- 1000 ## number of observations nx <- 2 ## number of covariates (1st one is an intercept) nv <- 3 ## number of points of support for exclusion rho <- 0.2 ## correlation between selection and outcome equation shocks (they are jointly normal) sig <- 1 ## standard deviation of shock to outcome b <- matrix(rep(1,nx),nrow=nx,ncol=1) ## coefficients for outcome g <- matrix(rep(1,nx+1),nrow=nx+1,ncol=1) ## coefficients for selection g[nx+1,1] <- 1.5 ## Simulation x <- cbind(rep(1,N),matrix(runif(N*(nx-1)),nrow=N,ncol=nx-1)) if (nx>2) { x[,3] <- round(x[,3]*3)/3 } v <- (t(rmultinom(N,1,rep(1/nv,nv))) %*% as.matrix(seq(1,nv)) - 1 - (nv-1)/2)/(nv-1) e <- rnorm(N) u <- (rho*e + sqrt(1-rho^2)*rnorm(N)) g[1] <- 0.3 d <- cbind(x,v)%*%g + u ystar <- x%*%b + e*sig y <- ystar y[d<0] <- NaN y0 <- y y0[d<0] <- -20 y1 <- y y1[d<0] <- 20 x <- as.matrix(x[,2:nx]) vf <- as.factor(1-v) ## Estimation require(quantreg) require(ggplot2) summary(rq(y ~ x)) if (nx>2) { fmla <- as.formula(paste("y~poly(", paste( paste("X",1:(nx-1),sep="",collapse=",") ,",degree=2)*vf"))) } else { fmla <- y~poly(x,degree=4)*vf } data <- data.frame(x,y0,y1,vf,y,d) ## Plot conditional quantile functions at #if (nx==2) { # ggplot(data, aes(x=x)) + geom_point(aes(y=rq(update(fmla,y1~.),data=data)\$fit)) + # geom_point(aes(y=rq(update(fmla,y0~.),data=data)\$fit),colour="red") #} else { # ggplot(data, aes(x=x[,1])) + geom_point(aes(y=rq(update(fmla,y1~.),data=data)\$fit)) + # geom_point(aes(y=rq(update(fmla,y0~.),data=data)\$fit),colour="red") #} source("selectionQuantiles.R",echo=FALSE) bnds <- selectionQuantileBounds(y~x,fmla,exclusion.levels=levels(vf),exclusion.name="vf", tau=c(0.2,0.4,0.5,0.6,0.8),y0=y0,y1=y1, data=data) #setplot <- plot(bnds,q0=matrix(rep(0,nx)),dims=c(nx-1,nx),lo=0.5,hi=1.5,nq=25) #setplot ## bootstrap boot.sfn <- selectionQuantileBounds(y~x,fmla,exclusion.levels=levels(vf),exclusion.name="vf", tau=c(0.2,0.4,0.5,0.6,0.8),y0=y0,y1=y1, data=data, bootstrap=TRUE, boot.rep=10, boot.save.interval=100000, boot.parallel=1, boot.intersect.method="CLR") ```