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

Given allocation space parameters (A, alpha) for log linear regression, test log-linear optima allocation line by line. The result is not included in the final version of the paper, because the results are not closed-form, but are based on a root search. The root search is very fast when initial condition is zero.

The objective of this file is to solve the log linear $N_i$ and $H_i$ problem without invoking any functions using bisection with the implicit solution derived in the paper.

File name ffv_opt_solog_bisec_allrw:

• opt: optimal allocation project
• solog: log linear production function
• bisec: implicit formula solved via bisection
• allrw: all code line by line raw original file

Input and Output

There is a dataset with child attributes, nutritional inputs and outputs. Run regression to estimate some input output relationship first. Then generate required inputs for code.

1. Required Input
2. @param df tibble data table including variables using svr names below each row is potentially an individual who will receive alternative allocations
3. @param svr_A_i string name of the A_i variable, dot product of covariates and coefficients, TFP from log linear
4. @param svr_alpha_i string name of the alpha_i variable, individual specific elasticity information, elasticity from log linear
5. @param svr_beta_i string name of the beta_i variable, relative preference weight for each child
6. @param svr_N_i string name of the vector of existing inputs, based on which to compute aggregate resource
7. @param fl_N_hat float total resource avaible for allocation, if not specific, sum by svr_N_i
8. @param fl_rho float preference for equality for the planner
9. @return a dataframe that expands the df inputs with additional results, including bisection solution.
10. The structure assumes some regression has already taken place to generate the i specific variables listed.

Load Packages and Data

Load Dependencies

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

library(tibble)
library(dplyr)
library(tidyr)
library(stringr)
library(broom)
library(ggplot2)

library(PrjOptiAlloc)

library(knitr)
library(kableExtra)

common parameters

some important shared, common fixed parameters specify them here:

# which params to test
it_test = 6

# parameters
fl_rho = 0.20
svr_id_var = 'indi_id'

# it_child_count = n, the number of children
it_n_child_cnt = 2
# it_heter_param = q, number of parameters that are heterogeneous across children
it_q_hetpa_cnt = 2

# choice grid for nutritional feasible choices for each
fl_n_agg = 100
fl_n_min = 0
it_n_choice_cnt_ttest = 3
it_n_choice_cnt_dense = 100

a and alpha input matrix

different types of alpha and a combinations. the algorithm, solutions should work for any alpha and a.

if (it_test == 1) {
  # A. test params 1
  # these are initially what i tested with, positively correlated a and alpha, some have both higher tfp as well as elasticity.
  ar_nn_a = seq(-2, 2, length.out = it_n_child_cnt)
  ar_nn_alpha = seq(0.1, 0.9, length.out = it_n_child_cnt)
  mt_nn_by_nq_a_alpha = cbind(ar_nn_a, ar_nn_alpha)
} else if (it_test == 2) {
  # B. test params 2
  ar_nn_a = c(0,0)
  ar_nn_alpha = c(0.1,0.2)
  mt_nn_by_nq_a_alpha = cbind(ar_nn_a, ar_nn_alpha)
} else if (it_test == 3) {
  # C. test params 3
  ar_nn_a = c(0,0)
  ar_nn_alpha = c(0.1,0.1)
  mt_nn_by_nq_a_alpha = cbind(ar_nn_a, ar_nn_alpha)
} else if (it_test == 4) {
  # D. test params 4
  ar_nn_a = c(0,0)
  ar_nn_alpha = c(0.006686704,0.009251407)
  mt_nn_by_nq_a_alpha = cbind(ar_nn_a, ar_nn_alpha)
} else if (it_test == 5) {
  # E. test params 5
  ar_nn_a = c(4.342501, 4.359987)
  ar_nn_alpha = c(0.006686704,0.009251407)
  mt_nn_by_nq_a_alpha = cbind(ar_nn_a, ar_nn_alpha)
} else if (it_test == 6) {
  # graph b
  # input group two, from the project
  ls_opti_alpha_a <- ffy_opt_dtgch_cbem4()
  df_esti <- ls_opti_alpha_a$df_esti
  df_dtgch_cbem4_bisec <- df_esti[, c('A_log', 'alpha_log')]
  ar_nn_a <- df_dtgch_cbem4_bisec %>% pull(A_log)
  ar_nn_alpha <- df_dtgch_cbem4_bisec %>% pull(alpha_log)
  mt_nn_by_nq_a_alpha = cbind(ar_nn_a, ar_nn_alpha)
}

