R/paf.R

In DER: Income Polarization Index

paf <- function(y, a) {
  y <- y / mean(y)
  n <- length(y)
  d <- Rfast::vecdist(y)
  if ( length(a) == 1 ) {
    h <- 4.7 / sqrt(n) * sd(y) * a^0.1  ## bandwidth
    fhat <- Rfast::rowmeans( exp( -0.5 * d^2 / h^2 ) ) / sqrt(2 * pi) / h
    fhata <- fhat^a
    paf <- sum( fhata * d ) / n^2
    alien <- mean(d)
    ident <- mean(fhata)
    rho <- paf / (alien * ident) - 1
    res <- c(paf, alien, ident, 1 + rho)
    names(res) <- c("paf", "alienation", "identification", "1 + rho")
  } else {
    com <- 4.7 / sqrt(n) * sd(y)
    d2 <-  -0.5 * d^2
    alien <- mean(d)
    lena <- length(a)
    paf <- ident <- numeric(lena)
    for ( i in 1:lena ) {
      h <- com * a[i]^0.1  ## bandwidth
      fhat <- Rfast::rowmeans( exp( d2 / h^2 ) ) / sqrt(2 * pi) / h
      fhata <- fhat^a[i]
      paf[i] <- sum( fhata * d ) / n^2
      ident[i] <- mean(fhata)
    }
    rho <- paf / ( alien * ident ) - 1
    res <- cbind(paf, alien, ident, 1 + rho)
    colnames(res) <- c("paf", "alienation", "identification", "1 + rho")
    rownames(res) <- paste( "alpha=", a, sep = "" )
  }
  res
}