code_r/cesloglin/fst_ces_plan_linlog.R

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

#' ---
#' title: "CES log lin formulas working with micro estimates"
#' output:
#'   html_document:
#'     df_print: paged
#'     number_sections: true
#'     toc: true
#'     toc_depth: 3
#'   html_notebook:
#'     number_sections: true
#'   word_document:
#'     number_sections: true
#'   pdf_document:
#'     number_sections: true
#'     toc: true
#'     toc_depth: 3
#' urlcolor: blue
#' ---
#' 
#' Back to **[Fan](https://fanwangecon.github.io/)**'s REconTools Homepage **[Table of Content](https://fanwangecon.github.io/REconTools/)**
#' 
#' # Outline
#' 
#' 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
#'   + @param df tibble data table including variables using svr names below each row is potentially an individual who will receive alternative allocations
#'   + @param svr_A_i string name of the A_i variable, dot product of covariates and coefficients
#'   + @param svr_alpha_i string name of the alpha_i variable, individual specific elasticity information
#'   + @param svr_beta_i string name of the beta_i variable, relative preference weight for each child
#'   + @param svr_N_i string name of the vector of existing inputs, based on which to compute aggregate resource
#'   + @param fl_N_hat float total resource avaible for allocation, if not specific, sum by svr_N_i
#'   + @param fl_rho float preference for equality for the planner
#'   + @return a dataframe that expands the df inputs with additional results.
#' 2. The structure assumes some regression has already taken place to generate the i specific variables listed. and
#' 
#' Doing this allows for lagged intereaction that are time specific in an arbitrary way.
#' 
#' ## Set Up
#' 
## ----GlobalOptions, echo = T, results = 'hide', message=F, warning=F---------------------------------------------------------------------------------------
rm(list = ls(all.names = TRUE))
options(knitr.duplicate.label = 'allow')

## ----loadlib, echo = T, results = 'hide', message=F, warning=F---------------------------------------------------------------------------------------------
library(tidyverse)
library(tidymodels)
library(REconTools)
library(knitr)
library(kableExtra)
# file name
st_file_name = 'fst_ces_plan_linlog'
# Generate R File
purl(paste0(st_file_name, ".Rmd"), output=paste0(st_file_name, ".R"), documentation = 2)
# Generate PDF and HTML
# rmarkdown::render("C:/Users/fan/REconTools/support/function/fs_funceval.Rmd", "pdf_document")
# rmarkdown::render("C:/Users/fan/REconTools/support/function/fs_funceval.Rmd", "html_document")

#' 
#' ## Get Data
#' 
## ----Load Packages and Process Data------------------------------------------------------------------------------------------------------------------------
# Load Library

# Select Cebu Only
df_hw_cebu_m24 <- df_hgt_wgt %>% filter(S.country == 'Cebu' & svymthRound == 24 & prot > 0 & hgt > 0) %>% drop_na()

# Generate Discrete Version of momEdu
df_hw_cebu_m24 <- df_hw_cebu_m24 %>%
    mutate(momEduRound = cut(momEdu,
                             breaks=c(-Inf, 10, Inf),
                             labels=c("MEduLow","MEduHigh"))) %>%
    mutate(hgt0med = cut(hgt0,
                             breaks=c(-Inf, 50, Inf),
                             labels=c("h0low","h0high")))

df_hw_cebu_m24$momEduRound = as.factor(df_hw_cebu_m24$momEduRound)
df_hw_cebu_m24$hgt0med = as.factor(df_hw_cebu_m24$hgt0med)

# Attach
attach(df_hw_cebu_m24)

#' 
#' # Regression with Data and Construct Input Arrays
#' 
#' ## Linear Regression
#' 
## ----Linear Regression-------------------------------------------------------------------------------------------------------------------------------------
# Input Matrix
mt_lincv <- model.matrix(~ hgt0 + wgt0)
mt_linht <- model.matrix(~ sex:hgt0med - 1)

# Regress Height At Month 24 on Nutritional Inputs with controls
rs_hgt_prot_lin = lm(hgt ~ prot:mt_linht + mt_lincv - 1)
print(summary(rs_hgt_prot_lin))
rs_hgt_prot_lin_tidy = tidy(rs_hgt_prot_lin)

#' 
#' ## Log-Linear Regression
#' 
## ----Log-Linear Regression---------------------------------------------------------------------------------------------------------------------------------
# Input Matrix
mt_logcv <- model.matrix(~ hgt0 + wgt0)
mt_loght <- model.matrix(~ sex:hgt0med - 1)

