Last active
September 2, 2020 02:25
-
-
Save potterzot/de27545f36f7ab8d606d7a398541bef2 to your computer and use it in GitHub Desktop.
Replication of Lewbel, Schennach, and Zhang (2020): Identification of a Triangular Two Equation System Without Instruments
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# This script attempts to replicate a potion of Lewbel et al (2020) | |
# Working paper: https://drive.google.com/file/d/1UHWtjK81Cd0f3ksg6pv9B-xg1TGWNRQ1/view | |
# We focus on the simulation design #2 described on pg. 15 | |
# See also Appendix B on pg. 28 for a clear description of estimation approach | |
# Results are in Table 2 on pg. 33 | |
# Number of observations and true coefficient values | |
set.seed(15) | |
N <- 100 | |
beta <- 1 | |
gamma <- 1 | |
### Function to generate data of the form: | |
# Y = \beta_1 U + V | |
# W = \gamma Y + \beta_2 U + R | |
# log(U) ~ N(-0.5, 1) | |
# V ~ Unif(-2, 2) | |
# R ~ Unif(-2, 2) | |
make_simdata <- function() { | |
dt <- data.frame( | |
u = exp(rnorm(N, -0.5, 1)), | |
v = runif(N, -2, 2), | |
r = runif(N, -2, 2)) %>% | |
dplyr::mutate( | |
y = u + v, | |
w = gamma * y + beta * u + r) %>% | |
dplyr::select(w, y, u, v, r) %>% | |
dplyr::mutate( | |
y = scale(y, scale = FALSE), | |
w = scale(w, scale = FALSE) | |
) | |
} | |
### Moment function | |
g <- function(theta, dt) { | |
# Extract the data | |
W <- data.matrix(dt[,1]) | |
Y <- data.matrix(dt[,2]) | |
# Extract parameters | |
gamma <- theta[1] | |
idx <- 1 | |
B <- exp(theta[idx + 1]) | |
s2U <- exp(theta[idx + 2]) | |
s2V <- exp(theta[idx + 3]) | |
s2R <- exp(theta[idx + 4]) | |
#uWW <- theta[idx + 5] # For overidentification, not used | |
# Convenience definitions to standardize moments | |
uYY <- s2U + s2V | |
uYW <- B * s2U + gamma * uYY | |
a <- gamma + B | |
Q <- W - Y * gamma | |
P <- W - Y * a | |
YY <- Y * Y | |
YW <- Y * W | |
# Moments. See Appendix B of Lewbel et al. (2020), eq. 74-77. | |
m1 <- YW - uYW | |
m2 <- YY - uYY | |
m3 <- Q * P * Y | |
m4 <- Q * P * m2 - 2 * B * s2U * P * Y | |
m5 <- Q^2 - B^2 * s2U - s2R | |
M <- cbind(m1, m2, m3, m4, m5) | |
return(M) | |
} | |
# Generate data and estimate 1000 times | |
res <- lapply(1:1000, function(i) { | |
dt <- make_simdata(method = 2) | |
# Starting values that match the true values and the mean values found in Lewbel et al. (2020) | |
# We use the exponentiated values to enforce > 0 requirement | |
theta <- c(gamma = 1, tB = log(1), tU = log(1.72), tV = log(1.64), tR = log(1.64)) | |
# GMM estimation | |
m <- gmm(g, dt, t0 = theta, | |
itermax = 500, | |
vcov = "iid", | |
control=list( | |
fnscale=1e-8 | |
)) | |
c(coef(m)[1], exp(coef(m)[-1])) # extract the coefficients | |
}) %>% | |
dplyr::bind_rows() | |
# Coefficient means. Should be similar to: | |
# Lewbel (2020): c(gamma = 1.09, beta = 1.02, s2U = 1.11, s2V = 1.90, s2R = 1.49) | |
# My results : c(gamma = 1.17, beta = 0.86, s2U = 1.32, s2V = 1.71, s2R = 1.57) | |
colMeans(res) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment