vignettes/ffv_opt_n2_fn_welfare.Rmd
ffv_opt_n2_fn_welfare.Rmd
spt_img_save <- '../_img/'
# spt_img_save_draft <- 'C:/Users/fan/Documents/Dropbox (UH-ECON)/repos/HgtOptiAlloDraft/_img/'
spt_img_save_draft <- 'G:/repos/HgtOptiAlloDraft/_img/'
# Draw data
fl_beta_1 = 0.5
fl_beta_2 = 1-fl_beta_1
fl_w_dollar = 100
it_w_units_avg = 5
# Specify more parameters
fl_dis_w <- (it_w_units_avg + 1)
fl_cts_w <- fl_dis_w # common price
ar_rho <- c(0.5, -5)
ar_rho <- c(0.99, -0.01, -100)
# ar_rho <- c(0.99, -0.01, -6)
# ar_rho <- 1 - (10^(c(seq(-2,2, length.out=20))))These parameters persist and are used throughout the rest of the file for allocation allocations. Parameters are more or less shared by discrete and continuous parameters.
it_rand_seed <- sample(10^7, 1, replace=T)
it_w_units_avg = 5
it_rand_seed <- 8135788
# it_w_units_avg = 5
# it_rand_seed <- 1234
# it_w_units_avg = 4
# it_rand_seed <- 449879812
it_N <- 2
ls_out_n2 <- suppressMessages(
ffd_draw_n_alpha(fl_w_dollar=fl_w_dollar,
it_w_units_avg=it_w_units_avg,
it_N = it_N, it_rand_seed=it_rand_seed))
# Translate to 2D budget indiff units
# fl_e = endowment points
fl_A_1 = ls_out_n2$ar_A_i[1]
fl_A_2 = ls_out_n2$ar_A_i[2]
# fl_alpha_1l = alpha_il
# fl_eh_alloc_il = outcome at alloations
ar_alpha_1l = ls_out_n2$ls_ar_alpha_il[[1]]$ar_levels
ar_alpha_2l = ls_out_n2$ls_ar_alpha_il[[2]]$ar_levels
ar_alpha_1l_zr = c(0,ls_out_n2$ls_ar_alpha_il[[1]]$ar_levels)
ar_alpha_2l_zr = c(0,ls_out_n2$ls_ar_alpha_il[[2]]$ar_levels)
ar_A_1l <- fl_A_1 + cumsum(head(ar_alpha_1l_zr, -1))
ar_A_2l <- fl_A_2 + cumsum(head(ar_alpha_2l_zr, -1))
it_D_1l <- length(ar_alpha_1l)
it_D_2l <- length(ar_alpha_2l)
# Max Discrete Choice After which budget outside of choice strike zone fully
it_w_max <- length(ar_alpha_1l) + length(ar_alpha_2l)
ar_w_solve_at_disc <- seq(1, it_w_max)
# print
print(paste0('it_rand_seed:',it_rand_seed))## [1] "it_rand_seed:8135788"## [1] "fl_A_1:1.21673425481374,fl_A_2:0.958163167726895"## [1] "ar_alpha_1l:1.10258826611402" "ar_alpha_1l:0.998256934864124"
## [3] "ar_alpha_1l:0.974318847935741" "ar_alpha_1l:0.832411817148815"
## [5] "ar_alpha_1l:0.809088128226687" "ar_alpha_1l:0.783336005710612"## [1] "ar_alpha_2l:0.862318390675336" "ar_alpha_2l:0.812810517766219"
## [3] "ar_alpha_2l:0.731285183772518" "ar_alpha_2l:0.680349434537137"
## [5] "ar_alpha_2l:0.579903139915456"Given \(A_i\) and \(\alpha_{il}\), solve for optimal discrete allocation
Build up discrete input frame
# ID and max Discrete Allocation
mt_i_D <- cbind(c(1,2), c(it_D_1l, it_D_2l))
colnames(mt_i_D) <- c('id_i', 'D_max_i')
tb_i_D <- as_tibble(mt_i_D)
# A_i and alpha_il as matrix
mt_A_alpha <- rbind(cbind(1, seq(1,it_D_1l), ar_alpha_1l, ar_A_1l, fl_beta_1),
cbind(2, seq(1,it_D_2l), ar_alpha_2l, ar_A_2l, fl_beta_2))
colnames(mt_A_alpha) <- c('id_i', 'D_il', 'alpha_il', 'A_il', 'beta_i')
# Combine to generate input_il matrix
df_input_il <- tb_i_D %>%
uncount(D_max_i) %>%
rowid_to_column(var = "id_il") %>%
left_join(tb_i_D, by='id_i') %>%
group_by(id_i) %>%
mutate(D_il = row_number()) %>%
left_join(as_tibble(mt_A_alpha), by=(c('id_i'='id_i', 'D_il'='D_il')))Solve for optimal choices.
ls_dis_solu <- suppressWarnings(suppressMessages(
ffp_opt_anlyz_rhgin_dis(ar_rho,
fl_dis_w,
df_input_il,
bl_df_alloc_il = TRUE,
bl_return_V = TRUE,
bl_return_allQ_V = TRUE,
bl_return_inner_V = TRUE)))
df_queue_il_long_n2 <-ls_dis_solu$df_queue_il_long
df_queue_il_wide_n2 <- ls_dis_solu$df_queue_il_wide
df_alloc_i_long_n2 <- ls_dis_solu$df_alloc_i_long
df_rho_gini_n2 <- ls_dis_solu$df_rho_gini
df_alloc_il_long_n2 <- ls_dis_solu$df_alloc_il_longSort Value along Resource Expansion Path:
df_queue_il_long_n2 %>%
arrange(rho_val, Q_il) %>%
select(rho_val, id_i, id_il, Q_il, D_Wbin_il, V_sum_l, V_inner_Q_il, V_star_Q_il)## # A tibble: 33 x 8
## rho_val id_i id_il Q_il D_Wbin_il V_sum_l V_inner_Q_il V_star_Q_il
## <dbl> <dbl> <int> <int> <dbl> <dbl> <dbl> <dbl>
## 1 -100 2 7 1 1 4.79e-27 1.51e- 9 1.23
## 2 -100 1 1 2 1 1.45e-37 4.79e-27 1.83
## 3 -100 2 8 3 1 4.46e-43 1.45e-37 2.34
## 4 -100 1 2 4 1 4.14e-53 4.46e-43 2.65
## 5 -100 2 9 5 1 1.01e-53 5.15e-53 3.33
## 6 -100 1 3 6 1 2.72e-64 1.01e-53 3.39
## 7 -100 2 10 7 0 1.02e-61 1.02e-61 4.07
## 8 -100 2 11 8 0 1.55e-67 2.72e-64 4.32
## 9 -100 1 4 9 0 5.44e-72 1.55e-67 4.66
## 10 -100 1 5 10 0 2.34e-78 1.55e-67 4.66
## # ... with 23 more rowsGenerate all outcome combinations, and expected outcome combinations for the two individuals at all combinations. If allocations are binary, then there are 4 possible combinations, but have to skip zero/zero.
Remember that the choice set is a box: (1) when resources availabe is 1, can either give \((1,0)\) or \((0,1)\) to the two individuals. (2) Given individual upper bound say 3 at most for individual 1, and 3 at most for individual 2, if we have 6 units of resources, to exhaust budget, only allocation that uses all budget is \((3,3)\). (3) If the resources available is equal to 3 units, can choose \((0,3)\), \((1,2)\), \((2,1)\), \((3,0)\). The box is such that when we drop negative slopes through it, at the bottom left and top right, there are limited choices. Note the box’s bounds are the allocation bounds/constraints.
First generate combination index frame:
# generate index frame
df_feasible_allocate <-
as_tibble(expand.grid(id_1l = seq(0, it_D_1l),
id_2l = seq(0, it_D_2l))) %>%
filter(!(id_1l == 0 & id_2l == 0)) %>%
rowid_to_column(var = "id_alloc") %>%
expand_grid(as_tibble(seq(1, 2)) %>% rename(id_i = value)) %>%
mutate(beta_i = 0.5) %>%
select(id_alloc, id_i, id_1l, id_2l, beta_i)
df_feasible_allocate$id_i = as.factor(df_feasible_allocate$id_i)
df_feasible_allocate$id_1l = as.factor(df_feasible_allocate$id_1l)
df_feasible_allocate$id_2l = as.factor(df_feasible_allocate$id_2l)
# print
head(df_feasible_allocate, 10)## # A tibble: 10 x 5
## id_alloc id_i id_1l id_2l beta_i
## <int> <fct> <fct> <fct> <dbl>
## 1 1 1 1 0 0.5
## 2 1 2 1 0 0.5
## 3 2 1 2 0 0.5
## 4 2 2 2 0 0.5
## 5 3 1 3 0 0.5
## 6 3 2 3 0 0.5
## 7 4 1 4 0 0.5
## 8 4 2 4 0 0.5
## 9 5 1 5 0 0.5
## 10 5 2 5 0 0.5Generate \(\alpha\) matrix for the two individuals to prep for merging:
# For individual 1
df_alpha_1 <- as_tibble(mt_A_alpha) %>% filter(id_i == 1) %>%
select(alpha_il, D_il) %>%
rename(id_1l = D_il) %>%
mutate(id_i = 1) %>%
mutate(alpha_il_cum = cumsum(alpha_il)) %>%
select(-alpha_il)
# For individual 2
df_alpha_2 <- as_tibble(mt_A_alpha) %>% filter(id_i == 2) %>%
select(alpha_il, D_il) %>%
rename(id_2l = D_il) %>%
mutate(id_i = 2) %>%
mutate(alpha_il_cum = cumsum(alpha_il)) %>%
select(-alpha_il)
# To factors
df_alpha_1$id_i = as.factor(df_alpha_1$id_i)
df_alpha_1$id_1l = as.factor(df_alpha_1$id_1l)
df_alpha_2$id_i = as.factor(df_alpha_2$id_i)
df_alpha_2$id_2l = as.factor(df_alpha_2$id_2l)
# Print
print(df_alpha_1)## # A tibble: 6 x 3
## id_1l id_i alpha_il_cum
## <fct> <fct> <dbl>
## 1 1 1 1.10
## 2 2 1 2.10
## 3 3 1 3.08
## 4 4 1 3.91
## 5 5 1 4.72
## 6 6 1 5.5
print(df_alpha_2)## # A tibble: 5 x 3
## id_2l id_i alpha_il_cum
## <fct> <fct> <dbl>
## 1 1 2 0.862
## 2 2 2 1.68
## 3 3 2 2.41
## 4 4 2 3.09
## 5 5 2 3.67Merging result for A_i_l0 and alpha_o_i:
df_input_ib_all <- df_feasible_allocate %>%
left_join(df_alpha_1, by=c('id_i'='id_i', 'id_1l'='id_1l')) %>%
left_join(df_alpha_2, by=c('id_i'='id_i', 'id_2l'='id_2l')) %>%
mutate(alpha_o_i = case_when(id_i == 1 ~ alpha_il_cum.x,
id_i == 2 ~ alpha_il_cum.y)) %>%
select(-alpha_il_cum.x, -alpha_il_cum.y) %>%
mutate(alpha_o_i = case_when(is.na(alpha_o_i) ~ 0,
TRUE ~ alpha_o_i)) %>%
mutate(A_i_l0 = case_when(id_i == 1 ~ fl_A_1,
id_i == 2 ~ fl_A_2)) %>%
select(id_alloc, id_i, id_1l, id_2l, A_i_l0, alpha_o_i, beta_i) %>%
mutate(w_agg = as.numeric(id_1l) + as.numeric(id_2l))
print(df_input_ib_all)## # A tibble: 82 x 8
## id_alloc id_i id_1l id_2l A_i_l0 alpha_o_i beta_i w_agg
## <int> <fct> <fct> <fct> <dbl> <dbl> <dbl> <dbl>
## 1 1 1 1 0 1.22 1.10 0.5 3
## 2 1 2 1 0 0.958 0 0.5 3
## 3 2 1 2 0 1.22 2.10 0.5 4
## 4 2 2 2 0 0.958 0 0.5 4
## 5 3 1 3 0 1.22 3.08 0.5 5
## 6 3 2 3 0 0.958 0 0.5 5
## 7 4 1 4 0 1.22 3.91 0.5 6
## 8 4 2 4 0 0.958 0 0.5 6
## 9 5 1 5 0 1.22 4.72 0.5 7
## 10 5 2 5 0 0.958 0 0.5 7
## # ... with 72 more rowsThere are only two individuals, pick any feasible point, and compare against optimal allocations, at one aggregate resource point. The welfare function only takes one set of allocations, and the optimal allocation results contain welfare along optimal allocation queues for all individuals.
In the code below, I loop over each possible set of observed allocation from df_input_ib_all (2 rows at a time, each row an individual, 2 row is a set of allocations). Within each loop, resource equivalent variation is computed comparing the welfare at this “observed allocations” set to welfare along the optimal allocation queue, and finding out, along the resource expansion path (allocation queue), how much less resources are needed to achieve the same welfare as achieved by the “observed allocations”.
# Optimal Allocation Results Along Queue with Value
df_queue_il_long_with_V <- df_queue_il_long_n2
# Get information and statistics
df_input_ib_agg_grp <- df_input_ib_all %>% filter(w_agg==11)
unique(df_input_ib_agg_grp %>% pull(id_alloc))## [1] 27 33 39## # A tibble: 2 x 8
## id_alloc id_i id_1l id_2l A_i_l0 alpha_o_i beta_i w_agg
## <int> <fct> <fct> <fct> <dbl> <dbl> <dbl> <dbl>
## 1 41 1 6 5 1.22 5.5 0.5 13
## 2 41 2 6 5 0.958 3.67 0.5 13
it_w_agg <- df_input_ib %>% slice(1L) %>% pull(w_agg)
# Call REV function
ffp_opt_anlyz_sodis_rev(
ar_rho[3],
it_w_agg,
df_input_ib,
df_queue_il_long_with_V %>% filter(rho_val == -100),
)## Warning in min(df_queue_il_with_V %>% filter(!!sym(svr_V_star_Q_il) >=
## fl_util_alter_alloc) %>% : no non-missing arguments to min; returning Inf## Warning: `cols` is now required when using unnest().
## Please use `cols = c(rev)`## Warning in min(df_queue_il_with_V %>% filter(!!sym(svr_V_star_Q_il) >=
## fl_util_alter_alloc) %>% : no non-missing arguments to min; returning Inf## Warning: `cols` is now required when using unnest().
## Please use `cols = c(fl_util_alter_alloc)`## # A tibble: 1 x 4
## rho rho_val REV AlterOutcome
## <int> <dbl> <dbl> <dbl>
## 1 1 -100 0 4.66Compute REV all all feasible allocation combinations, using (MxP by N) to (MxQ by N+Z-1):
# Temp function
ffi_opt_anlyz_sodis_rev_wrapper <- function(df_input_ib_sub){
it_w_agg <- df_input_ib_sub %>%slice(1L) %>% pull(w_agg)
df_rev <- ffp_opt_anlyz_sodis_rev(
ar_rho=ar_rho,
it_w_agg=it_w_agg,
df_input_ib=df_input_ib_sub,
df_queue_il_long_with_V=df_queue_il_long_n2,
)
return(df_rev)
}
# Function over Groups
df_rev_allrho <- suppressWarnings(
df_input_ib_all %>%
group_by(id_alloc) %>%
do(df_rev = ffi_opt_anlyz_sodis_rev_wrapper(.)) %>%
unnest() %>%
rowid_to_column(var = "id_alloc_rho")) %>%
left_join(df_input_ib_all %>%
group_by(id_alloc) %>% slice(1L) %>% ungroup() %>%
select(id_alloc, id_1l, id_2l, w_agg),
join='id_alloc') %>%
arrange(w_agg, rho, id_1l, id_2l)## Joining, by = "id_alloc"
# show results
head(df_rev_allrho, 50)## # A tibble: 50 x 9
## id_alloc_rho id_alloc rho rho_val REV AlterOutcome id_1l id_2l w_agg
## <int> <int> <int> <dbl> <dbl> <dbl> <fct> <fct> <dbl>
## 1 19 7 1 0.99 0.667 1.52 0 1 3
## 2 1 1 1 0.99 0.667 1.64 1 0 3
## 3 20 7 2 -0.01 0.667 1.49 0 1 3
## 4 2 1 2 -0.01 0.667 1.49 1 0 3
## 5 21 7 3 -100 0.667 1.23 0 1 3
## 6 3 1 3 -100 0.667 0.965 1 0 3
## 7 40 14 1 0.99 0.5 1.92 0 2 4
## 8 22 8 1 0.99 0.5 2.07 1 1 4
## 9 4 2 1 0.99 0.5 2.13 2 0 4
## 10 41 14 2 -0.01 0.5 1.79 0 2 4
## # ... with 40 more rowsREV grouped by Rho:
df <- df_rev_allrho
vars.group <- c('rho_val')
var.numeric <- 'REV'
str.stats.group <- 'allperc'
ar.perc <- c(0.01, 0.05, 0.10, 0.25, 0.5, 0.75, 0.9, 0.95, 0.99)
ls_summ_by_group <- REconTools::ff_summ_bygroup(df, vars.group, var.numeric, str.stats.group, ar.perc)## Warning: `funs()` was deprecated in dplyr 0.8.0.
## Please use a list of either functions or lambdas:
##
## # Simple named list:
## list(mean = mean, median = median)
##
## # Auto named with `tibble::lst()`:
## tibble::lst(mean, median)
##
## # Using lambdas
## list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.## Warning: attributes are not identical across measure variables;
## they will be dropped
df_table_grp_stats <- ls_summ_by_group$df_table_grp_stats
df_row_grp_stats <- ls_summ_by_group$df_row_grp_stats
df_overall_stats <- ls_summ_by_group$df_overall_stats
df_row_stats_all <- ls_summ_by_group$df_row_stats_all
print(df_table_grp_stats)## # A tibble: 3 x 20
## rho.val mean median sd IQR mad `1%` `5%` `10%` `25%` `50%` `75%`
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 -100 0.458 0.4 0.247 0.394 0.264 0.0333 0.0909 0.2 0.273 0.4 0.667
## 2 -0.01 0.346 0.286 0.155 0.278 0.132 0.159 0.167 0.182 0.222 0.286 0.5
## 3 0.99 0.297 0.25 0.132 0.133 0.101 0.117 0.154 0.167 0.2 0.25 0.333
## # ... with 8 more variables: 90% <dbl>, 95% <dbl>, 99% <dbl>, min <dbl>,
## # max <dbl>, first <dbl>, last <dbl>, n.distinct <int>REV grouped by Aggregate Resources:
df <- df_rev_allrho
vars.group <- c('w_agg')
var.numeric <- 'REV'
str.stats.group <- 'allperc'
ar.perc <- c(0.01, 0.05, 0.10, 0.25, 0.5, 0.75, 0.9, 0.95, 0.99)
ls_summ_by_group <- REconTools::ff_summ_bygroup(df, vars.group, var.numeric, str.stats.group, ar.perc)## Warning: attributes are not identical across measure variables;
## they will be dropped
df_table_grp_stats <- ls_summ_by_group$df_table_grp_stats
df_row_grp_stats <- ls_summ_by_group$df_row_grp_stats
df_overall_stats <- ls_summ_by_group$df_overall_stats
df_row_stats_all <- ls_summ_by_group$df_row_stats_all
print(df_table_grp_stats)## # A tibble: 11 x 20
## w.agg mean median sd IQR mad `1%` `5%` `10%` `25%` `50%`
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 3 0.667 0.667 0 0 0 0.667 0.667 0.667 0.667 0.667
## 2 4 0.556 0.5 0.110 0 0 0.5 0.5 0.5 0.5 0.5
## 3 5 0.5 0.4 0.160 0.2 0 0.4 0.4 0.4 0.4 0.4
## 4 6 0.422 0.333 0.188 0.167 0 0.19 0.283 0.333 0.333 0.333
## 5 7 0.413 0.286 0.195 0.25 0 0.286 0.286 0.286 0.286 0.286
## 6 8 0.347 0.25 0.180 0.125 0 0.146 0.231 0.25 0.25 0.25
## 7 9 0.274 0.222 0.132 0.0556 0 0.127 0.189 0.222 0.222 0.222
## 8 10 0.233 0.2 0.0888 0 0 0.2 0.2 0.2 0.2 0.2
## 9 11 0.202 0.182 0.0758 0 0 0.0982 0.127 0.164 0.182 0.182
## 10 12 0.181 0.167 0.0819 0 0 0.0875 0.104 0.125 0.167 0.167
## 11 13 0.103 0.154 0.0888 0.0769 0 0.00308 0.0154 0.0308 0.0769 0.154
## # ... with 9 more variables: 75% <dbl>, 90% <dbl>, 95% <dbl>, 99% <dbl>,
## # min <dbl>, max <dbl>, first <dbl>, last <dbl>, n.distinct <int>