Nothing
#' @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)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.