R/IdentifiedSet.R

In fditraglia/ivdoctr: Ensures Mutually Consistent Beliefs When Using IVs

#' @import stats
NULL

#' Given data and function specification, returns the relevant correlations and
#'   covariances with any exogenous controls projected out.
#' @param y_name Name of the dependent variable
#' @param T_name Name(s) of the preferred regressor(s)
#' @param z_name Name(s) of the instrumental variable(s)
#' @param data Data to be analyzed
#' @param controls Exogenous regressors to be included
#'
#' @return List of correlations, covariances, and R^2 of first and second stage
#'   regressions after projecting out any exogenous control regressors
get_observables <- function(y_name, T_name, z_name, data, controls = NULL) {

  # Project out control regressors if present
  if (!is.null(controls)) {
    y <- resid(lm(reformulate(controls, response = y_name), data))
    T_reg <- lm(reformulate(controls, response = T_name), data)
    z_reg <- lm(reformulate(controls, response = z_name), data)
    Tobs <- resid(T_reg)
    T_Rsq <- summary(T_reg)$r.squared
    z <- resid(z_reg)
    z_Rsq <- summary(z_reg)$r.squared
  } else {
    y <- get(y_name, data)
    Tobs <- get(T_name, data)
    z <- get(z_name, data)
    T_Rsq <- z_Rsq <- 0 # No controls is the same as controls that are
                        # uncorrelated with both Tobs and z
  }

  Sigma <- cov(cbind(Tobs, y, z))
  s2_T <- Sigma["Tobs", "Tobs"]
  s2_y <- Sigma["y", "y"]
  s2_z <- Sigma["z", "z"]
  s_Ty <- Sigma["Tobs", "y"]
  s_Tz <- Sigma["Tobs", "z"]
  s_zy <- Sigma["z", "y"]

  Rho <- cov2cor(Sigma)
  r_Ty <- Rho["Tobs", "y"]
  r_Tz <- Rho["Tobs", "z"]
  r_zy <- Rho["z", "y"]

  list(n = nrow(data),
       T_Rsq = T_Rsq,
       z_Rsq = z_Rsq,
       s2_T = s2_T,
       s2_y = s2_y,
       s2_z = s2_z,
       s_Ty = s_Ty,
       s_Tz = s_Tz,
       s_zy = s_zy,
       r_Ty = r_Ty,
       r_Tz = r_Tz,
       r_zy = r_zy)
}

#' Given observables from the data, generates unrestricted bounds for kappa.
#'   Vectorized
#'
#' @param obs Observables generated by get_observables
#' @param tilde Boolean of whether or not kappa_tilde or kappa is desired
#'
#' @return List of upper bounds and lower bounds for kappa
get_k_bounds_unrest <- function(obs, tilde) {
  k_tilde <- with(obs, (r_Ty ^ 2 + r_Tz ^ 2 - 2 * r_Ty * r_Tz * r_zy) / (1 - r_zy ^ 2))
  lower_bound <- k_tilde * tilde + ((1 - obs$T_Rsq) * k_tilde + obs$T_Rsq) * (1 - tilde)
  ans <- list(Lower = lower_bound, Upper = rep(1, length(k_tilde)))
  return(ans)
}

#' Given observables from the data, generates the unrestricted bounds for rho_uz.
#'   Vectorized
#'
#' @param obs Observables generated by get_observables
#'
#' @return List of upper and lower bounds for rho_uz
get_r_uz_bounds_unrest <- function(obs) {
  k_tilde_lower <- get_k_bounds_unrest(obs, tilde = TRUE)$Lower
  bound <- abs(obs$r_Tz) / sqrt(k_tilde_lower)
  nontrivial_lower_bound <- with(obs, r_Ty * r_Tz < k_tilde_lower * r_zy)
  upper_bound <- ifelse(nontrivial_lower_bound, 1, bound)
  lower_bound <- ifelse(nontrivial_lower_bound, -1 * bound, -1)
  ans <- list(Lower = lower_bound, Upper = upper_bound)
  return(ans)
}

#' Given observables from the data, generates the unrestricted bounds for rho_TstarU.
#'   Data does not impose any restrictions on r_TstarU
#'   Vectorized
#'
#' @param obs Observables generated by get_observables
#'
#' @return List of upper and lower bounds for r_TstarU
get_r_TstarU_bounds_unrest <- function(obs) {
  return(list(Lower = rep(-1, length(obs[[1]])), Upper = rep(1, length(obs[[1]]))))
}

#' Wrapper function combines all unrestricted bounds together. Vectorized
#'
#' @param obs Observables generated by get_observables
#'
#' @return List of unrestricted bounds for r_TstarU, r_uz, and kappa
get_bounds_unrest <- function(obs) {
  ans <- list(r_TstarU = get_r_TstarU_bounds_unrest(obs),
              r_uz = get_r_uz_bounds_unrest(obs),
              k = get_k_bounds_unrest(obs, tilde = FALSE),
              k_tilde = get_k_bounds_unrest(obs, tilde = TRUE))
  return(ans)
}

