In FanWangEcon/PrjOptiAlloc: The Optimal Allocation of Resources Among Heterogeneous Individuals

# try(source("../.Rprofile"))
rm(list = ls(all.names = TRUE))
# Load Libraries
library(stats)
library(tibble)
library(tidyr)
library(stringr)
library(broom)
library(haven)
library(forcats)
library(ggplot2)
library(dplyr)

library(PrjOptiAlloc)

# source('C:/Users/fan/Documents/Dropbox (UH-ECON)/repos/PrjOptiAlloc/R/ffd_opt_datagen.R')

bl_save_img = FALSE

Input and Output

Specify Parameters

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))))

Generate A and alpha

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))
print(paste0('fl_A_1:', fl_A_1, ',fl_A_2:', fl_A_2))
print(paste0('ar_alpha_1l:', ar_alpha_1l))
print(paste0('ar_alpha_2l:', ar_alpha_2l))

Optimal Discrete Choices

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_long

Sort 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)

The Rev Function, compare aginst some Alternative Allocation

Generate All Outcome Combinations

Generate 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)

Generate $\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)
print(df_alpha_2)

Merging 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)

REV Test for One Observed Set of Allocations at One Resource Point

There 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))
df_input_ib <- df_input_ib_all %>% filter(id_alloc==41)
df_input_ib
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),
)

REV for all Possible Allocation Points Compared to What is Optimal

Compute 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)

# show results
head(df_rev_allrho, 50)

REV Summary Statistics

REV 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)
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)

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)
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)