Welfare Test Resource Equivalent Variation (REV) Computation (Line by Line) • Optimal Allocation

Line by line, calculate welfare distances, Resource Equivalent Variation (REV), for the Binary Allocation problem based on Theorem 1. Hand-input A and alpha values.

Load Dependencies

rm(list = ls(all.names = TRUE))
options(knitr.duplicate.label = 'allow')

library(dplyr)
library(tidyr)
library(tibble)

library(PrjOptiAlloc)

library(knitr)
library(kableExtra)

Implement Algorithm Line by Line

Generate Testing Rank Data

Generate data needed.

# 1. Generate A and alpha Testing Vectors, random list no meanings
ar_alpha <- c(0.1,1.5,2.5,4  ,9,1.2,3,2,8)
ar_A <-     c(0.5,1.5,3.5,6.5,1.9,3,4,6,4)
ar_beta <-  rep(0, length(ar_A)) + 1/length(ar_A)
mt_alpha_A <- cbind(ar_alpha, ar_A, ar_beta)
ar_st_varnames <- c('alpha', 'A', 'beta')
tb_alpha_A <- as_tibble(mt_alpha_A) %>% rename_all(~c(ar_st_varnames))
tb_alpha_A

## # A tibble: 9 x 3
##   alpha     A  beta
##   <dbl> <dbl> <dbl>
## 1   0.1   0.5 0.111
## 2   1.5   1.5 0.111
## 3   2.5   3.5 0.111
## 4   4     6.5 0.111
## 5   9     1.9 0.111
## 6   1.2   3   0.111
## 7   3     4   0.111
## 8   2     6   0.111
## 9   8     4   0.111

ar_rho <- c(0.99)
for (it_rho_ctr in seq(1,length(ar_rho))) {
  fl_rho <- ar_rho[it_rho_ctr]

  ls_ranks <- ffp_opt_sobin_target_row(tb_alpha_A[1,], fl_rho, ar_A, ar_alpha, ar_beta)
  it_rank <- ls_ranks$it_rank
  ar_it_rank <- ls_ranks$ar_it_rank
  ar_fl_rank_val <- ls_ranks$ar_fl_rank_val

  cat('fl_rho:', fl_rho, 'it_rank:', it_rank, '\n')
  cat('ar_it_rank:', ar_it_rank, '\n')
  cat('ar_fl_rank_val:', ar_fl_rank_val, '\n')

}

## fl_rho: 0.99 it_rank: 9 
## ar_it_rank: 9 7 5 3 1 8 4 6 2 
## ar_fl_rank_val: 1 14.79282 24.46935 38.92284 87.90466 11.77704 29.32046 19.49807 77.92111

Compute REV

Compute REV giving observable and ranking

# Use the just generated optimal ranking above
ar_it_obs <- c(0,1,0,1,0,0,1,1,0)
ar_opti_r <- ar_it_rank
mt_rev <- cbind(seq(1,length(ar_opti_r)), ar_it_obs, ar_opti_r, ar_A, ar_alpha, ar_beta)
ar_st_varnames <- c('id', 'observed', 'optimal', 'A', 'alpha', 'beta')
tb_onevar <- as_tibble(mt_rev) %>% rename_all(~c(ar_st_varnames)) %>%
                arrange(optimal)

## Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if `.name_repair` is omitted as of tibble 2.0.0.
## Using compatibility `.name_repair`.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.

it_w_res <- sum(ar_it_obs)
fl_rho <- ar_rho

# Generate util observed and util unobserved columns
tb_onevar <- tb_onevar %>% mutate(utility = beta*((A+alpha)^fl_rho - A^fl_rho)) %>%
                           mutate(util_ob = case_when(observed == 1 ~ utility,
                                                      TRUE ~ 0 )) %>%
                           mutate(util_notob = case_when(observed == 0 ~ utility,
                                                         TRUE ~ 0 ))

# Display
kable(tb_onevar) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