# Log and log regression for month 24
rs_hgt_prot_log = lm(log(hgt) ~ log(prot):mt_loght + mt_logcv - 1)
print(summary(rs_hgt_prot_log))
rs_hgt_prot_log_tidy = tidy(rs_hgt_prot_log)

#' 
#' ## Construct Input Arrays $A_i$ and $\alpha_i$
#' 
## ----Post Regression Input Processing----------------------------------------------------------------------------------------------------------------------

# Generate A_i
ar_Ai_lin <- mt_lincv %*% as.matrix(rs_hgt_prot_lin_tidy %>% filter(!str_detect(term, 'prot')) %>% select(estimate))
ar_Ai_log <- mt_logcv %*% as.matrix(rs_hgt_prot_log_tidy %>% filter(!str_detect(term, 'prot')) %>% select(estimate))

# Generate alpha_i
# Generate A_i
ar_alphai_lin <- mt_linht %*% as.matrix(rs_hgt_prot_lin_tidy %>% filter(str_detect(term, 'prot')) %>% select(estimate))
ar_alphai_log <- mt_loght %*% as.matrix(rs_hgt_prot_log_tidy %>% filter(str_detect(term, 'prot')) %>% select(estimate))


#' 
#' # Optimal Allocations
#' 
#' ## Common Parameters for Optimal Allocation
#' 
## ----Set Allocation Parameters-----------------------------------------------------------------------------------------------------------------------------

# Child Count
df_hw_cebu_m24_full <- df_hw_cebu_m24
it_obs = dim(df_hw_cebu_m24)[1]

# Total Resource Count
ar_prot_data = df_hw_cebu_m24$prot
fl_N_agg = sum(ar_prot_data)

# Vector of Planner Preference
ar_rho = c(seq(-200, -100, length.out=5), seq(-100, -25, length.out=5), seq(-25, -5, length.out=5), seq(-5, -1, length.out=5), seq(-1, -0.01, length.out=5), seq(0.01, 0.25, length.out=5), seq(0.25, 0.99, length.out=5))
ar_rho = unique(ar_rho)

#' 
#' ## Optimal Linear Allocation (CRS)
#' 
#' This also works with any CRS CES.
#' 
#' ### Optimal Linear Allocation Hard-Coded
#' 
## ----Optimal Linear Allocation Hard Code All Rho-----------------------------------------------------------------------------------------------------------
# Optimal Linear Equation
# Planner Inputs
mt_hev_lin = matrix(, nrow = length(ar_rho), ncol = 2)
mt_opti_N = matrix(, nrow = it_obs, ncol = length(ar_rho))

for (it_rho_ctr in seq(1,length(ar_rho))) {
  rho = ar_rho[it_rho_ctr]

  ar_term_b = ar_Ai_lin*(ar_alphai_lin*(1/(rho - 1)))
  ar_term_c = ar_Ai_lin*(ar_alphai_lin*(1/(rho - 1)))
  ar_term_d = (ar_alphai_lin*(rho/(rho - 1)))

  # Child Specific Optimal Allocation Array to Store
  ar_opti_lin = matrix(, nrow = it_obs, ncol = 1)
  for (m in seq(1:it_obs)) {
    fl_topright_q = sum((ar_term_b[m] - ar_term_c)/ar_term_d)
    fl_bottom_q = sum((ar_alphai_lin[m]/ar_alphai_lin)^(rho/(rho-1)))
    fl_opti_q = (fl_N_agg - fl_topright_q)/fl_bottom_q
    ar_opti_lin[m] = fl_opti_q
  }

  # Min and Max 
  ar_opti_lin = pmin(fl_N_agg, pmax(0, ar_opti_lin))
  mt_opti_N[,it_rho_ctr] = ar_opti_lin
  df_hw_cebu_m24_full = cbind(df_hw_cebu_m24_full, ar_opti_lin)

  # Utilities
  fl_v_data_lin = sum((ar_Ai_lin + ar_prot_data*ar_alphai_lin - 70)^rho)^(1/rho)
  fl_v_opti_lin = sum((ar_Ai_lin + ar_opti_lin*ar_alphai_lin - 70)^rho)^(1/rho)

  ## HEV
  fl_hev = (fl_v_opti_lin/fl_v_data_lin - 1)
  mt_hev_lin[it_rho_ctr,1] = rho;
  mt_hev_lin[it_rho_ctr,2] = fl_hev;
}

#' 
#' ### Optimal Linear Allocation Pseudo-Function
#' 
## ----Optimal Linear Allocation Equation Line by Line Test one Rho Selected Individuals---------------------------------------------------------------------
# Randomly test 10 individuals/rho combinations to see if the same results produced by the functional (vectorized) code below as the looped code above.

