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

The objective of this file is to solve the linear $N_i\in{0,1}$ and $H_i$ problem without invoking any functions.

File name ffv_opt_sobin_rkone_allrw:

• opt: optimal allocation project
• sobin: binary provision solution
• rkone: rank at one
• allrw: all code line by line raw original file

Algorithm Outline

Given $N$ individuals, if each individual could or could not receive the provision, given finite resource, and common cost, if total resource avaible could finance input for $M$ individuals, there would be $\frac{N!}{(M-N)!\cdot N!)}$ number of possible choices. Even with 10 individuals,

1. Solve unconstrained relative optimal allocation solutions from the continuous optimal allocation problem.
2. Evaluate relative allocations when allocation for $q$ is equal to $1$
3. This generates a new function where the y-axis values show outcomes for each other individual when individual $q$ allocation is at $1$.
4. Treat the $m$ components of the function as slope and intercept, the $q$ remaining component is the $x$, for each individual, there is a unique $x$
5. The randk order fully determin the sequence of provisions. This is similar to how rank order is determined in the continuous problem, with just a tiny difference of $\alpha$ added in the numerator.

Load Packages and Data

Load Dependencies

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

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

library(broom)
library(ggplot2)
library(REconTools)

library(PrjOptiAlloc)

library(knitr)
library(kableExtra)

Load Data

Generate four categories by initial height and mother's education levels combinations.

# Load Data
tb_mtcars <- as_tibble(rownames_to_column(mtcars, var = "carname")) %>% rowid_to_column()
attach(tb_mtcars)

# Summarize
str(tb_mtcars)
summary(tb_mtcars)

# Generate Discrete Version of continuous variables
tb_mtcars <- tb_mtcars %>%
    mutate(wgtLowHigh = cut(wt,
                             breaks=c(-Inf, 3.1, Inf),
                             labels=c("LowWgt","HighWgt"))) %>%
    mutate(dratLowHigh = cut(drat,
                             breaks=c(-Inf, 3.59, Inf),
                             labels=c("lowRearAxleRatio","highRearAxleRatio")))

# Relabel some factors
tb_mtcars$vs <- factor(tb_mtcars$vs,
                        levels = c(0,1),
                        labels = c("engineVShaped", "engineStraight"))
tb_mtcars$am <- factor(tb_mtcars$am,
                        levels = c(0,1),
                        labels = c("automatic", "manual"))

Exploration

Tabulation

# tabulate
tb_mtcars %>%
  group_by(am, vs) %>%
  summarize(freq = n()) %>%
  pivot_wider(names_from = vs, values_from = freq)

Summarize all, conditional on shift status. am variable is Transmission (0 = automatic, 1 = manual):

# summarize all
REconTools::ff_summ_percentiles(tb_mtcars, bl_statsasrows = FALSE)
# summarize auto transition
REconTools::ff_summ_percentiles(tb_mtcars %>% filter(am == 0), bl_statsasrows = FALSE)
# summarize manual transition
REconTools::ff_summ_percentiles(tb_mtcars %>% filter(am == 1), bl_statsasrows = FALSE)

Exploration Regressions

# A. Regree MPG on horse power and binary if Manual or not (manual = 1)
# 1. more horse power less MPG
# 2. manual larger MPG
print(summary(lm(mpg ~ hp + wt + factor(am) - 1)))

# B. Also incorporate now engine shape, vs = 0 if v-shaped
print(summary(lm(mpg ~ hp + factor(am) - 1)))
print(summary(lm(mpg ~ hp + carb + factor(am):factor(vs) - 1)))
print(summary(lm(mpg ~ hp + carb + factor(vs) + factor(am):factor(vs))))

Regression with Data and Construct Input Arrays

Regress MPG on Manual Transmission

# regress and store results
mt_model <- model.matrix(~ hp + qsec + factor(vs) + factor(am):factor(vs))
rs_mpg_on_auto = lm(mpg ~ mt_model - 1)
print(summary(rs_mpg_on_auto))
rs_mpg_on_auto_tidy = tidy(rs_mpg_on_auto)
rs_mpg_on_auto_tidy

Construct Input Arrays $A_i$ and $\alpha_i$

Multiply coefficient vector by covariate matrix to generate A vector that is child/individual specific.

