library("LARF") library("Rlab") library("MSBVAR") library("systemfit") set.seed(1111) sim <- function(J = 200, N = 1000, outcome = "binary") { tab <- matrix(NA, J, 6) em <- 0 for (j in 1:J) { # Covariates X <- matrix(NA, N, 2) for (i in 1:2) { X[, i] <- rnorm(N, 0, i) } # Binary instrument Z <- rbern(N, 0.5) # Correlated errors error <- rmultnorm(N, matrix(c(0, 0), 2, 1), vmat = matrix(c(1, 0.8, 0.8, 1), 2, 2)) # Binary treatment variable B1 <- c(0.3, -0.5) B2 <- 0.7 D.star <- -0.1 + X %*% as.matrix(B1) + B2 * Z + error[, 1] D <- ifelse(D.star > 0, 1, 0) # sum(D) / N # with large N, we should be able to recover all the true parameters in the first # equation eq1 <- glm(D ~ Z + X, family = binomial(probit)) # summary(eq1) # Binary outcome B3 <- c(0.4, -0.6) B4 <- 0.8 Y.star <- -0.2 + X %*% as.matrix(B3) + B4 * D + error[, 2] Y <- ifelse(Y.star > 0, 1, 0) # sum(Y) / N if (outcome == "binary") { # Even with large N, we should NOT be able to recover the true parameter D in the # outcome equation eq2 <- glm(Y ~ D + X, family = binomial(probit)) tab[j, 1:2] <- summary(eq2)$coefficients[2, 1:2] ### Test LARF NX <- cbind(1, X) lf <- try(larf(Y ~ X, D, Z), TRUE) if (class(lf) == "larf") { tab[j, 3:4] <- c(lf$coefficients[1], lf$StdErr[1]) } else { tab[j, 3:4] <- NA em <- em + 1 } ### 2sls fit2sls <- systemfit(Y ~ D + X, method = "2SLS", inst = ~Z + X) tab[j, 5] <- fit2sls$coefficients[2] tab[j, 6] <- sqrt(fit2sls$coefCov[2, 2]) # print(c('Current Iteration:', j), quote = F) } # For continuous outcomes if (outcome == "continuous") { eq2 <- lm(Y.star ~ D + X) tab[j, 1:2] <- summary(eq2)$coefficients[2, 1:2] ### Test LARF NX <- cbind(1, X) lf <- try(larf(Y.star ~ X, D, Z), TRUE) if (class(lf) == "larf") { tab[j, 3:4] <- c(lf$coefficients[1], lf$StdErr[1]) } else { tab[j, 3:4] <- NA em <- em + 1 } ### 2sls fit2sls <- systemfit(Y.star ~ D + X, method = "2SLS", inst = ~Z + X) tab[j, 5] <- fit2sls$coefficients[2] tab[j, 6] <- sqrt(fit2sls$coefCov[2, 2]) # print(c('Current Iteration:', j), quote = F) } } table <- t(as.matrix(apply(cbind(tab[, 1] - 0.8, (tab[, 1] - 0.8)^2, tab[, 3] - 0.8, (tab[, 3] - 0.8)^2, tab[, 5] - 0.8, (tab[, 5] - 0.8)^2), 2, mean, trim = 0.1, na.rm = T))) colnames(table) <- c("b", "mse", "LARF", "mse", "2SLS", "mse") output <- list(table = table, em = em, J = J, N = N, tab = tab) return(output) } sim1 <- sim(J = 500, N = 2000, outcome = "continuous") tab1 <- sim1$table round(tab1, digits = 2) # Row 1 in Table 2 sim2 <- sim(J = 500, N = 2000) tab2 <- sim2$table round(tab2, digits = 2) # Row 2 in Table 2 ### Empirical exmaple. To reproduce the LARF estimates in Table 4. library("LARF") data("c401k", package = "LARF") attach(c401k) ## Continuous outcome. Treatment effects on net family financial assest est1 <- larf(nettfa ~ inc + age + agesq + marr + fsize, treatment = p401k, instrument = e401k, data = c401k) summary(est1) ### Nonparametric power series estimates of the probability of receiving the ### treatment inducement ps <- npse(e401k ~ inc + age + agesq + marr + fsize, order = 5) ps$Lambda est2 <- larf(nettfa ~ inc + age + agesq + marr + fsize, treatment = p401k, instrument = e401k, data = c401k, zProb = ps$fitted) summary(est2) # Binary outcome. Treatment effects on participation in IRA est3 <- larf(pira ~ inc + age + agesq + marr + fsize, p401k, e401k, data = c401k) summary(est3) # Average marginal effects est4 <- larf(pira ~ inc + age + agesq + marr + fsize, p401k, e401k, data = c401k, AME = TRUE) summary(est4) # MLE est5 <- larf(pira ~ inc + age + agesq + marr + fsize, p401k, e401k, data = c401k, method = "ML") summary(est5) ### The remaining results are computed in Stata: ### => see v71c01.do