set.seed(123)
for (it_ctr in seq(1, 10)) {
  # Which Individual to Test
  it_indi_ctr_test = sample(it_obs, 1)
  it_rho_ctr_test = sample(length(ar_rho), 1)
  
  # Get Inputs for Function
  fl_A = ar_Ai_lin[it_indi_ctr_test]
  fl_alpha = ar_alphai_lin[it_indi_ctr_test]
  fl_N = df_hw_cebu_m24$prot[it_indi_ctr_test]
  ar_A = ar_Ai_lin
  ar_alpha = ar_alphai_lin
  fl_rho = ar_rho[it_rho_ctr_test]
  # Existing optimal choice
  fl_opti_loop = mt_opti_N[it_indi_ctr_test, it_rho_ctr_test]
  
  # From Paper Proposition 1
  # top of fraction
  fl_p1_s1 = (fl_A * ((fl_alpha) * (1 / (fl_rho - 1))))
  ar_p1_s2 = (ar_A * ((ar_alpha) * (1 / (fl_rho - 1))))
  ar_p1_s3 = ((ar_alpha) * (fl_rho / (fl_rho - 1)))
  fl_p1_s3 = fl_N_agg - sum((fl_p1_s1 - ar_p1_s2) / ar_p1_s3)
  
  # bottom of fraction
  ar_p2 = (fl_alpha / ar_alpha) ^ (fl_rho / (fl_rho - 1))
  
  # overall
  fl_opti_equa = fl_p1_s3 / sum(ar_p2)
  fl_opti_equa = pmin(fl_N_agg, pmax(0, fl_opti_equa))
  
  print(
    paste0(
      'it_ctr:',
      it_ctr,
      ', fl_rho:',
      fl_rho,
      ', i:',
      it_indi_ctr_test,
      ', fl_opti_equa:',
      fl_opti_equa,
      ', fl_opti_loop:',
      fl_opti_loop
    )
  )
}

#' 
#' ### Optimal Linear Allocation Function DPLYR
#' 
## ----Optimal Linear Allocation Functional Equation DPLYR one rho all individuals---------------------------------------------------------------------------

# Define Explicit Optimal Choice Function
ffi_linear_dplyrdo <- function(fl_A, fl_alpha, fl_rho, ar_A, ar_alpha, fl_N_agg){

  # Apply Function From Paper Proposition 1
  fl_p1_s1 = (fl_A * ((fl_alpha) * (1 / (fl_rho - 1))))
  ar_p1_s2 = (ar_A * ((ar_alpha) * (1 / (fl_rho - 1))))
  ar_p1_s3 = ((ar_alpha) * (fl_rho / (fl_rho - 1)))
  fl_p1_s3 = fl_N_agg - sum((fl_p1_s1 - ar_p1_s2) / ar_p1_s3)
  
  # bottom of fraction
  ar_p2 = (fl_alpha / ar_alpha) ^ (fl_rho / (fl_rho - 1))
  
  # overall
  fl_opti_equa = fl_p1_s3 / sum(ar_p2)
  fl_opti_equa = pmin(fl_N_agg, pmax(0, fl_opti_equa))


  return(fl_opti_equa)
}

# Child Heterogeneity as Matrix
mt_nN_by_nQ_A_alpha = cbind(ar_Ai_lin, ar_alphai_lin, ar_prot_data)
ar_st_col_names = c('fl_A', 'fl_alpha', 'ar_prot_data')
tb_nN_by_nQ_A_alpha <- as_tibble(mt_nN_by_nQ_A_alpha) %>% rename_all(~c(ar_st_col_names))

# fl_A, fl_alpha are from columns of tb_nN_by_nQ_A_alpha
fl_rho = ar_rho[5]
tb_nN_by_nQ_A_alpha = tb_nN_by_nQ_A_alpha %>% rowwise() %>% 
                        mutate(dplyr_eval_opti = ffi_linear_dplyrdo(fl_A, fl_alpha, fl_rho,
                                                                    ar_Ai_lin, ar_alphai_lin,
                                                                    fl_N_agg))
# Show
kable(tb_nN_by_nQ_A_alpha[1:10,]) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "responsive"))


#' 
#' ### Optimal Linear Allocation Function DPLYR All Rho
#' 
## ----Optimal Linear Allocation Functional Equation DPLYR all rho all individuals---------------------------------------------------------------------------

