Simulations of Variance Estimators for Rank-Rank Instrumental Variable Settings"
In csranks: Statistical Tools for Ranks

Data Generating Process

We use a data generating process given by \begin{equation}\label{eq:multinorm} \begin{pmatrix} Y \ X \ Z \end{pmatrix} \sim \mathcal{N}\left(\begin{pmatrix} 2 \ 1 \ 0 \end{pmatrix} , \begin{pmatrix} 1& \kappa_1 & \kappa_2 \ \kappa_1 & 1 & \kappa_3 \ \kappa_2 & \kappa_3 & 1 \end{pmatrix} \right). \end{equation}

We simulate data using the multivariate normal distribution given by \eqref{eq:multinorm} with parameters controlling correlations among variables. The parameters will take the following values.

library(MASS)
library(AER)
library(dplyr)
library(csranks)
theta_vals <- c(0.3, 0.7)
kappa_3_vals <- 2 * sin(pi / 6 * theta_vals)

# Compute corresponding k2 values to render rho = 0.5
kappa_2_vals <- 2 * sin(0.5 * asin(kappa_3_vals / 2))

# Compute k1 accordingly to acchieve different levels of wished Endogeneity
# Solve analytically when Endogeneity = 0
get_k1_zero <- function(kappa_2, kappa_3) {
  rho <- asin(kappa_2 / 2) / asin(kappa_3 / 2)
  return(2 * sin((pi / 6) * rho))
}

# Solve via optimization for other Endogeneity targets
get_k1_for_endog <- function(kappa_2, kappa_3, target_endog) {
  rho <- asin(kappa_2 / 2) / asin(kappa_3 / 2)

  f <- function(k1) {
    var_eps <- 1 / 12 + rho^2 / 12 - rho / pi * asin(k1 / 2)
    var_nu <- 1 / 12 + (6 / pi * asin(kappa_3 / 2))^2 / 12 - (6 / pi * asin(kappa_3 / 2)) / pi * asin(kappa_3 / 2)
    cov_eps_nu <- 1 / (2 * pi) * asin(k1 / 2) - 1 / 12 * rho

    if (var_eps <= 0 || var_nu <= 0) {
      return(Inf)
    }

    return(abs(target_endog - cov_eps_nu / (sqrt(var_eps) * sqrt(var_nu))))
  }

  return(optimize(f, interval = c(-1, 1))$minimum)
}

# Construct the grid
kappa_grid <- data.frame()
for (i in 1:length(kappa_3_vals)) {
  k3 <- kappa_3_vals[i]
  k2 <- kappa_2_vals[i]

  k1_zero <- get_k1_zero(k2, k3)
  k1_03 <- get_k1_for_endog(k2, k3, 0.3)
  k1_07 <- get_k1_for_endog(k2, k3, 0.7)

  kappa_grid <- rbind(
    kappa_grid,
    data.frame(
      kappa_1 = c(k1_zero, k1_03, k1_07),
      kappa_2 = k2,
      kappa_3 = k3,
      Endogeneity_Target = c(0, 0.3, 0.7)
    )
  )
}

save(kappa_grid, file = "kappa_grid.RData")

knitr::opts_chunk$set(echo = TRUE, eval = TRUE, message = FALSE, warning = FALSE)

params <- expand.grid(n = c(500, 1000, 5000)) %>%
  merge(kappa_grid, by = NULL) # Merge to replicate for each sample size

params$Instrument_Strength <- c(rep(0.3, 9), rep(0.7, 9))
params$mu <- list(c(2, 1, 0))

We generate the ranked data $F_Y(Y), F_X(X), F_Z(Z)$ through multinormal samples that are then ranked via the Gaussian cumulative distribution functions $\Phi(x)$.

gen_data <- function(n, kappa_1, kappa_2, kappa_3, mu) {
  # Define covariance matrix
  Sigma <- matrix(
    c(
      1, kappa_1, kappa_2,
      kappa_1, 1, kappa_3,
      kappa_2, kappa_3, 1
    ),
    nrow = 3, byrow = TRUE
  )

  samples <- mvrnorm(n = n, mu = mu, Sigma = Sigma)

  # Create dataframe and transform to ranks in [0,1]
  df <- data.frame(Y = samples[, 1], X = samples[, 2], Z = samples[, 3])
  df$rank_Y <- pnorm(df$Y - mu[1])
  df$rank_X <- pnorm(df$X - mu[2])
  df$rank_Z <- pnorm(df$Z - mu[3])

  df$emp_rank_Y <- rank(df$Y) / length(df$Y)
  df$emp_rank_X <- rank(df$X) / length(df$X)
  df$emp_rank_Z <- rank(df$Z) / length(df$Z)
  return(df)
}

Simulations

We now established the estimation of the asymptotic variance and can thus proceed to the simulations. For every combination of parameters we generate $S=1000$ data frames according to the data generating process described above. The following two functions then obtain the IV parameter estimate $\hat \rho$ and the estimated variance from the vcov-function:

S <- 10000
# Function to compute IV estimate
compute_iv <- function(df) {
  iv_model <- ivreg::ivreg(emp_rank_Y ~ emp_rank_X | emp_rank_Z, data = df)
  return(coef(iv_model))
}

For every parameter combination and every data frame generated within, we store the empirical variance of the estimated parameters $\hat \rho$ and the mean of the estimated variances according to vcov, together with the combination of parameters in a list:

# Run simulations
results <- list()
estm_list <- numeric(nrow(params))

j <- 1
estm <- matrix(NA, nrow = S, ncol = 2)

for (s in 1:S) {
  df <- gen_data(
    n = params$n[j],
    kappa_1 = params$kappa_1[j],
    kappa_2 = params$kappa_2[j],
    kappa_3 = params$kappa_3[j],
    mu = params$mu[[j]]
  )

  estm[s, ] <- compute_iv(df) # Store IV estimate of rho
}

# Store empirical variance
results <- list(
  params = params[j, ],
  var_iv_emp = params$n[j] * cov(estm)
)
estm_list <- estm