Nothing
#' WRE13, WRE33, WRU13 and WRU33 bootstraps as in MNW (2022) "Fast and
#' reliable" and MacKinnon "Fast Cluster Bootstraps"
#' (Econometrics & Statistics, 2021)
#'
#' @param preprocessed_object A list: output of the preprocess2 function.
#' @param B number of bootstrap iterations
#' @param r Shifts the null hypothesis.
#' @param sign_level The significance level.
#' @param param name of the test parameter.
#' @param p_val_type type Type of p-value. By default "two-tailed".
#' Other options: "equal-tailed", ">", "<"
#' @param nthreads The number of threads. Can be: a) an integer lower than,
#' or equal to, the maximum number of threads; b) 0: meaning
#' all available threads will be used; c) a number strictly
#' between 0 and 1 which represents the fraction of all
#' threads to use. The default is to use 50\% of all
#' threads. You can set permanently the number of threads
#' used within this package using the function ...
#' @param type character or function. The character string specifies the
#' type of boostrap to use: One of "rademacher", "mammen", "norm"
#' and "webb". Alternatively, type can be a function(n) for drawing
#' wild bootstrap factors. "rademacher" by default.
#' @param full_enumeration Is full enumeration employed? Full enum.
#' is used if N_G^2 < B for Mammen and Rademacher weights
#' @param small_sample_correction The small sample correction to be applied.
#' See ssc().
#' If FALSE, run wild cluster bootstrap.
#' @param object the regression object
#' @param impose_null logical scalar. Should the null be imposed on the
#' bootstrap dgp or not?
#' @importFrom MASS ginv
#' @importFrom summclust vcov_CR3J
#' @return A list of bootstrap results.
#' @noRd
boot_algo_fastnreliable <- function(
preprocessed_object,
B,
bootstrap_type,
r = 0,
sign_level,
param,
p_val_type,
nthreads,
type,
full_enumeration,
small_sample_correction,
object,
impose_null,
sampling){
#here for debugging
#preprocessed_object <- preprocess
if(substr(bootstrap_type, 2, 2) == 1){
crv_type <- "crv1"
} else {
crv_type <- "crv3"
}
X <- preprocessed_object$X
y <- preprocessed_object$Y
R <- preprocessed_object$R0
cluster_df <- preprocessed_object$clustid
clustid <- names(cluster_df)
fe <- preprocessed_object$fe
cluster <- as.factor(cluster_df[,1])
bootcluster <- preprocessed_object$bootcluster
G <- N_G_bootcluster <- length(unique(bootcluster[[1]]))
k <- length(R)
bootstrap_type_x <- paste0(substr(bootstrap_type, 1, 1), "x")
if(impose_null){
# WCR1x or WCR3x
bootstrap_type <- paste0("WCR", bootstrap_type_x)
} else {
# WCU1x or WCU3x
bootstrap_type <- paste0("WCU", bootstrap_type_x)
}
v <- get_weights(
type = type,
full_enumeration = full_enumeration,
N_G_bootcluster = N_G_bootcluster,
boot_iter = B,
sampling = sampling
)
# create X_g's, X1_g's, y_g's etc
X_list <- matrix_split(X, cluster, "row")
y_list <- split(y, cluster, drop = FALSE)
# precompute a range of other objects
tXgXg <- lapply(
seq_along(1:G),
function(g) crossprod(X_list[[g]])
)
tXgyg <- lapply(
seq_along(1:G),
function(g) t(X_list[[g]]) %*% y_list[[g]]
)
tXX <- Reduce("+", tXgXg) # crossprod(X)
tXy <- Reduce("+", tXgyg) # t(X) %*% y
tXXinv <- solve(tXX)
RtXXinv <- R %*% tXXinv
tXgX1g <- NULL
beta_hat <- NULL
beta_tilde <- NULL
beta_g_hat <- NULL
beta_1g_tilde <- NULL
inv_tXX_tXgXg <- NULL
if(bootstrap_type %in% c("WCR3x", "WCU3x")){
# X1: X without parameter beta for which hypothesis beta = 0 is tested
X1 <- X[, which(R == 0)]
X1_list <- matrix_split(X1, cluster, "row")
tX1gX1g <- lapply(
seq_along(1:G),
function(g) crossprod(X1_list[[g]])
)
tX1gyg <- lapply(
seq_along(1:G),
function(g) t(X1_list[[g]]) %*% y_list[[g]]
)
tXgX1g <- lapply(
seq_along(1:G),
function(g) t(X_list[[g]]) %*% X1_list[[g]]
)
tX1X1 <- Reduce("+", tX1gX1g) # crossprod(X1)
tX1y <- Reduce("+", tX1gyg) #t(X1) %*% y
tX1X1inv <- solve(tX1X1)
}
if(bootstrap_type == "WCR1x"){
beta_hat <- tXXinv %*% tXy
beta_tilde <- beta_hat -
tXXinv %*% R %*% solve(t(R) %*% tXXinv %*% R) %*% (R %*% beta_hat - 0)
} else if (bootstrap_type == "WCU1x"){
beta_hat <- tXXinv %*% tXy
} else if (bootstrap_type == "WCR3x"){
inv_tXX_tXgXg <- lapply(
1:G,
function(x) MASS::ginv(tXX - tXgXg[[x]])
)
beta_1g_tilde <- lapply(
1:G,
function(g) MASS::ginv(tX1X1 - tX1gX1g[[g]]) %*% (tX1y - tX1gyg[[g]])
)
} else if(bootstrap_type == "WCU3x"){
beta_g_hat <- lapply(
1:G,
function(g) MASS::ginv(tXX - tXgXg[[g]]) %*% (tXy - tXgyg[[g]])
)
}
if(crv_type == "crv1"){
if(is.null(beta_hat)){
beta_hat <- tXXinv %*% tXy
}
} else {
if(is.null(inv_tXX_tXgXg)){
inv_tXX_tXgXg <- lapply(
1:G,
function(x) MASS::ginv(tXX - tXgXg[[x]])
)
}
}
# compute scores
scores_list <- get_scores(
bootstrap_type = bootstrap_type,
G = G,
tXgyg = tXgyg,
tXgXg = tXgXg,
tXgX1g = tXgX1g,
beta_hat = beta_hat,
beta_tilde = beta_tilde,
beta_g_hat = beta_g_hat,
beta_1g_tilde = beta_1g_tilde
)
# pre-allocate space for bootstrap
# start the bootstrap loop
t_boot <- t_boot2 <- vector(mode = "numeric", B + 1)
se <- se2 <- vector(mode = "numeric", B + 1)
dim(R) <- c(1, k) # turn R into matrix
Cg <- R %*% tXXinv %*% Reduce("cbind", scores_list)
numer <- Cg %*% v
if(crv_type == "crv1"){
H <- matrix(NA, G, G)
for(g in 1:G){
for(h in 1:G){
H[g,h] <- R %*% tXXinv %*% tXgXg[[g]] %*% tXXinv %*% scores_list[[h]]
}
}
denom <- boot_algo3_crv1_denom(
B = B,
G = G,
ssc = small_sample_correction,
H = H,
Cg = Cg,
v = v,
cores = nthreads
)
t_boot <- c(numer / sqrt(c(denom)))
} else if (crv_type == "crv3"){
for(b in 1:(B + 1)){
# Step 1: get bootstrapped scores
scores_g_boot <- matrix(NA,G,k)
v_ <- v[,b]
for(g in 1:G){
scores_g_boot[g,] <- scores_list[[g]] * v_[g] #* v[g, b]
}
# numerator (both for WCR, WCU)
scores_boot <- colSums(scores_g_boot)
delta_b_star <- tXXinv %*% scores_boot
delta_diff <- matrix(NA, G, k)
for(g in 1:G){
score_diff <- scores_boot - scores_g_boot[g,]
delta_diff[g,] <-
(
(inv_tXX_tXgXg[[g]] %*% score_diff) - delta_b_star
)^2
}
se[b] <-
sqrt(
((G-1) / G) *
colSums(
delta_diff
)
)[which(R == 1)]
t_boot[b] <- c(delta_b_star)[which(R == 1)] / se[b]
}
}
# get original t-stat.
if(crv_type == "crv1"){
score_all <- lapply(
1:G, function(g)
tcrossprod(
crossprod(X_list[[g]], y_list[[g]] - X_list[[g]] %*% beta_hat)
)
)
meat <- Reduce("+", score_all)
vcov <- tXXinv %*% meat %*% tXXinv
} else if(crv_type == "crv3"){
vcov3 <- quote(
summclust::vcov_CR3J(
obj = object,
cluster = clustid
)
)
if(inherits(object, "fixest")){
vcov3$absorb_cluster_fixef <- FALSE
}
vcov <- eval(vcov3)
}
se0 <- sqrt(small_sample_correction * R %*% vcov %*% t(R))
se0 <- as.vector(se0)
t_stat <- as.vector(
coef(object)[which(R == 1)] / se0
)
t_boot <- t_boot[-1]
p_val <-
get_bootstrap_pvalue(
p_val_type = p_val_type,
t_stat = t_stat,
t_boot = t_boot
)
res <- list(
p_val = p_val,
t_stat = t_stat,
t_boot = t_boot,
B = B,
R0 = R,
param = param,
clustid = clustid,
invalid_t = NULL,
ABCD = NULL,
small_sample_correction = small_sample_correction
)
class(res) <- "boot_algo3"
invisible(res)
}
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.