R/bpbounds_biv_x2y2z3.R

In bpbounds: Nonparametric Bounds for the Average Causal Effect Due to Balke and Pearl and Extensions

A_biv_x2y2z3 <- matrix(nrow = 60, ncol = 12)
A_biv_x2y2z3[1, ] <- c(0, -1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0)
A_biv_x2y2z3[2, ] <- c(0, 0, 0, -1, 0, 0, 0, 0, 0, -1, 0, 0)
A_biv_x2y2z3[3, ] <- c(0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, -1)
A_biv_x2y2z3[4, ] <- c(0, -1, 0, 0, 0, 1, 0, -1, 0, 0, 0, -1)
A_biv_x2y2z3[5, ] <- c(0, -1, 0, 1, 0, 0, 0, -1, 0, -1, 0, 0)
A_biv_x2y2z3[6, ] <- c(0, 0, 0, -1, 0, 1, 0, 0, 0, -1, 0, -1)
A_biv_x2y2z3[7, ] <- c(0, 0, 0, 1, 0, -1, 0, 0, 0, -1, 0, -1)
A_biv_x2y2z3[8, ] <- c(0, 1, 0, -1, 0, 0, 0, -1, 0, -1, 0, 0)
A_biv_x2y2z3[9, ] <- c(0, 1, 0, 0, 0, -1, 0, -1, 0, 0, 0, -1)
A_biv_x2y2z3[10, ] <- c(0, -2, 0, 0, 0, 1, 0, -2, 0, 0, 0, 0)
A_biv_x2y2z3[11, ] <- c(0, -2, 0, 1, 0, 0, 0, -2, 0, 0, 0, 0)
A_biv_x2y2z3[12, ] <- c(0, 0, 0, -2, 0, 1, 0, 0, 0, -2, 0, 0)
A_biv_x2y2z3[13, ] <- c(0, 0, 0, 1, 0, -2, 0, 0, 0, 0, 0, -2)
A_biv_x2y2z3[14, ] <- c(0, 1, 0, -2, 0, 0, 0, 0, 0, -2, 0, 0)
A_biv_x2y2z3[15, ] <- c(0, 1, 0, 0, 0, -2, 0, 0, 0, 0, 0, -2)
A_biv_x2y2z3[16, ] <- c(0, -1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0)
A_biv_x2y2z3[17, ] <- c(0, 0, 0, -1, 0, 0, 0, 0, 0, 1, 0, 0)
A_biv_x2y2z3[18, ] <- c(0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 1)
A_biv_x2y2z3[19, ] <- c(0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, -1)
A_biv_x2y2z3[20, ] <- c(0, 0, 0, 1, 0, 0, 0, 0, 0, -1, 0, 0)
A_biv_x2y2z3[21, ] <- c(0, 1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0)
A_biv_x2y2z3[22, ] <- c(0, -1, 0, 0, 0, 1, 0, 1, 0, 0, 0, -1)
A_biv_x2y2z3[23, ] <- c(0, -1, 0, 1, 0, 0, 0, 1, 0, -1, 0, 0)
A_biv_x2y2z3[24, ] <- c(0, 0, 0, -1, 0, 1, 0, 0, 0, 1, 0, -1)
A_biv_x2y2z3[25, ] <- c(0, 0, 0, 1, 0, -1, 0, 0, 0, -1, 0, 1)
A_biv_x2y2z3[26, ] <- c(0, 1, 0, -1, 0, 0, 0, -1, 0, 1, 0, 0)
A_biv_x2y2z3[27, ] <- c(0, 1, 0, 0, 0, -1, 0, -1, 0, 0, 0, 1)
A_biv_x2y2z3[28, ] <- c(0, -1, 0, 0, 0, 1, 0, -1, 0, 0, 0, 1)
A_biv_x2y2z3[29, ] <- c(0, -1, 0, 1, 0, 0, 0, -1, 0, 1, 0, 0)
A_biv_x2y2z3[30, ] <- c(0, 0, 0, -1, 0, 1, 0, 0, 0, -1, 0, 1)
A_biv_x2y2z3[31, ] <- c(0, 0, 0, 1, 0, -1, 0, 0, 0, 1, 0, -1)
A_biv_x2y2z3[32, ] <- c(0, 1, 0, -1, 0, 0, 0, 1, 0, -1, 0, 0)
A_biv_x2y2z3[33, ] <- c(0, 1, 0, 0, 0, -1, 0, 1, 0, 0, 0, -1)
A_biv_x2y2z3[34, ] <- c(0, -1, 0, 0, 0, 2, 0, 0, 0, 0, 0, -2)
A_biv_x2y2z3[35, ] <- c(0, -1, 0, 2, 0, 0, 0, 0, 0, -2, 0, 0)
A_biv_x2y2z3[36, ] <- c(0, 0, 0, -1, 0, 2, 0, 0, 0, 0, 0, -2)
A_biv_x2y2z3[37, ] <- c(0, 0, 0, 2, 0, -1, 0, 0, 0, -2, 0, 0)
A_biv_x2y2z3[38, ] <- c(0, 2, 0, -1, 0, 0, 0, -2, 0, 0, 0, 0)
A_biv_x2y2z3[39, ] <- c(0, 2, 0, 0, 0, -1, 0, -2, 0, 0, 0, 0)
A_biv_x2y2z3[40, ] <- c(0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1)
A_biv_x2y2z3[41, ] <- c(0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0)
A_biv_x2y2z3[42, ] <- c(0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0)
A_biv_x2y2z3[43, ] <- c(0, -1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1)
A_biv_x2y2z3[44, ] <- c(0, -1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0)
A_biv_x2y2z3[45, ] <- c(0, 0, 0, -1, 0, 1, 0, 0, 0, 1, 0, 1)
A_biv_x2y2z3[46, ] <- c(0, 0, 0, 1, 0, -1, 0, 0, 0, 1, 0, 1)
A_biv_x2y2z3[47, ] <- c(0, 1, 0, -1, 0, 0, 0, 1, 0, 1, 0, 0)
A_biv_x2y2z3[48, ] <- c(0, 1, 0, 0, 0, -1, 0, 1, 0, 0, 0, 1)
A_biv_x2y2z3[49, ] <- c(0, -2, 0, 0, 0, 1, 0, 2, 0, 0, 0, 0)
A_biv_x2y2z3[50, ] <- c(0, -2, 0, 1, 0, 0, 0, 2, 0, 0, 0, 0)
A_biv_x2y2z3[51, ] <- c(0, 0, 0, -2, 0, 1, 0, 0, 0, 2, 0, 0)
A_biv_x2y2z3[52, ] <- c(0, 0, 0, 1, 0, -2, 0, 0, 0, 0, 0, 2)
A_biv_x2y2z3[53, ] <- c(0, 1, 0, -2, 0, 0, 0, 0, 0, 2, 0, 0)
A_biv_x2y2z3[54, ] <- c(0, 1, 0, 0, 0, -2, 0, 0, 0, 0, 0, 2)
A_biv_x2y2z3[55, ] <- c(0, -1, 0, 0, 0, 2, 0, 0, 0, 0, 0, 2)
A_biv_x2y2z3[56, ] <- c(0, -1, 0, 2, 0, 0, 0, 0, 0, 2, 0, 0)
A_biv_x2y2z3[57, ] <- c(0, 0, 0, -1, 0, 2, 0, 0, 0, 0, 0, 2)
A_biv_x2y2z3[58, ] <- c(0, 0, 0, 2, 0, -1, 0, 0, 0, 2, 0, 0)
A_biv_x2y2z3[59, ] <- c(0, 2, 0, -1, 0, 0, 0, 2, 0, 0, 0, 0)
A_biv_x2y2z3[60, ] <- c(0, 2, 0, 0, 0, -1, 0, 2, 0, 0, 0, 0)

