potterzot · September 2, 2020 02:25
diff --git a/simulation.R b/simulation.R
 # 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)
	# 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)