Honore-style fixed effects estimators for censored data in R

Below is some R code to implement fixed effects Tobit, a la Honore (1992) "Trimmed LAD and Least Squares Estimation of Truncated and Censored Regression Models with Fixed Effects", and also fixed effects double-censored Tobit a la Alan, Honore and Petersen (2008).

# FE tobit a la Honore (1992). See Honore (2002) Port Econ J for a simple guide

# y1 partner y2 other

honore <- function (b, x1, x2, y1, y2) {

dxb <- (x1 - x2) %*% b

sum(

(pmax(y1, dxb) - pmax(y2, - dxb) - dxb)^2 +

2*(y1 < dxb)*(dxb-y1)*y2 +

2*(y1 < -dxb)* (-dxb-y2)*y1

)

}

fetobit <- function (x1, x2, y1, y2) {

# x <- model.matrix(form, dataset)

initvals <- lm.fit(x1 - x2, y1 - y2)$coefficients

res <- optim(initvals, fn=honore, x1=x1, x2=x2, y1=y1, y2=y2,

method="BFGS", control=list(reltol=1e-8, maxit=1000))

if (res$convergence != 0) warning("Didn't converge")

res$par

}

# estimate the standard error, asymptotically.

fetstderr <- function(x1, x2, y1, y2, bhat) {

N <- nrow(x1)

x <- x1 - x2

xb <- x %*% bhat

G <- matrix(0, nrow=length(bhat), ncol=length(bhat))

V <- G

# this should be done faster...

for (i in 1:N) {

xtx <- outer(x[i,], x[i,])

G <- G + ( - y2[i] < xb[i] & xb[i] < y1[i] ) * xtx

V <- V + ((y2[i]^2) * (y1[i] <= xb[i]) +

(y1[i]^2) * (xb[i] <= - y2[i]) +

(y1[i] - y2[i] - xb[i])^2 * (- y2[i] < xb[i] & xb[i] < y1[i])

) * xtx

}

G <- G / N

V <- V / N

Ginv <- solve(G)

sqrt(diag(Ginv %*% V %*% Ginv / N)) # right?

}

# two-sided censoring

# y's must be censored at 0 and 1

alan <- function (b, x1, x2, y1, y2) {

d <- pmax(-1, pmin(as.matrix(x1-x2) %*% as.matrix(b), 1))

# sum of R(y1, y2, max(-1, min(,(x1-x2) %*% b, 1)))

r1 <- r2 <- rep(NA, length(d))

r1 <- ifelse(d <= y1 -1, d + d^2/2 +(y1-1)^2/2,

ifelse(y1-1 < d & d <= 0, d*y1,

ifelse(0 < d & d <= y1, d*y1 - d^2/2,

ifelse(y1 < d, 2,<="" div="" y1^2="">

NA))))

r2 <- ifelse(d <= -y2, -y2^2/2,

ifelse(-y2 < d & d <= 0, d^2/2 + d*y2,

ifelse(0 < d & d <= 1-y2, d*y2,

ifelse(1-y2 < d, (y2-1)^2="" -="" 2,<="" 2="" d="" d^2="" div="">

NA))))

- sum(r1 - r2) # Alan et al say minimize, but I say maximize, ie minimize the minus sum

}

fe2tobit <- function (x1, x2, y1, y2) {

initvals <- lm.fit(x1 - x2, y1 - y2)$coefficients

res <- optim(initvals, fn=alan, x1=x1, x2=x2, y1=y1, y2=y2,

method="BFGS", control=list(reltol=1e-8, maxit=1000))

if (res$convergence != 0) warning("Didn't converge")

res$par

}

fet2stderr <- function (x1, x2, y1, y2, bhat) {

N <- nrow(x1)

x <- x1 - x2

xb <- x %*% bhat

V <- G <- matrix(0, nrow=length(bhat), ncol=length(bhat))

u1 <- ifelse(xb > 0, pmax(y1-xb,0), pmin(y1,1+xb))

u2 <- ifelse(xb > 0, pmin(y2, 1-xb), pmax(y2+xb, 0))

u <- u1 - u2

gwt <- ifelse(xb < -1 | xb > 1, 0,

(-1 < xb & xb < y1-1) - (0 < xb & xb < y1) -

(-y1 < xb & xb < 0) + (1-y2 < xb & xb < 1))

for (i in 1:N) {

xtx <- outer(x[i,], x[i,])

G <- G + gwt[i] * xtx

if (-1 < xb[i] & xb[i] < 1) {

v <- u[i]*x[i,]/2

V <- V + outer(v, v)

}

}

G <- G / N

V <- V / N

Ginv <- solve(G)

sqrt(diag(Ginv %*% V %*% Ginv / N))

}

Points to note:

• fetobit expects data to be censored at 0; fe2tobit expects censoring at 0 and 1. There should be no NAs. 
• Warning! I have done only a very quick test on these. They seem to work but I am not a proper econometrician. Use at your own risk. 
• The standard errors may be inaccurate for low N. Consider bootstrap estimation. 
• Warning! Google searches for "Trimmed LAD" may lead to unexpected results.