# Estimates Table
head(rs_mpg_on_auto_tidy, 6)
# Covariates
head(mt_model, 5)

# Covariates coefficients from regression (including constant)
ar_fl_cova_esti <- as.matrix(rs_mpg_on_auto_tidy %>% filter(!str_detect(term, 'am')) %>% select(estimate))
ar_fl_main_esti <- as.matrix(rs_mpg_on_auto_tidy %>% filter(str_detect(term, 'am')) %>% select(estimate))
head(ar_fl_cova_esti, 5)
head(ar_fl_main_esti, 5)

# Select Matrix subcomponents
mt_cova <- as.matrix(as_tibble(mt_model) %>% select(-contains("am")))
mt_intr <- model.matrix(~ factor(vs) - 1)

# Generate A_i, use mt_cova_wth_const
ar_A_m <- mt_cova %*% ar_fl_cova_esti
head(ar_A_m, 5)

# Generate alpha_i
ar_alpha_m <- mt_intr %*% ar_fl_main_esti
head(ar_alpha_m, 5)

Individual Weight

# Child Weight
ar_beta_m <- rep(1/length(ar_A_m), times=length(ar_A_m))

Matrix with Required Inputs for Allocation

# Initate Dataframe that will store all estimates and optimal allocation relevant information
# combine identifying key information along with estimation A, alpha results
# note that we only need indi.id as key
mt_opti <- cbind(ar_alpha_m, ar_A_m, ar_beta_m)
ar_st_varnames <- c('alpha', 'A', 'beta')
df_esti_alpha_A_beta <- as_tibble(mt_opti) %>% rename_all(~c(ar_st_varnames))
tb_key_alpha_A_beta <- bind_cols(tb_mtcars, df_esti_alpha_A_beta) %>%
              select(one_of(c('rowid', 'carname', 'mpg', 'hp', 'qsec', 'vs', 'am', ar_st_varnames)))

# Unique beta, A, and alpha check
tb_opti_unique <- tb_key_alpha_A_beta %>% group_by(!!!syms(ar_st_varnames)) %>%
                    arrange(!!!syms(ar_st_varnames)) %>%
                    summarise(n_obs_group=n())


# Only include currently automatic cars
tb_key_alpha_A_beta <-tb_key_alpha_A_beta %>% filter(am == 'automatic')

# Show cars
head(tb_key_alpha_A_beta, 32)

Optimal Linear Allocations

Parameters for Optimal Allocation

# Child Count
it_obs = dim(tb_opti_unique)[1]

# Vector of Planner Preference
ar_rho <- 1 - (10^(c(seq(-0.5, 0.5, length.out=18))))
ar_rho <- unique(ar_rho)
print(ar_rho)

Optimal binary Allocation (CRS)

This also works with any CRS CES.

# Optimal Linear Equation

# Pull arrays out
ar_A <- tb_key_alpha_A_beta %>% pull(A)
ar_alpha <- tb_key_alpha_A_beta %>% pull(alpha)
ar_beta <- tb_key_alpha_A_beta %>% pull(beta)

# Define Function for each individual m, with hard coded arrays
# This function is a function in the R folder's ffp_opt_sobin.R file as well
ffi_binary_dplyrdo_func <- function(ls_row, fl_rho,
                                    bl_old=FALSE){
  # @param bl_old, weather to use old incorrect solution
  # hard coded inputs are:
  # 1, ar_A
  # 2, ar_alpha
  # 3, ar_beta
  # note follow https://fanwangecon.github.io/R4Econ/support/function/fs_applysapplymutate.html

  svr_A_i = 'A'
  svr_alpha_i = 'alpha'
  svr_beta_i = 'beta'

  fl_alpha <- ls_row[[svr_alpha_i]]
  fl_A <- ls_row[[svr_A_i]]
  fl_beta <- ls_row[[svr_beta_i]]

  ar_left <- (
              ((ar_A + ar_alpha)^fl_rho - (ar_A)^fl_rho)
              /
              ((fl_A + fl_alpha)^fl_rho - (fl_A)^fl_rho)
             )
  ar_right <- ((ar_beta)/(fl_beta))
  ar_rank_val <- ar_left*ar_right

  ar_bl_rank_val <- (ar_rank_val >= 1)

  # Note need the negative multiple in front of array
  ar_rank <- rank((-1)*ar_rank_val, ties.method='min')
  it_rank <- sum(ar_bl_rank_val)

  return(list(it_rank=it_rank,
              ar_rank=ar_rank,
              ar_rank_val=ar_rank_val))
  # return(it_rank)
}