# Generate Data, all individuals specific parameters
mt_nN_by_nQ_A_alpha = cbind(ar_Ai_lin, ar_alphai_lin, ar_prot_data)
ar_st_col_names = c('fl_A', 'fl_alpha', 'ar_prot_data')
tb_nN_by_nQ_A_alpha <- as_tibble(mt_nN_by_nQ_A_alpha) %>% rename_all(~c(ar_st_col_names))
# Duplicate to rhos
tb_nN_by_nQ_A_alpha_mesh_rho <- tb_nN_by_nQ_A_alpha %>% expand_grid(fl_rho = ar_rho) %>%
                                  arrange(fl_A, fl_alpha, fl_rho) %>% 
                                  select(fl_rho, !!!syms(ar_st_col_names))

# fl_A, fl_alpha are from columns of tb_nN_by_nQ_A_alpha
tb_nN_by_nQ_A_alpha_mesh_rho = tb_nN_by_nQ_A_alpha_mesh_rho %>% rowwise() %>% 
                              mutate(dplyr_eval_opti = ffi_linear_dplyrdo(fl_A, fl_alpha, fl_rho,
                                                                    ar_Ai_lin, ar_alphai_lin,
                                                                    fl_N_agg)) %>%
                              ungroup()

# Check if Total Allocations sum Up to Same Level for Each RHO
tb_nN_by_nQ_A_alpha_mesh_rho %>%
  group_by(fl_rho) %>%
  summarise(N_opti_all_sum = sum(dplyr_eval_opti))%>%
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "responsive"))

# Show
kable(tb_nN_by_nQ_A_alpha_mesh_rho[1:50,]) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "responsive"))

#' 
#' 
#' ### Graphical Illustration of Optimal Allocation
#' 
## ----------------------------------------------------------------------------------------------------------------------------------------------------------
tb_hev_lin <- as_tibble(mt_hev_lin) %>% mutate(id = row_number())

lineplot_lin <- tb_hev_lin %>%
    ggplot(aes(x=id, y=V2)) +
        geom_line() +
        geom_point() +
        labs(title = 'HEV and Preference',
             x = 'pref',
             y = 'HEV',
             caption = 'Linear') +
        scale_x_continuous(labels = as.character(tb_hev_lin$V1),
                           breaks = tb_hev_lin$V1)
print(lineplot_lin)


#' 
#' ## Optimal LogLinear Allocation
#' 
#' This also works with any CRS CES.
#' 
#' ### Optimal LogLinear Allocation Hard-Coded
#' 
## ----Optimal Loglinear Allocation--------------------------------------------------------------------------------------------------------------------------

mt_hev_log = matrix(, nrow = length(ar_rho), ncol = 2)

for (it_rho_ctr in seq(1,length(ar_rho))) {
  rho = ar_rho[it_rho_ctr]
  fl_N_hat = sum(df_hw_cebu_m24$prot)

  ar_term_b = ar_Ai_lin*(ar_alphai_lin*(1/(rho - 1)))
  ar_term_c = ar_Ai_lin*(ar_alphai_lin*(1/(rho - 1)))
  ar_term_d = (ar_alphai_lin*(rho/(rho - 1)))

  # Child Specific Optimal Allocation Array to Store
  ar_opti_lin = matrix(, nrow = it_obs, ncol = 1)
  for (m in seq(1:it_obs)) {
    fl_topright_q = sum((ar_term_b[m] - ar_term_c)/ar_term_d)
    fl_bottom_q = sum((ar_alphai_lin[m]/ar_alphai_lin)^(rho/(rho-1)))
    fl_opti_q = (fl_N_hat - fl_topright_q)/fl_bottom_q
    ar_opti_lin[m] = fl_opti_q
  }

  # Min and Max
  ar_opti_lin = pmin(fl_N_hat, pmax(0, ar_opti_lin))
  df_hw_cebu_m24_full = cbind(df_hw_cebu_m24_full, ar_opti_lin)

  # Utilities
  fl_v_data_lin = sum((ar_Ai_lin + prot*ar_alphai_lin - 70)^rho)^(1/rho)
  fl_v_opti_lin = sum((ar_Ai_lin + ar_opti_lin*ar_alphai_lin - 70)^rho)^(1/rho)

  ## HEV
  fl_hev = (fl_v_opti_lin/fl_v_data_lin - 1)
  mt_hev_log[it_rho_ctr,1] = rho;
  mt_hev_log[it_rho_ctr,2] = fl_hev;
}

tb_hev_log <- as_tibble(mt_hev_log) %>% mutate(id = row_number())

lineplot_log <- tb_hev_log %>%
    ggplot(aes(x=id, y=V2)) +
        geom_line() +
        geom_point() +
        labs(title = 'HEV and Preference',
             x = 'pref',
             y = 'HEV',
             caption = 'Linear') +
        scale_x_continuous(labels = as.character(tb_hev_log$V1),
                           breaks = tb_hev_log$V1)
print(lineplot_log)