Friday, 15 July 2011

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.