expand matrix

# choice grid for nutritional feasible choices for each
it_n_choice_cnt_ttest = 3
it_n_choice_cnt_dense = 100
ar_n_choices_ttest = seq(fl_n_min, fl_n_agg, length.out = it_n_choice_cnt_ttest)
ar_n_choices_dense = seq(fl_n_min, fl_n_agg, length.out = it_n_choice_cnt_dense)

# mesh expand
tb_states_choices <- as_tibble(mt_nn_by_nq_a_alpha) %>% rowid_to_column(var=svr_id_var)
tb_states_choices_ttest <- tb_states_choices %>% sample_n(9, replace = FALSE) %>%
                              expand_grid(choices = ar_n_choices_ttest)
tb_states_choices_dense <- tb_states_choices %>% sample_n(9, replace = FALSE) %>%
                              expand_grid(choices = ar_n_choices_dense)

# display
summary(tb_states_choices_dense)
kable(tb_states_choices_ttest) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

Apply Same Function All Rows, some inputs row-specific, other shared

3 points and denser dataframs and define function

# convert matrix to tibble
ar_st_col_names = c(svr_id_var,'fl_a', 'fl_alpha')
tb_states_choices <- tb_states_choices %>% rename_all(~c(ar_st_col_names))
ar_st_col_names = c(svr_id_var,'fl_a', 'fl_alpha', 'fl_n')
tb_states_choices_ttest <- tb_states_choices_ttest %>% rename_all(~c(ar_st_col_names))
tb_states_choices_dense <- tb_states_choices_dense %>% rename_all(~c(ar_st_col_names))

# define implicit function
ffi_nonlin_dplyrdo <- function(fl_a, fl_alpha, fl_n, ar_a, ar_alpha, fl_n_agg, fl_rho){
  # scalar value that are row-specific, in dataframe already: *fl_a*, *fl_alpha*, *fl_n*
  # array and scalars not in dataframe, common all rows: *ar_a*, *ar_alpha*, *fl_n_agg*, *fl_rho*

  # test parameters
  # ar_a = ar_nn_a
  # ar_alpha = ar_nn_alpha
  # fl_n = 100
  # fl_rho = -1
  # fl_n_q = 10

  # apply function
  ar_p1_s1 = exp((fl_a - ar_a)*fl_rho)
  ar_p1_s2 = (fl_alpha/ar_alpha)
  ar_p1_s3 = (1/(ar_alpha*fl_rho - 1))
  ar_p1 = (ar_p1_s1*ar_p1_s2)^ar_p1_s3
  ar_p2 = fl_n^((fl_alpha*fl_rho-1)/(ar_alpha*fl_rho-1))
  ar_overall = ar_p1*ar_p2
  fl_overall = fl_n_agg - sum(ar_overall)

  return(fl_overall)
}

evaluate at three choice points and show table

in the example below, just show results evaluating over three choice points and show table.

# fl_a, fl_alpha are from columns of tb_nn_by_nq_a_alpha
tb_states_choices_ttest_eval = tb_states_choices_ttest %>% rowwise() %>%
                        mutate(dplyr_eval = ffi_nonlin_dplyrdo(fl_a, fl_alpha, fl_n,
                                                               ar_nn_a, ar_nn_alpha,
                                                               fl_n_agg, fl_rho))
# show
kable(tb_states_choices_ttest_eval) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

evaluate at many choice points and show graphically

same as above, but now we evaluate the function over the individuals at many choice points so that we can graph things out.

# fl_a, fl_alpha are from columns of tb_nn_by_nq_a_alpha
tb_states_choices_dense_eval = tb_states_choices_dense %>% rowwise() %>%
                        mutate(dplyr_eval = ffi_nonlin_dplyrdo(fl_a, fl_alpha, fl_n,
                                                               ar_nn_a, ar_nn_alpha,
                                                               fl_n_agg, fl_rho))