#' Solves for r_uz given observables, r_TstarU, and kappa
#'
#' This function solves for r_uz given r_TstarU and kappa. It handles 3
#'   potential cases when r_uz must be evaluated:
#'   1. Across multiple simulations, but given the same r_TstarU and k
#'   2. For multiple simulations, each with a value of r_TstarU and k
#'   3. For one simulation across a grid of r_TstarU and k
#' @param r_TstarU Vector of r_TstarU values
#' @param k Vector of kappa values
#' @param obs Observables generated by get_observables
#'
#' @return Vector of r_uz values.
get_r_uz <- function(r_TstarU, k, obs) {
  if ((length(r_TstarU) == 1 & length(k) == 1) |
      (length(r_TstarU) == length(obs$r_Ty) & length(r_TstarU == length(k))) |
      (length(obs$r_Ty) == 1 & length(r_TstarU) == length(k))) {
    if (length(r_TstarU) == 1 & length(k) == 1) {
      k <- rep(k, length(obs$r_Ty))
      r_TstarU <- rep(r_TstarU, length(obs$r_Ty))
    }
    A <- with(obs, r_TstarU * r_Tz / sqrt(k))
    B1 <- with(obs, r_Ty * r_Tz - k * r_zy)
    B2 <- with(obs, sqrt((1 - r_TstarU ^ 2) / (k * (k - r_Ty ^ 2))))
    A - B1 * B2
  } else {
    stop("Function is only designed to evaluate r_uz across multiple simulations
         for the same r_TstarU and k, across multiple simulations each with
         a different r_TstarU and k, or for the same simulation across a grid
         of r_TstarU and k values. Please fix your input r_TstarU and k to
         either only take 1 value or match the number of simulations run.")
  }
}

#' Solves for the magnification factor
#'
#' This function solves for the magnification factor given r_TstarU and kappa.
#'   It handles 3 potential cases when the magnification factor must be evaluated:
#'   1. Across multiple simulations, but given the same r_TstarU and k
#'   2. For multiple simulations, each with a value of r_TstarU and k
#'   3. For one simulation across a grid of r_TstarU and k
#' @param r_TstarU Vector of r_TstarU values
#' @param k Vector of kappa values
#' @param obs Observables generated by get_observables
#'
#' @return Vector of magnification factors
get_M <- function(r_TstarU, k, obs) {
  if ((length(r_TstarU) == 1 & length(k) == 1) |
      (length(r_TstarU) == length(obs$r_Ty) & length(r_TstarU == length(k))) |
      (length(obs$r_Ty) == 1 & length(r_TstarU) == length(k))) {
    if (length(r_TstarU) == 1 & length(k) == 1) {
      k <- rep(k, length(obs$r_Ty))
      r_TstarU <- rep(r_TstarU, length(obs$r_Ty))
    }
  B1 <- with(obs, r_Ty * r_Tz - k * r_zy)
  B2 <- with(obs, sqrt((1 - r_TstarU ^ 2) / (k * (k - r_Ty ^ 2))))
  dr_TstarU <- with(obs, r_Tz / sqrt(k) + r_TstarU * B1 /
                    sqrt(k * (k - r_Ty ^ 2) * (1 - r_TstarU ^ 2)))
  dk <- with(obs, -r_TstarU * r_Tz / (2 * k ^ (3 / 2)) + B2 *
                  (r_zy + B1 / 2.0 * (1 / k + 1 / (k - r_Ty ^ 2))))
  sqrt(1 + dr_TstarU ^ 2 + dk ^ 2)
  } else {
    stop("Function is only designed to evaluate the magnification factor across
          multiple simulations for the same r_TstarU and k, across multiple
          simulations each with a different r_TstarU and k, or for the same
          simulation across a grid of r_TstarU and k values. Please fix your
          input r_TstarU and k to either only take 1 value or match the number
          of simulations run.")
  }
}

#' Solves for the variance of the error term u
#'
#' This function solves for the variance of u given r_TstarU and kappa.
#'   It handles 3 potential cases when the variance of u must be evaluated:
#'   1. Across multiple simulations, but given the same r_TstarU and k
#'   2. For multiple simulations, each with a value of r_TstarU and k
#'   3. For one simulation across a grid of r_TstarU and k
#' @param r_TstarU Vector of r_TstarU values
#' @param k Vector of kappa values
#' @param obs Observables generated by get_observables
#'
#' @return Vector of variances of u
get_s_u <- function(r_TstarU, k, obs) {
  if ((length(r_TstarU) == 1 & length(k) == 1) |
      (length(r_TstarU) == length(obs$r_Ty) & length(r_TstarU == length(k))) |
      (length(obs$r_Ty) == 1 & length(r_TstarU) == length(k))) {
    if (length(r_TstarU) == 1 & length(k) == 1) {
      k <- rep(k, length(obs$r_Ty))
      r_TstarU <- rep(r_TstarU, length(obs$r_Ty))
    }
    with(obs, sqrt(s2_y * (k - r_Ty ^ 2) / (k * (1 - r_TstarU ^ 2))))
  } else {
    stop("Function is only designed to evaluate s_u across multiple simulations
         for the same r_TstarU and k, across multiple simulations each with
         a different r_TstarU and k, or for the same simulation across a grid
         of r_TstarU and k values. Please fix your input r_TstarU and k to
         either only take 1 value or match the number of simulations run.")
  }
}

#' Solves for beta
#'
#' This function solves for beta given r_TstarU and kappa.
#'   It handles 3 potential cases when beta must be evaluated:
#'   1. Across multiple simulations, but given the same r_TstarU and k
#'   2. For multiple simulations, each with a value of r_TstarU and k
#'   3. For one simulation across a grid of r_TstarU and k
#' @param r_TstarU Vector of r_TstarU values
#' @param k Vector of kappa values
#' @param obs Observables generated by get_observables
#'
#' @return Vector of betas
get_beta <- function(r_TstarU, k, obs) {
  if ((length(r_TstarU) == 1 & length(k) == 1) |
      (length(r_TstarU) == length(obs$r_Ty) & length(r_TstarU == length(k))) |
      (length(obs$r_Ty) == 1 & length(r_TstarU) == length(k))) {
    if (length(r_TstarU) == 1 & length(k) == 1) {
      k <- rep(k, length(obs$r_Ty))
      r_TstarU <- rep(r_TstarU, length(obs$r_Ty))
    }
    r_uz <- get_r_uz(r_TstarU, k, obs)
    s_u <- get_s_u(r_TstarU, k, obs)
    with(obs, (r_zy * sqrt(s2_y) - r_uz * s_u) / (r_Tz * sqrt(s2_T)))
  } else {
    stop("Function is only designed to evaluate beta across multiple simulations
         for the same r_TstarU and k, across multiple simulations each with
         a different r_TstarU and k, or for the same simulation across a grid
         of r_TstarU and k values. Please fix your input r_TstarU and k to
         either only take 1 value or match the number of simulations run.")
  }
}

# This function is vectorized wrt k_min, k_max and obs.
get_beta_lower <- function(r_TstarU_max, k_min, k_max, obs) {
  # min(beta) occurs at max(r_TstarU) but could be a corner value for kappa
  beta_corner <- pmin(get_beta(r_TstarU_max, k_min, obs),
                      get_beta(r_TstarU_max, k_max, obs))

  # If r_TstarU_max = 0, no need to check the interior solution
  if (identical(r_TstarU_max, 0)) {
    out <- beta_corner
  } else {
    C <- 1 / sqrt(1 + (r_TstarU_max ^ 2 / (1 - r_TstarU_max ^ 2)))
    k_interior <- with(obs, 2 * r_Ty ^ 2 / (1 - sign(r_TstarU_max * r_Ty) * C))

    # It only makes sense to calculate beta_interior if k_interior is between
    # k_min and k_max. If not, set it to Inf so it always exceeds beta_corner.
    k_interior_is_valid <- (k_interior > k_min) & (k_interior < k_max)
    beta_interior <- ifelse(k_interior_is_valid,
                            get_beta(r_TstarU_max, k_interior, obs), Inf)
    out <- pmin(beta_corner, beta_interior)
  }
  return(out)
}

# This function is vectorized wrt k_min, k_max, and obs.
get_beta_upper <- function(r_TstarU_min, k_min, k_max, obs) {
  # max(beta) occurs at min(r_TstarU) but could be a corner value for kappa
  beta_corner <- pmax(get_beta(r_TstarU_min, k_min, obs),
                      get_beta(r_TstarU_min, k_max, obs))

  # If r_TstarU_min = 0, no need to check interior solution
  if (identical(r_TstarU_min, 0)) {
    out <- beta_corner
  } else {
    C <- 1 / sqrt(1 + (r_TstarU_min ^ 2 / (1 - r_TstarU_min ^ 2)))
    k_interior <- with(obs, 2 * r_Ty ^ 2 / (1 + sign(r_TstarU_min * r_Ty) * C))

    # It only makes sense to calculate beta_interior if k_interior is between
    # k_min and k_max. If not, set it to -Inf so it never exceeds beta_corner.
    k_interior_is_valid <- (k_interior > k_min) & (k_interior < k_max)
    beta_interior <- ifelse(k_interior_is_valid,
                            get_beta(r_TstarU_min, k_interior, obs), -Inf)
    out <- pmax(beta_corner, beta_interior)
  }
  return(out)
}