id	observed	optimal	A	alpha	beta	utility	util_ob	util_notob
5	0	1	1.9	9.0	0.1111111	0.9727629	0.0000000	0.9727629
9	0	2	4.0	8.0	0.1111111	0.8622839	0.0000000	0.8622839
4	1	3	6.5	4.0	0.1111111	0.4307246	0.4307246	0.0000000
7	1	4	4.0	3.0	0.1111111	0.3244636	0.3244636	0.0000000
3	0	5	3.5	2.5	0.1111111	0.2707806	0.0000000	0.2707806
8	1	6	6.0	2.0	0.1111111	0.2157678	0.2157678	0.0000000
2	1	7	1.5	1.5	0.1111111	0.1636991	0.1636991	0.0000000
6	0	8	3.0	1.2	0.1111111	0.1303261	0.0000000	0.1303261
1	0	9	0.5	0.1	0.1111111	0.0110661	0.0000000	0.0110661

Generate Util observed um and not observed Reverse Sum

# Generate directional cumulative sums
tb_onevar <- tb_onevar %>%
                arrange(-optimal) %>%
                mutate(util_ob_sum = cumsum(util_ob)) %>%
                arrange(optimal) %>%
                mutate(util_notob_sum = cumsum(util_notob)) %>%
                select(id, observed, optimal, alpha, A, beta,
                       utility,
                       util_ob, util_ob_sum,
                       util_notob, util_notob_sum)

# Display
kable(tb_onevar) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

id	observed	optimal	alpha	A	beta	utility	util_ob	util_ob_sum	util_notob	util_notob_sum
5	0	1	9.0	1.9	0.1111111	0.9727629	0.0000000	1.1346551	0.9727629	0.9727629
9	0	2	8.0	4.0	0.1111111	0.8622839	0.0000000	1.1346551	0.8622839	1.8350468
4	1	3	4.0	6.5	0.1111111	0.4307246	0.4307246	1.1346551	0.0000000	1.8350468
7	1	4	3.0	4.0	0.1111111	0.3244636	0.3244636	0.7039305	0.0000000	1.8350468
3	0	5	2.5	3.5	0.1111111	0.2707806	0.0000000	0.3794669	0.2707806	2.1058273
8	1	6	2.0	6.0	0.1111111	0.2157678	0.2157678	0.3794669	0.0000000	2.1058273
2	1	7	1.5	1.5	0.1111111	0.1636991	0.1636991	0.1636991	0.0000000	2.1058273
6	0	8	1.2	3.0	0.1111111	0.1303261	0.0000000	0.0000000	0.1303261	2.2361534
1	0	9	0.1	0.5	0.1111111	0.0110661	0.0000000	0.0000000	0.0110661	2.2472195

# Following Theorem 1, simply compare util_ob_sum vs util_notob_sum
tb_onevar <- tb_onevar %>%
      mutate(opti_better_obs =
               case_when(util_notob_sum/util_ob_sum > 1 ~ 1,
                         TRUE ~ 0)
             ) %>%
      select(id, observed, optimal,
             utility,
             util_ob, util_ob_sum,
             util_notob, util_notob_sum,
             opti_better_obs,
             alpha, A, beta)

# Display
kable(tb_onevar) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

id	observed	optimal	utility	util_ob	util_ob_sum	util_notob	util_notob_sum	opti_better_obs	alpha	A	beta
5	0	1	0.9727629	0.0000000	1.1346551	0.9727629	0.9727629	0	9.0	1.9	0.1111111
9	0	2	0.8622839	0.0000000	1.1346551	0.8622839	1.8350468	1	8.0	4.0	0.1111111
4	1	3	0.4307246	0.4307246	1.1346551	0.0000000	1.8350468	1	4.0	6.5	0.1111111
7	1	4	0.3244636	0.3244636	0.7039305	0.0000000	1.8350468	1	3.0	4.0	0.1111111
3	0	5	0.2707806	0.0000000	0.3794669	0.2707806	2.1058273	1	2.5	3.5	0.1111111
8	1	6	0.2157678	0.2157678	0.3794669	0.0000000	2.1058273	1	2.0	6.0	0.1111111
2	1	7	0.1636991	0.1636991	0.1636991	0.0000000	2.1058273	1	1.5	1.5	0.1111111
6	0	8	0.1303261	0.0000000	0.0000000	0.1303261	2.2361534	1	1.2	3.0	0.1111111
1	0	9	0.0110661	0.0000000	0.0000000	0.0110661	2.2472195	1	0.1	0.5	0.1111111

