Instrumental variable Tobit in R I have a data generating process of the form:
res1 <- rnorm(N);
res2 <- res1*0.5 + rnorm(N)
x    <-  z[,1]*2 + res1;
ys   <- x*b + res2;
d    <- (ys>0); #dummy variable
y    <- d*ys;

(blatantly stolen from an old RHelp thread, where the question remained unanswered back in 2006)
That is, my independent variable x is endogenous but there's an exogenous instrument z, and my dependent variable y is censored below at zero.
Are there any ways I can consistently estimate this model in R? What I'm looking for is essentially an instrumental variable Tobit (as implemented in Stata's ivtobit command).
 A: You ask for two things in the question.  First you ask for a way to consistently estimate this model in R.  Second, you ask for the particular ways implemented in Stata's ivtobit (these ways are full-information maximum likelihood and the Newey, 1987 two-step estimator).  Doing the second would be a nice service for the R community which I don't have time to do.  It seems strange that nobody has done it, though.  However, I can help you with the first.
There are several ways to get a consistent estimator for the model.  Here is an easy and intuitive way, based on the estimator of Rivers and Vuong (J of Econometrics, 1988).  The steps in the estimator are:


*

*Estimate the reduced form for the endogenous variable (i.e. regress x on z in your notation)

*Collect the predicted values and the residuals from that regression

*Use Tobit to estimate the equation of interest, substituting predicted endogenous variables for endogenous variables and including the residuals collected in 2 (i.e. regress y on x-predicted and residuals using Tobit)


The coefficient on the predicted x is the estimator of the coefficient of interest.  
The usual standard errors are wrong --- that is, the standard errors Tobit spits out in step 3 are not the right standard errors.  This is because the Tobit routine does not "know" that the variables you are handing it are fitted variables instead of exact variables.  You could react to this by "doing it right:" looking up the formula for the variance matrix and then implementing it (booooring).  Or you could just bootstrap.  I favor the latter.  It is easier, and if you are really going to "do it right," then you would not use this estimator anyway, since it is not efficient.  Rather, you would actually implement the Newey (1987) estimator or maximum likelihood.
If you have additional right-hand-side variables which are exogenous, just include them along the way in all the commands (both in the tobit and in the reduced form regression).  If you have additional right-hand-side variables which are endogenous, then you treat them like the first endogenous RHS variable.  Run reduced form regressions for each of the endogenous variables.  Include predicted values and residuals in the tobit at the end for each endogenous variable.  As with any proper correction for endogeneity, you have to have as many excluded exogenous variables as you have included endogenous variables for this all to work.
The code below implements both the estimator and the bootstrapped standard errors for the simple example you gave.  I also fixed up the example so that it should work turn-key:
require(censReg)
require(boot)

a    <- 2    # structural parameter of interest
b    <- 1    # strength of instrument
rho  <- 0.5  # degree of endogeneity

N    <- 1000
z    <- rnorm(N)
res1 <- rnorm(N)
res2 <- res1*rho + sqrt(1-rho*rho)*rnorm(N)
x    <- z*b + res1
ys   <- x*a + res2
d    <- (ys>0) #dummy variable
y    <- d*ys

inconsistent.tobit <- censReg(y~x)
summary(inconsistent.tobit)

reduced.form <- lm(x~z)
summary(reduced.form)


consistent.tobit <- censReg(y~fitted(reduced.form)+residuals(reduced.form))
summary(consistent.tobit)


# I'd like bootstrapped standard errors, please!
my.data <- data.frame(y,x,z)
tobit_2siv_coef <- function(data,indices){
  d <- data[indices,]
  reduced.form <- lm(x~z,data=d)
  consistent.tobit <- censReg(d[,"y"]~fitted(reduced.form)+residuals(reduced.form))
  return(summary(consistent.tobit)$estimate["fitted(reduced.form)",1])
}

boot.results <- boot(data=my.data,statistic=tobit_2siv_coef,R=100)
boot.results

