R/paf2.boot.R

In DER: Income Polarization Index

paf2.boot <- function(y, a, R = 1000) {
  index <- DER::paf2(y, a)
  boot <- matrix(0, R, 3)
  colnames(boot) <- c("paf", "deprivation", "surplus")
  n <- length(y)
  Y <- y

  for (i in 1:R) {
    ind <- Rfast2::Sample.int(n, n, replace = TRUE)
    y <- Y[ind]
    y <- y / mean(y)
    h <- 4.7 / sqrt(n) * sd(y) * a^0.1  ## bandwidth
    dD <- dS <- outer(y, y, "-")
    fhat <- Rfast::rowmeans( exp( -0.5 * dD^2 / h^2 ) ) / sqrt(2 * pi) / h
    fhata <- fhat^a
    dD[dD > 0] <- 0
    dS[dS < 0] <- 0
    D <- sum( fhata * abs(dD) ) / n^2
    S <- sum( fhata * dS ) / n^2
    boot[i, ] <- c(D + S, D, S)
  }

  mesoi <- Rfast::colmeans(boot)
  bias <- index - mesoi
  se <- Rfast::colVars(boot, std = TRUE)
  ci <- Rfast2::colQuantile( boot, probs = c(0.025, 0.975) )
  info <- rbind(mesoi, bias, se, ci)
  colnames(info) <- colnames(boot)
  rownames(info) <- c("mesoi", "bias", "se", "2.5%", "97.5%" )
  list(boot = boot, index = index, info = info)
}