alpha_biv_x2y2z3 <- c(-1, -1, -1, 0, 0, 0, 0, 0, 0,
                      -1, -1, -1, -1, -1, -1, 1, 1, 1,
                      1, 1, 1, 1, 1, 1, 1, 1, 1, -1,
                      -1, -1, -1, -1, -1, 1, 1, 1, 1,
                      1, 1, -1, -1, -1, 0, 0, 0, 0, 0,
                      0, 1, 1, 1, 1, 1, 1, -1, -1, -1,
                      -1, -1, -1)

cons_biv_x2y2z3 <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                     0, 0, 0, 0, 0, -1, -1, -1, -1, -1,
                     -1, -1, -1, -1, -1, -1, -1, -1, -1,
                     -1, -1, -1, -1, -1, -1, -1, -1, -1,
                     -1, -2, -2, -2, -2, -2, -2, -2, -2,
                     -2, -2, -2, -2, -2, -2, -2, -3, -3,
                     -3, -3, -3, -3)

bpbounds_biv_x2y2z3 <- function(p,
                                A = A_biv_x2y2z3,
                                alpha = alpha_biv_x2y2z3,
                                cons = cons_biv_x2y2z3) {
  prod <- A %*% p + cons
  ivinequality <- prod[alpha == 0]
  low <- prod[alpha == -1]
  upp = -1 * prod[alpha == 1]
  inequality <- max(ivinequality) <= 0
  bplow <- max(low)
  bpupp <- min(upp)

  return(
    list(
      "inequality" = inequality,
      "bplow" = bplow,
      "bpupp" = bpupp,
      "bplower" = low,
      "bpupper" = upp
    )
  )
}