# Find the Rank number tha tmatches the First 1 in opti_better_obs
tb_onevar <- tb_onevar %>%
  mutate(opti_better_obs_cumu = cumsum(opti_better_obs)) %>%
  mutate(min_w_hat =
           case_when(opti_better_obs_cumu == 1 ~ optimal,
                     TRUE ~ 0)) %>%
      select(id, observed, optimal,
             utility,
             util_ob, util_ob_sum,
             util_notob, util_notob_sum,
             opti_better_obs,
             min_w_hat,
             alpha, A, beta)

# Display
  kable(tb_onevar) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

id	observed	optimal	utility	util_ob	util_ob_sum	util_notob	util_notob_sum	opti_better_obs	min_w_hat	alpha	A	beta
5	0	1	0.9727629	0.0000000	1.1346551	0.9727629	0.9727629	0	0	9.0	1.9	0.1111111
9	0	2	0.8622839	0.0000000	1.1346551	0.8622839	1.8350468	1	2	8.0	4.0	0.1111111
4	1	3	0.4307246	0.4307246	1.1346551	0.0000000	1.8350468	1	0	4.0	6.5	0.1111111
7	1	4	0.3244636	0.3244636	0.7039305	0.0000000	1.8350468	1	0	3.0	4.0	0.1111111
3	0	5	0.2707806	0.0000000	0.3794669	0.2707806	2.1058273	1	0	2.5	3.5	0.1111111
8	1	6	0.2157678	0.2157678	0.3794669	0.0000000	2.1058273	1	0	2.0	6.0	0.1111111
2	1	7	0.1636991	0.1636991	0.1636991	0.0000000	2.1058273	1	0	1.5	1.5	0.1111111
6	0	8	0.1303261	0.0000000	0.0000000	0.1303261	2.2361534	1	0	1.2	3.0	0.1111111
1	0	9	0.0110661	0.0000000	0.0000000	0.0110661	2.2472195	1	0	0.1	0.5	0.1111111

# Compute REV
it_min_w_hat <- max(tb_onevar %>% pull(min_w_hat))
it_w_hat_obs <- sum(ar_it_obs)
delta_rev <- 1-(it_min_w_hat/it_w_hat_obs)

cat('it_min_w_hat:', it_min_w_hat, '\n')

## it_min_w_hat: 2

cat('it_w_hat_obs:', it_w_hat_obs, '\n')

## it_w_hat_obs: 4

cat('delta_rev:', delta_rev, '\n')

## delta_rev: 0.5

Function test

Test the function written based on the details above.

# ar_alpha <- c(0.1,1.5,2.5,4  ,9,1.2,3,2,8)
# ar_A <-     c(0.5,1.5,3.5,6.5,1.9,3,4,6,4)
# ar_beta <-  rep(0, length(ar_A)) + 1/length(ar_A)
# mt_alpha_A <- cbind(ar_alpha, ar_A, ar_beta)
# ar_st_varnames <- c('alpha', 'A', 'beta')
# tb_alpha_A <- as_tibble(mt_alpha_A) %>% rename_all(~c(ar_st_varnames))
# tb_alpha_A
# ar_rho <- c(-1)
for (it_rho_ctr in seq(1,length(ar_rho))) {
  fl_rho <- ar_rho[it_rho_ctr]

  ls_ranks <- ffp_opt_sobin_target_row(tb_alpha_A[1,], fl_rho, ar_A, ar_alpha, ar_beta)
  it_rank <- ls_ranks$it_rank
  ar_it_rank <- ls_ranks$ar_it_rank
  ar_fl_rank_val <- ls_ranks$ar_fl_rank_val

  cat('fl_rho:', fl_rho, 'it_rank:', it_rank, '\n')
  cat('ar_it_rank:', ar_it_rank, '\n')
  cat('ar_fl_rank_val:', ar_fl_rank_val, '\n')

}

## fl_rho: 0.99 it_rank: 9 
## ar_it_rank: 9 7 5 3 1 8 4 6 2 
## ar_fl_rank_val: 1 14.79282 24.46935 38.92284 87.90466 11.77704 29.32046 19.49807 77.92111

ar_bin_observed <- ar_it_obs
ar_queue_optimal <- ar_it_rank
ffp_opt_sobin_rev(ar_queue_optimal, ar_bin_observed, ar_A, ar_alpha, ar_beta, fl_rho)

## [1] 0.5