ls_ranks <- ffi_binary_dplyrdo_func(tb_key_alpha_A_beta[2,], 0.1)
it_rank <- ls_ranks$it_rank
ar_rank <- ls_ranks$ar_rank
ar_rank_val <- ls_ranks$ar_rank_val

# accumulate allocation results
tb_opti_alloc_all_rho <- tb_key_alpha_A_beta

# A. First Loop over Planner Preference
# Generate Rank Order
for (it_rho_ctr in seq(1,length(ar_rho))) {
  rho = ar_rho[it_rho_ctr]

  tb_with_rank <- tb_key_alpha_A_beta %>%
    mutate(queue_rank = ffi_binary_dplyrdo_func(tb_key_alpha_A_beta[1,], rho)$ar_rank)

  # m. Keep for df collection individual key + optimal allocation
  # _on stands for optimal nutritional choices
  # _eh stands for expected height
  tb_opti_allocate_wth_key <- tb_with_rank %>% select(one_of('rowid','queue_rank')) %>%
                                rename(!!paste0('rho_c', it_rho_ctr, '_rk') := !!sym('queue_rank'))

  # n. merge optimal allocaiton results from different planner preference
  tb_opti_alloc_all_rho <- tb_opti_alloc_all_rho %>% left_join(tb_opti_allocate_wth_key, by='rowid')
}

# o. print results
print(summary(tb_opti_alloc_all_rho))
str(tb_opti_alloc_all_rho)

# Make Longer
st_bisec_prefix <- 'rho_c'
svr_abfafb_long_name <- 'rho'
svr_bisect_iter <- 'nothing'
svr_number_col <- 'rank'
  tb_opti_alloc_all_rho_long <- tb_opti_alloc_all_rho %>%
  pivot_longer(
    cols = starts_with(st_bisec_prefix),
    names_to = c(svr_abfafb_long_name, svr_bisect_iter),
    names_pattern = paste0(st_bisec_prefix, "(.*)_(.*)"),
    values_to = svr_number_col
  )

# rho as numeric
tb_opti_alloc_all_rho_long <- tb_opti_alloc_all_rho_long %>% mutate(rho_num = as.numeric(rho))
print(summary(tb_opti_alloc_all_rho_long))
str(tb_opti_alloc_all_rho_long)

Bump Plot for Optimal Binary Allocations

st_title = paste0("Binary Allocation Rank, Miles Per Gallon and Manual Shift")
st_subtitle = paste0("Which automatic car to convert to manual shift first?\n",
                     "The expected outcome of interest is miles per gallon\n",
                     "Switch from Auto to Manual shift to improve mpg\n",
                     "Rank = 1, conver this car first from auto to manual shift")
st_caption = paste0("Linear regression of mpg on shift-type. mtcars dataset (R Core Team 2019).")
st_x_lab = paste0("More Efficient --> Cobb-Douglas --> More Equal")
st_y_lab = paste0("Rank Along Binary Allocation Queue ( 1 = highest ranked)")

tb_opti_alloc_all_rho_long %>%
  ggplot(aes(x = rho_num, y = rank, group = carname)) +
    geom_line(aes(color = carname, alpha = 1), size = 2) +
    geom_point(aes(color = carname, alpha = 1), size = 4) +
    scale_x_discrete(expand = c(0.85,0))+
    scale_y_reverse(breaks = 1:nrow(tb_opti_alloc_all_rho_long))+
    theme(legend.position = "none") +
    labs(x = st_x_lab,
         y = st_y_lab,
         title = st_title,
         subtitle = st_subtitle,
         caption = st_caption) +
    ffy_opt_ghthm_dk() +
    geom_text(data = tb_opti_alloc_all_rho, aes(y=rho_c1_rk, x=0.6, label=carname),hjust="right")

FanWangEcon/PrjOptiAlloc documentation built on Jan. 25, 2022, 6:55 a.m.