# show
dim(tb_states_choices_dense_eval)
summary(tb_states_choices_dense_eval)
lineplot <- tb_states_choices_dense_eval %>%
    ggplot(aes(x=fl_n, y=dplyr_eval)) +
        geom_line() +
        facet_wrap( . ~ indi_id, scales = "free") +
        geom_hline(yintercept=0, linetype="dashed",
                color = "red", size=1)
        labs(title = 'evaluate non-linear functions to search for roots',
             x = 'x values',
             y = 'f(x)',
             caption = 'evaluating the function')
print(lineplot)

dplyr implementation of bisection

initialize matrix

first, initialize the matrix with $a_0$ and $b_0$, the initial min and max points:

# common prefix to make reshaping easier
st_bisec_prefix <- 'bisec_'
svr_a_lst <- paste0(st_bisec_prefix, 'a_0')
svr_b_lst <- paste0(st_bisec_prefix, 'b_0')
svr_fa_lst <- paste0(st_bisec_prefix, 'fa_0')
svr_fb_lst <- paste0(st_bisec_prefix, 'fb_0')

# add initial a and b
tb_states_choices_bisec <- tb_states_choices %>%
                            mutate(!!sym(svr_a_lst) := fl_n_min, !!sym(svr_b_lst) := fl_n_agg)

# evaluate function f(a_0) and f(b_0)
tb_states_choices_bisec <- tb_states_choices_bisec %>% rowwise() %>%
                            mutate(!!sym(svr_fa_lst) := ffi_nonlin_dplyrdo(fl_a, fl_alpha, !!sym(svr_a_lst),
                                                                          ar_nn_a, ar_nn_alpha,
                                                                          fl_n_agg, fl_rho),
                                   !!sym(svr_fb_lst) := ffi_nonlin_dplyrdo(fl_a, fl_alpha, !!sym(svr_b_lst),
                                                                          ar_nn_a, ar_nn_alpha,
                                                                          fl_n_agg, fl_rho))
# summarize
dim(tb_states_choices_bisec)
summary(tb_states_choices_bisec)

iterate and solve for f(p), update f(a) and f(b)

implement the dplyr based concurrent bisection algorithm.

# fl_tol = float tolerance criteria
# it_tol = number of interations to allow at most
fl_tol <- 10^-6
it_tol <- 100
# fl_p_dist2zr = distance to zero to initalize
fl_p_dist2zr <- 1000
it_cur <- 0
while (it_cur <= it_tol && fl_p_dist2zr >= fl_tol ) {

  it_cur <- it_cur + 1

  # new variables
  svr_a_cur <- paste0(st_bisec_prefix, 'a_', it_cur)
  svr_b_cur <- paste0(st_bisec_prefix, 'b_', it_cur)
  svr_fa_cur <- paste0(st_bisec_prefix, 'fa_', it_cur)
  svr_fb_cur <- paste0(st_bisec_prefix, 'fb_', it_cur)

  # evaluate function f(a_0) and f(b_0)
  # 1. generate p
  # 2. generate f_p
  # 3. generate f_p*f_a
  tb_states_choices_bisec <- tb_states_choices_bisec %>% rowwise() %>%
                              mutate(p = ((!!sym(svr_a_lst) + !!sym(svr_b_lst))/2)) %>%
                              mutate(f_p = ffi_nonlin_dplyrdo(fl_a, fl_alpha, p,
                                                              ar_nn_a, ar_nn_alpha,
                                                              fl_n_agg, fl_rho)) %>%
                              mutate(f_p_t_f_a = f_p*!!sym(svr_fa_lst))
  # fl_p_dist2zr = sum(abs(p))
  fl_p_dist2zr <- mean(abs(tb_states_choices_bisec %>% pull(f_p)))

  # update a and b
  tb_states_choices_bisec <- tb_states_choices_bisec %>%
                              mutate(!!sym(svr_a_cur) :=
                                       case_when(f_p_t_f_a < 0 ~ !!sym(svr_a_lst),
                                                 TRUE ~ p)) %>%
                              mutate(!!sym(svr_b_cur) :=
                                       case_when(f_p_t_f_a < 0 ~ p,
                                                 TRUE ~ !!sym(svr_b_lst)))
  # update f(a) and f(b)
  tb_states_choices_bisec <- tb_states_choices_bisec %>%
                              mutate(!!sym(svr_fa_cur) :=
                                       case_when(f_p_t_f_a < 0 ~ !!sym(svr_fa_lst),
                                                 TRUE ~ f_p)) %>%
                              mutate(!!sym(svr_fb_cur) :=
                                       case_when(f_p_t_f_a < 0 ~ f_p,
                                                 TRUE ~ !!sym(svr_fb_lst)))
  # save from last
  svr_a_lst <- svr_a_cur
  svr_b_lst <- svr_b_cur
  svr_fa_lst <- svr_fa_cur
  svr_fb_lst <- svr_fb_cur

  # summar current round
  print(paste0('it_cur:', it_cur, ', fl_p_dist2zr:', fl_p_dist2zr))
  summary(tb_states_choices_bisec %>% select(one_of(svr_a_cur, svr_b_cur, svr_fa_cur, svr_fb_cur)))
}

