R/wald_test.R

In sdpdth: M-Estimator for Threshold Spatial Dynamic Panel Data Model

#' Wald test for threshold spatial dynamic panel data model
#' 
#' @importFrom stats pnorm
#' 
#' @param th_res msdpdth class, generated by function msdpdth()
#' 
#' 
#' @details Two sided Wald test for testing whether two estimated parameters for each group are equal
#' \itemize{
#' \item{"h_0"} {\eqn{\theta_1 = \theta_2}}
#' \item{"h_1"} {\eqn{\theta_1 \neq \theta_2}}
#' }
#' 
#' @return A list of p-values for each parameter.
#' 
#' @examples
#' \donttest{
#' data(data_th, data_w)
#' result <- msdpdth(y = data_th$y, x = data_th$x, w1 = data_w, th = data_th$th)
#' wald_test(result)
#' }
#' 
#' @export
#' 
#' 

wald_test = function(th_res){
  tt_p_value = function(x){
    temp_vec = x > 0.5
    return(2*abs(temp_vec - x))
  }
  output = list()
  k = length(th_res$coefficient$beta)
  var_vec = numeric(k/2)
  for (i in 1:(k/2)){
    var_vec[i] = th_res$vc_mat[i,i] + th_res$vc_mat[i+k/2,i+k/2] - 2*th_res$vc_mat[i, i+k/2]
  }
  output[["beta"]] = tt_p_value(pnorm(th_res$coefficient$beta[1:(k/2)] - th_res$coefficient$beta[(k/2+1):k], mean = rep(0,k/2), sd = sqrt(var_vec)))
  switch (th_res$model,
          "full" = {
            output[["lambda1"]] = tt_p_value(pnorm(th_res$coefficient$lambda1_1 - th_res$coefficient$lambda1_2, mean = 0, sd = sqrt(th_res$coefficient$lambda1_1_se^2 + th_res$coefficient$lambda1_2_se^2 - 2*th_res$vc_mat[k+3, k+7])))
            output[["lambda2"]] = tt_p_value(pnorm(th_res$coefficient$lambda2_1 - th_res$coefficient$lambda2_2, mean = 0, sd = sqrt(th_res$coefficient$lambda2_1_se^2 + th_res$coefficient$lambda2_2_se^2 - 2*th_res$vc_mat[k+4, k+8])))
            output[["lambda3"]] = tt_p_value(pnorm(th_res$coefficient$lambda3_1 - th_res$coefficient$lambda3_2, mean = 0, sd = sqrt(th_res$coefficient$lambda3_1_se^2 + th_res$coefficient$lambda3_2_se^2 - 2*th_res$vc_mat[k+5, k+9])))
            output[["rho"]] = tt_p_value(pnorm(th_res$coefficient$rho_1 - th_res$coefficient$rho_2, mean = 0, sd = sqrt(th_res$coefficient$rho_1_se^2 + th_res$coefficient$rho_2_se^2 - 2*th_res$vc_mat[k+2, k+6])))
          },
          "slm" = {
            output[["lambda1"]] = tt_p_value(pnorm(th_res$coefficient$lambda1_1 - th_res$coefficient$lambda1_2, mean = 0, sd = sqrt(th_res$coefficient$lambda1_1_se^2 + th_res$coefficient$lambda1_2_se^2 - 2*th_res$vc_mat[k+3, k+5])))
            output[["rho"]] = tt_p_value(pnorm(th_res$coefficient$rho_1 - th_res$coefficient$rho_2, mean = 0, sd = sqrt(th_res$coefficient$rho_1_se^2 + th_res$coefficient$rho_2_se^2 - 2*th_res$vc_mat[k+2, k+4])))
          },
          "sem" = {
            output[["lambda3"]] = tt_p_value(pnorm(th_res$coefficient$lambda3_1 - th_res$coefficient$lambda3_2, mean = 0, sd = sqrt(th_res$coefficient$lambda3_1_se^2 + th_res$coefficient$lambda3_2_se^2 - 2*th_res$vc_mat[k+3, k+5])))
            output[["rho"]] = tt_p_value(pnorm(th_res$coefficient$rho_1 - th_res$coefficient$rho_2, mean = 0, sd = sqrt(th_res$coefficient$rho_1_se^2 + th_res$coefficient$rho_2_se^2 - 2*th_res$vc_mat[k+2, k+4])))
          },
          "sltl" = {
            output[["lambda1"]] = tt_p_value(pnorm(th_res$coefficient$lambda1_1 - th_res$coefficient$lambda1_2, mean = 0, sd = sqrt(th_res$coefficient$lambda1_1_se^2 + th_res$coefficient$lambda1_2_se^2 - 2*th_res$vc_mat[k+3, k+6])))
            output[["lambda2"]] = tt_p_value(pnorm(th_res$coefficient$lambda2_1 - th_res$coefficient$lambda2_2, mean = 0, sd = sqrt(th_res$coefficient$lambda1_1_se^2 + th_res$coefficient$lambda2_2_se^2 - 2*th_res$vc_mat[k+4, k+7])))
            output[["rho"]] = tt_p_value(pnorm(th_res$coefficient$rho_1 - th_res$coefficient$rho_2, mean = 0, sd = sqrt(th_res$coefficient$rho_1_se^2 + th_res$coefficient$rho_2_se^2 - 2*th_res$vc_mat[k+2, k+5])))
          },
          stop("undefined model")
  )
  return(output)
}