#' General Inference
#' Testing general hypothesis regarding the qunatile process of the endogenous
#' variable.
#' @param object An ivqr object returned from the function \code{ivqr()}
#' @param variable A number indicates which endogenous variable to test. Since
#' at most two endongenous variables can be included in the function \code{ivqr},
#' this argument should be either 1 or 2.
#' @param trim A vector of two numbers indicating the lower and upper bounds
#' of the quantiles to consider.
#' @param B Number of sub-sampling in the bootstrap. Default is 2000.
#' @param size A vector indicates the desired size of the test. Critical values
#' will be reported accordingly.
#' @param nullH The null hypothesis to test. The four options are: No_Effect,
#' Dominance, Location_Shift, and Exogeneity as defined in
#' Chernozhukov and Hansen (2006).
#' @return An ivqr_ks object which contains information regarding test statistics,
#' critical value, sub-sampling block size, ...etc.
#' @references Chernozhukov, V., & Hansen, C. (2006). Instrumental quantile regression inference for structural and treatment effect models. Journal of Econometrics, 132(2), 491-525.
#' @examples
#' data(ivqr_eg)
#' fit <- ivqr(y ~ d | z | x, seq(0.15,0.85,0.02), grid = seq(-2,2,0.2), data = ivqr_eg) # taus should be a fine grid
#' ivqr(fit,nullH=No_Effect) # Test of no effect.
#' ivqr(fit,nullH=Dominance) # Test of dominance.
#' ivqr(fit,nullH=Location_Shift) # Test of location shift.
#' ivqr(fit,nullH=Exogeneity) # Test of exogeneity.
#' @export
ivqr.ks <- function(object, variable = NULL, trim = c(0.05,0.95), B = 2000,
size = c(0.1,0.05,0.01,0.001), nullH="No_Effect", ...){
if (any(object$error_tau_flag)) Stop("Error occurred for some tau. Re-specify taus
and run ivqr() again")
dim_d <- object$dim_d_d_k[1]
if (dim_d > 1 & is.null(variable)) {
warning("Please specify which (endogenous) variable to test. Default
is the first endogenous variable.")
}
if (dim_d == 1 & is.null(variable)) variable = 1
ks <- switch(nullH,
No_Effect = ivqr.ks.no(object, trim, B, variable, size),
Dominance = ivqr.ks.dom(object, trim, B, variable, size),
Location_Shift = ivqr.ks.const(object, trim, B, variable, size),
Exogeneity = ivqr.ks.exog(object, trim, B, variable, size),
"This test is not implemented")
return(ks)
}
ivqr.ks.no <- function(object, trim, B, variable, size) {
taus <- object$taus
data <- object$data
coef <- object$coef
fitted <- object$fitted
residuals <- object$residuals
endg_var_se <- object$se[variable]
dim_dk <- object$dim_d_d_k[1] + object$dim_d_d_k[3]
n <- object$n
PSI <- object$PSI
J <- object$vc$J
# nrow = n, ncol = length(taus)
L <- matrix(rep(taus,n), n, length(taus), byrow = TRUE) - as.numeric(residuals < 0)
Z <- matrix(NA,n,length(taus))
for (tau_index in 1:length(taus)) {
invJ <- solve(J[,,tau_index])
R_invJ <- invJ[variable,] # = R %*% invJ, R = c(0,0,,,0,1,0,..0)
Z[,tau_index] <- as.vector(L[,tau_index]) * (PSI %*% R_invJ)
}
sd_z <- colMeans(Z ^ 2) ^ (1 / 2)
block_size <- as.integer(5 * n ^ (2/5))
V <- matrix(NA, B, length(taus))
for (b in 1:B) {
resample_indexes <- sample(n, block_size, replace=FALSE)
V[b,] <- colSums(Z[resample_indexes,]) / block_size
}
# Trim the extreme tails
tl <- which(taus >= trim[1])[1]
th <- which(taus <= trim[2])[length(which(taus <= trim[2]))]
s <- block_size ^ (1 / 2) * apply(t(t(abs(V[,(tl:th)])) / sd_z[tl:th]), 1, max)
s <- as.vector(s)
critical_value <- quantile(s,1 - size)
process <- rbind(coef$endg_var,coef$exog_var)[variable,(tl:th)]
process <- process / sd_z[tl:th]
ks_stat <- max(abs(process)) * n ^ (1/2)
ks <- list()
class(ks) <- "ivqr_ks"
ks$ks_stat <- ks_stat
ks$critical_value <- critical_value
ks$block_size <- block_size
ks$B <- B
ks$s <- s
return(ks)
}
ivqr.ks.dom <- function(object, trim, B, variable, size) {
taus <- object$taus
data <- object$data
coef <- object$coef
fitted <- object$fitted
residuals <- object$residuals
endg_var_se <- object$se[variable]
dim_dk <- object$dim_d_d_k[1] + object$dim_d_d_k[3]
n <- object$n
PSI <- object$PSI
J <- object$vc$J
# nrow = n, ncol = length(taus)
L <- matrix(rep(taus,n), n, length(taus), byrow=TRUE) - as.numeric(residuals < 0)
Z <- matrix(NA,n,length(taus))
for (tau_index in 1:length(taus)) {
invJ <- solve(J[,,tau_index])
R_invJ <- invJ[variable,] # = R %*% invJ, R = c(0,0,,,0,1,0,..0)
Z[,tau_index] <- as.vector(L[,tau_index]) * (PSI %*% R_invJ)
}
sd_z <- colMeans(Z ^ 2) ^ (1 / 2)
block_size <- as.integer(5 * n ^ (2/5))
V <- matrix(NA, B, length(taus))
for (b in 1:B) {
resample_indexes <- sample(n, block_size, replace=FALSE)
V[b,] <- colSums(Z[resample_indexes,]) / block_size
}
# Trim the extreme tails
tl <- which(taus >= trim[1])[1]
th <- which(taus <= trim[2])[length(which(taus < trim[2]))]
s <- block_size ^ (1 / 2) * apply(t(t(abs(V[,(tl:th)])) / sd_z[tl:th]), 1, max)
s <- as.vector(s)
critical_value <- quantile(s,1 - size)
process <- rbind(coef$endg_var,coef$exog_var)[variable,(tl:th)]
ks_stat <- max(max(-1 * process / sd_z[tl:th]), 0) * n ^ (1/2)
ks <- list()
class(ks) <- "ivqr_ks"
ks$ks_stat <- ks_stat
ks$critical_value <- critical_value
ks$block_size <- block_size
ks$B <- B
ks$s <- s
return(ks)
}
ivqr.ks.const <- function(object, trim, B, variable, size) {
taus <- object$taus
data <- object$copy_data
coef <- object$coef
grid <- object$grid
fitted <- object$fitted
gridMethod <- object$gridMethod
ivqrMethod <- object$ivqrMethod
qrMethod <- object$qrMethod
residuals <- object$residuals
formula <- formula(object$formula)
# endg_var_se <- object$se[variable]
dim_dk <- object$dim_d_d_k[1] + object$dim_d_d_k[3]
n <- object$n
PSI <- object$PSI
J <- object$vc$J
ivqr_median_fit <- ivqr(formula, taus = 0.5, data = data , grid = grid,
gridMethod = gridMethod, ivqrMethod = ivqrMethod, qrMethod = qrMethod)
coef_m <- ivqr_median_fit$coef
coef_m <- rbind(coef_m$endg_var, coef_m$exog_var)[variable,1]
residuals_m <- ivqr_median_fit$residuals[,1]
J_m <- ivqr_median_fit$vc$J[,,1]
L_m <- matrix(rep(0.5,n), n, length(0.5), byrow=TRUE) - as.numeric(residuals_m < 0)
# nrow = n, ncol = length(taus)
L <- matrix(rep(taus,n), n, length(taus), byrow=TRUE) - as.numeric(residuals < 0)
Z <- matrix(NA,n,length(taus))
for (tau_index in 1:length(taus)) {
invJ <- solve(J[,,tau_index])
invJ_m <- solve(J_m)
R_invJ_m <- invJ_m[variable,]
R_invJ <- invJ[variable,] # = R %*% invJ, R = c(0,0,,,0,1,0,..0)
Z[,tau_index] <- (as.vector(L[,tau_index]) * (PSI %*% R_invJ)
- as.vector(L_m[,1]) * (PSI %*% R_invJ_m))
}
sd_z <- colMeans(Z ^ 2) ^ (1 / 2)
block_size <- as.integer(5 * n ^ (2/5))
V <- matrix(NA, B, length(taus))
for (b in 1:B) {
resample_indexes <- sample(n, block_size, replace=FALSE)
V[b,] <- colSums(Z[resample_indexes,]) / block_size
}
# Trim the extreme tails cut out [1/2 - eps, 1/2 + eps]
tl <- which(taus >= trim[1])[1]
th <- which(taus <= trim[2])[length(which(taus < trim[2]))]
left_to_median <- which(taus < 0.5)[length(which(taus < 0.5 ))]
right_to_median <- which(taus > 0.5)[1]
indexes <- c(tl:left_to_median,right_to_median:th)
s <- (block_size) ^ (1 / 2) * apply(t(t(abs(V[,indexes])) / sd_z[indexes]),1,max)
s <- as.vector(s)
critical_value <- quantile(s,1 - size)
process <- rbind(coef$endg_var,coef$exog_var)[variable,indexes]
ks_stat <- n ^ (1/2) * max(abs(process - coef_m) / sd_z[indexes])
# print(process)
# print(coef_m)
# print(sd_z[indexes])
# print(ks_stat)
# stop()
ks <- list()
class(ks) <- "ivqr_ks"
ks$ks_stat <- ks_stat
ks$critical_value <- critical_value
ks$block_size <- block_size
ks$B <- B
ks$s <- s
return(ks)
}
ivqr.ks.exog <- function(object, trim, B, variable, size, bd_rule = "Silver") {
taus <- object$taus
data <- object$copy_data
coef <- object$coef
grid <- object$grid
fitted <- object$fitted
gridMethod <- object$gridMethod
ivqrMethod <- object$ivqrMethod
qrMethod <- object$qrMethod
residuals <- object$residuals
formula <- formula(object$formula)
dim_dk <- object$dim_d_d_k[1] + object$dim_d_d_k[3]
dim_d <- object$dim_d_d_k[1]
n <- object$n
PSI <- object$PSI
J <- object$vc$J
formula_rq <- formula(Formula(formula), lhs = 1, rhs = c(1,3), collapse = TRUE)
fit_rq <- withCallingHandlers(
quantreg::rq(formula_rq, tau = taus, data = data, method = qrMethod),
warning = Suppress_Sol_not_Unique
)
coef_rq <- fit_rq$coef
rq_residuals <- fit_rq$residuals
model_matrix <- model.matrix(formula_rq,data)
X <- cbind(model_matrix[,1],model_matrix[,(2 + dim_d):(dim(model_matrix)[2])])
D <- model_matrix[,2:(1 + dim_d)]
X_tilde <- cbind(D,X)
# nrow = n, ncol = length(taus)
A <- matrix(rep(taus,n), n, length(taus), byrow=TRUE) - as.numeric(rq_residuals < 0)
L <- matrix(rep(taus,n), n, length(taus), byrow=TRUE) - as.numeric(residuals < 0)
Z <- matrix(NA,n,length(taus))
for (tau_index in 1:length(taus)) {
e <- rq_residuals[,tau_index]
if ( bd_rule == "Silver" ) {
h <- 1.364 * ( (2*sqrt(pi)) ^ (-1/5) ) * sd(e) * ( n ^ (-1/5) )
}
kernel <- c(as.numeric( abs(e) < h ))
H <- (1 / (2 * n * h)) * t(kernel * X_tilde) %*% X_tilde
invH <- solve(H)
invJ <- solve(J[,,tau_index])
R_invH <- invH[variable,]
R_invJ <- invJ[variable,] # = R %*% invJ, R = c(0,0,,,0,1,0,..0)
Z[,tau_index] <- (as.vector(L[,tau_index]) * (PSI %*% R_invJ)
- as.vector(A[,tau_index]) * (X_tilde %*% R_invH))
}
sd_z <- colMeans(Z ^ 2) ^ (1 / 2)
block_size <- as.integer(5 * n ^ (2/5))
V <- matrix(NA, B, length(taus))
for (b in 1:B) {
resample_indexes <- sample(n, block_size, replace=FALSE)
V[b,] <- colSums(Z[resample_indexes,]) / block_size
}
# Trim the extreme tails cut out [1/2 - eps, 1/2 + eps]
tl <- which(taus >= trim[1])[1]
th <- which(taus <= trim[2])[length(which(taus < trim[2]))]
s <- block_size ^ (1 / 2) * apply(t(t(abs(V[,(tl:th)])) / sd_z[tl:th]), 1, max)
s <- as.vector(s)
critical_value <- quantile(s,1 - size)
process <- rbind(
coef$endg_var - coef_rq[2:(1 + dim_d),],
coef$exog_var[1,] - coef_rq[1,],
coef$exog_var[2:dim(coef$exog_var)[1],] - coef_rq[(2 + dim_d):(dim(coef_rq))[1],]
)[variable,(tl:th)]
ks_stat <- n ^ (1/2) * max(abs(process) / sd_z[(tl:th)])
ks <- list()
class(ks) <- "ivqr_ks"
ks$ks_stat <- ks_stat
ks$critical_value <- critical_value
ks$block_size <- block_size
ks$B <- B
ks$s <- s
return(ks)
}
#' @export
print.ivqr_ks <- function(x, ...) {
print(x$ks_stat)
print(x$critical_value)
print(paste("Block size:",x$block_size))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.