Reshape wide to long to wide

Table One--very wide

head(tb_states_choices_bisec, 10)
print(paste0('input allocated total sum: ', sum(tb_states_choices_bisec %>% pull(p))))

Table Two--very wide to very long

we want to treat the iteration count information that is the suffix of variable names as a variable by itself. additionally, we want to treat the a,b,fa,fb as a variable. structuring the data very long like this allows for easy graphing and other types of analysis. rather than dealing with many many variables, we have only 3 core variables that store bisection iteration information.

here we use the very nice pivot_longer function. note that to achieve this, we put a common prefix in front of the variables we wanted to convert to long. this is helpful, because we can easily identify which variables need to be reshaped.

# new variables
svr_bisect_iter <- 'biseciter'
svr_abfafb_long_name <- 'varname'
svr_number_col <- 'value'
svr_id_bisect_iter <- paste0(svr_id_var, '_bisect_ier')

# pivot wide to very long
set.seed(123)
tb_states_choices_bisec_long <- 
  tb_states_choices_bisec[sample(dim(tb_states_choices_bisec)[1], 9, replace=FALSE),] %>%
  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
  )

# print
summary(tb_states_choices_bisec_long)
head(tb_states_choices_bisec_long %>% select(-one_of('p','f_p','f_p_t_f_a')), 30)
tail(tb_states_choices_bisec_long %>% select(-one_of('p','f_p','f_p_t_f_a')), 30)

Table Two--very very long to wider again

but the previous results are too long, with the a, b, fa, and fb all in one column as different categories, they are really not different categories, they are in fact different types of variables. so we want to spread those four categories of this variable into four columns, each one representing the a, b, fa, and fb values. the rows would then be uniquly identified by the iteration counter and individual id.

# pivot wide to very long to a little wide
tb_states_choices_bisec_wider <- tb_states_choices_bisec_long %>%
  pivot_wider(
    names_from = !!sym(svr_abfafb_long_name),
    values_from = svr_number_col
  )

# print
summary(tb_states_choices_bisec_wider)
head(tb_states_choices_bisec_wider %>% select(-one_of('p','f_p','f_p_t_f_a')), 30)
tail(tb_states_choices_bisec_wider %>% select(-one_of('p','f_p','f_p_t_f_a')), 30)

Graph Bisection Iteration results

actually we want to graph based on the long results, not the wider. wider easier to view in table.

# graph results
lineplot <- tb_states_choices_bisec_long %>%
    mutate(!!sym(svr_bisect_iter) := as.numeric(!!sym(svr_bisect_iter))) %>%
    filter(!!sym(svr_abfafb_long_name) %in% c('a', 'b')) %>%
    ggplot(aes(x=!!sym(svr_bisect_iter), y=!!sym(svr_number_col),
               colour=!!sym(svr_abfafb_long_name),
               linetype=!!sym(svr_abfafb_long_name),
               shape=!!sym(svr_abfafb_long_name))) +
        facet_wrap( ~ indi_id) +
        geom_line() +
        geom_point() +
        labs(title = 'bisection iteration over individuals until convergence',
             x = 'bisection iteration',
             y = 'a (left side point) and b (right side point) values',
             caption = 'dplyr concurrent bisection nonlinear multple individuals') +
      theme(axis.text.x = element_text(angle = 90, hjust = 1))
print(lineplot)

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