#-----------------------------------------------#
# Title: 02_altruism
#-----------------------------------------------#
# Project: Ordered Logit Empirical Applications
# Author: Carlos Gonzalez
# Date: November 2025
#-----------------------------------------------#
# Description:
# - Creates simulated data for the altruism application
# - Includes all functions related to the analysis of
# the altruism application 
#----------------------------------------------#

#--------------------------------------#
# SIMULATED DATA
#--------------------------------------#


# Step 1: Generate and Order Menu data #
#-------------------------------------#
# Generate menu data
data_gen_alt = function(J = 20, 
                        pi_vector = c(1/100, 1/10),
                        seed = 12345){
  
  # J: Number of menus
  # pi_vector: min and max in interval where prices are uniformly sampled from
  # seed for replicability
  
  # Set seed
  set.seed(seed)
  
  # Prices
  menu_data = replicate(2, runif(J, min = min(pi_vector), 
                                 max = max(pi_vector))) |> as_tibble()
  colnames(menu_data) = c("pi_1", "pi_2")
  
  # Order menus according to \pi_j ratios
  menu_data = menu_data |> mutate(pi_j = pi_2 / pi_1) |>
                           arrange(-pi_j) |>
                           mutate(menu = 1:J)
  
  menu_data
}

# LaTeX export function TABLE 1B
export_table_1b = function(df){
  
  # df: An object created via data_gen_alt()
  
  df |> select(menu, pi_j, pi_1, pi_2) |>
        kbl(caption = "Altruism: Menus indexed by \\pi_j (Table 1B Appendix D)",
            format="latex",
            col.names = c("Menu","\\pi_j", "\\pi_1", "\\pi_2"),
            align="c",
            digits = 3) |>
        kable_minimal(full_width = F,  html_font = "Source Sans Pro")
}

# Step 2: Generate Choice Data
#-----------------------------#
# (CES utility - Logistic distribution)

choice_data_alt_j = function(j, df, N,
                             tau = 0.2, sigma = 0.2, alpha = 0.4){
  
  # j index for menu
  # df: an object created via data_gen_alt()
  # N: number of observations per menu (paper considers {20, 80, 160, 480})
  # tau: location parameter in Logit (main specification sets \tau = 0.2)
  # sigma: scale parameter in Logit (main specification sets \sigma = 0.2)
  # alpha: convexity parameter in altruism CES (main specification sets \sigma = 0.4)
  
  # Recover data for menu j \in J
  pi_1 = df$pi_1[j]
  pi_2 = df$pi_2[j]
  
  # Sample N risk_aversion parameters from logistic distribution
  t = rlogis(N, location = tau, scale = sigma)
  
  # Select optimally under CES utility considerations
  # Noting that: If t <= 0 \implies c_2 = 0
  c_1 = ifelse(t <= 0, 1 / pi_1, 1 / (pi_1 + pi_2 * (t * pi_1 / pi_2)^(1 / (1 - alpha))))
  c_2 = ifelse(t <=0, 0, (1 - pi_1 * c_1) / pi_2)
  
  # Bind data
  cbind(menu = rep(j, N),
        t_real = t,
        c_1 = c_1,
        c_2 = c_2)
}

# Generate data for all menus
choice_data_alt_gen = function(df, N = 160,
                               tau = .2, sigma = .2,
                               alpha = .4, seed = 12345){
  
  # df: an object created via data_gen_alt()
  # N: number of observations per menu (paper considers {20, 80, 160, 480})
  # tau: location parameter in Logit (main specification sets \tau = 0.2)
  # sigma: scale parameter in Logit (main specification sets \sigma = 0.2)
  # alpha: convexity parameter in altruism CES (main specification sets \sigma = 0.4)
  # seed for replicability  
  
  # Set seed
  set.seed(seed)
  
  # Number of menus
  J = nrow(df)
  
  # lapply over choice_data_j
  cons_data = do.call(rbind, 
                      lapply(1:J, function(j) 
                        choice_data_alt_j(j, df, N,
                                          tau = tau, sigma = sigma,
                                          alpha = alpha))) |> as_tibble()
  colnames(cons_data) = c("menu", "t_real", "c_1", "c_2")
  
  # Merge in menu data and recover quantile info based on consumption ratio r (or v)
  cons_data |> left_join(df, by = "menu") |>
               mutate(r = c_2 / c_1,
                      v_j = 1 - exp(-c_2 /c_1),
                      n_obs = N) |>
               group_by(menu) |>
               mutate(quantile = signif(row_number(r) / N * 100)) |>
               ungroup()
  
}

# FIGURE 1B
# Plot CDF of v for each menu
plot_cdf_v = function(df_list, menu_selection, n_menus_random = 10, seed = 12345){
  
  # df_list: A list of one or more datasets created via choice_data_alt_gen()
  # menu_selection: one of c("all", "random", "num_vector")
  #   - all: selects all menus in unique(df$menu)
  #   - random: selects n_menus_random at random from unique(df$menu)
  #             if menu_selection != random, n_menus_random is ignored
  #   - num_vector: a vector of menus to be plotted i.e. c(1, 7, 14, 15)
  # seed for replicability (only relevant when menus are sampled at random)
  
  # Bind datasets together
  df = do.call(rbind, df_list)
  
  # Select menus according to user specification
  df = menu_selection_verif(df, menu_selection, n_menus_random)
  
  # Plot
  cdf_plot =
    df |> mutate(menu = as.factor(menu)) |>
    rename(Menu = menu) |>
    ggplot() +
    geom_line(aes(x = v_j, y = quantile, 
                  color = Menu,
                  linetype = Menu,
                  group = interaction(Menu))) +
    labs(group = "Menu", x = TeX(r"($v_j = 1 - e^{-r_j}$)"), 
         y = "Percentile") +
    theme_minimal() +
    #ggtitle("Empirical CDF for each Menu") +
    facet_wrap(~n_obs, nrow = 1) +
    theme(legend.position = "bottom",           
          legend.title = element_text(size = 16),
          legend.text = element_text(size = 16),
          axis.title.x = element_text(size = 16, 
                                      margin = margin(t = 10, r = 0, b = 0, l = 0)),
          axis.title.y = element_text(size = 16),
          strip.text = element_text(size = 16)) +
    guides(color = guide_legend(nrow = 1))
  
  plot(cdf_plot)
  ggsave("output_final/figure1b_altruism.png", cdf_plot, height = 5 , width = 10)
}


#--------------------------------------#
# NON-PARAMETRIC ANALYSIS
#--------------------------------------#

# Step 3: WARP for Altruism
#--------------------------#

# TABLE 4
# Verify WARP property for each pair of menus
altr_warp = function(menu_df, big_df){
  
  # menu_df: An object created via data_gen_alt()
  # big_df: An object created via choice_data_alt_gen()
  
  # Number of menus
  J = nrow(menu_df)
  
  # Find crossing points of budget line to identify x_{j j'}
  menu_df = menu_df |> mutate(y_point = 1 / pi_2,
                              x_point = 1 / pi_1)
  
  # Initialize store matrices
  inter_matrix = matrix(NA, J, J) # Intersection matrix
  viol_matrix = matrix(NA, J, J)  # Violations matrix
  
  # Find x_{j, j'} (NA when budgets lines do not intersect)
  for (i in 1:J){
    for (j in i:J){
      if((menu_df$y_point[i] > menu_df$y_point[j] & 
          menu_df$x_point[i] > menu_df$x_point[j]) |
         (menu_df$y_point[i] < menu_df$y_point[j] &
          menu_df$x_point[i] < menu_df$x_point[j]) | (i == j)) # Also do not compare a menu to itself
      {
      } else{
        
        # Coordinates of x_{j j'}
        c_1_cross = (menu_df$y_point[i] - menu_df$y_point[j]) / (1 / menu_df$pi_j[i] - 1 / menu_df$pi_j[j])
        c_2_cross = menu_df$y_point[i] - 1 / menu_df$pi_j[i] * c_1_cross
        
        # Find v_j(x_{j j'})
        v_j = 1 - exp(-c_2_cross / c_1_cross)
        
        if(menu_df$pi_j[i] > menu_df$pi_j[j]){ # ONLY a potential violation when F_j < F_j' given pi_j > pi_j'
          inter_matrix[i,j] = v_j              # Hence, we only assign v(x_{j j'}) when pi_j > pi_j'
        } else{inter_matrix[j,i] = v_j} # One must be bigger than the other provided they intersect
      }
    }
  }
  
  # Check if F_j(v(x_{j,j'})) < F_j'(v(x_{j,j'})) [given pi_j < pi_j']
  for (i in 1:J){
    for (j in 1:J){
      if (!is.na(inter_matrix[i,j])){ # Only check comparable menus
        
        # F_j(v(x_{j,j'}))
        F_j = big_df[(big_df$menu == i & big_df$v_j <= inter_matrix[i,j]), ]$quantile |>
          max()
        # F_j'(v(x_{j,j'}))
        F_j_prime = big_df[(big_df$menu == j & big_df$v_j <= inter_matrix[i,j]), ]$quantile |>
          max()
        
        if (F_j < F_j_prime){
          viol_matrix[i,j] = 1
        } else{viol_matrix[i,j] = 0}
      }
    }
  }
  
  # Output
  list(viol_matrix = viol_matrix,
       summary = tibble(menu = 1:J,
                        comp_cases = rowSums(!is.na(viol_matrix)),
                        share_viols = ifelse(comp_cases == 0, 0,
                                             rowSums(viol_matrix, na.rm = TRUE) / comp_cases)),
       total_violations = sum(viol_matrix, na.rm = T))
}

#--------------------------------------#
# SEMI-PARAMETRIC ANALYSIS
#--------------------------------------#

# Step 4. Recover Alpha
#----------------------#

# Recover alpha function
recover_alpha = function(df){
  
  # df: an object created via choice_data_alt_gen()
  
  # Identify smallest quantile such that c_2 > 0 across ALL menus
  min_quantile = df |> group_by(quantile) |>
                       summarize(all_positive = all(c_2 > 0), .groups = "drop") |>
                       filter(all_positive) |>
                       summarize(min_quantile = min(quantile)) |>
                       pull(min_quantile)
  
  # Select quantiles > min_quantile
  df_alpha =  df |> filter(quantile >= min_quantile) |>
                    mutate(r = c_2 / c_1) |>
                    select(menu, pi_j, r, quantile)
  
  # Initialize dataset for computing all pairs of menus
  store = tibble(menu_1 = NA, pi_j1 = NA, r1 = NA, quantile_1 = NA, 
                 menu_2 = NA, pi_j2 = NA, r2 = NA, quantile_2 = NA)
  
  # Construct dataset with all pairs of menus
  for (p in unique(df_alpha$quantile)){
    df_p = df_alpha |> filter(quantile == p) |> as.matrix()
    for (i in 1:(nrow(df_p) - 1)){
      for (j in (i+1):nrow(df_p)){
        store = rbind(store, c(df_p[i,], df_p[j,]))
      }
    }
  }
  
  # Compute alpha for each pair of menus (belonging to the same quantile)
  store = store[-1,] |> mutate(alpha_inputed = log(pi_j1/pi_j2) / log(r1 / r2) + 1,
                               n_obs = df$n_obs[1])
  store
}

# TABLE 5
# Median (induced) alpha
median_alpha = function(df_list, true_alpha){
  
  # df_list: A list of one or more objects created via recover_alpha()
  
  # Bind datasets together and recover median alpha
  df = do.call(rbind, df_list) |> group_by(n_obs) |>
       summarize(alpha = median(alpha_inputed))
  
}

# Step 5. Induced altruism parameters (given alpha)
#-------------------------------------------------#

# Recover induced parameters based on consumption, menu data and recovered alpha
altruism_params = function(df, median_alpha){
  
  # df: an object created via choice_data_alt_gen()
  # median_alpha: median alpha
  
  # Recover induced t parameter and compute quantile
  df |> left_join(median_alpha, by = "n_obs") |>
        mutate(t = pi_j * (c_2 / c_1)^(1-alpha)) |>
        group_by(menu, n_obs) |>
        mutate(quantile_t = signif(row_number(t) / n_obs * 100) ) |>
        ungroup()
}

# FIGURE 3B
# Plot ECDF of induced parameters for a selection of menus
plot_induced_param_alt = function(df_list, menu_selection = "random",
                                  n_menus_random = 10, seed = 12345,
                                  recover_params, median_alpha){
  
  # df_list: a list of one or more objects created via choice_data_alt_gen()
  # menu_selection: one of c("all", "random", "num_vector")
  #   - all: selects all menus in unique(df$menu)
  #   - random: selects n_menus_random at random from unique(df$menu)
  #             if menu_selection != random, n_menus_random is ignored
  #   - num_vector: a vector of menus to be selected i.e. c(1, 7, 14, 15)
  # recover_params: a function which outputs an induced parameter t (like altruism_params)
  # seed: for replicability (only relevant for menu_selection = random)
  # median_alpha: median alpha
  
  
  # Bind datasets together
  df = do.call(rbind, df_list)
  
  # Recover t-parameter via recover_params
  df = recover_params(df, median_alpha = median_alpha)
  
  # Set seed
  set.seed(seed)
  
  # Select menus for plotting according to user specification
  df = menu_selection_verif(df, menu_selection, n_menus_random)
  
  # Plot
  plot_cdf_induced_param =
    df |> filter(t < 1.5) |> # For plotting considerations only
          mutate(Menu = as.factor(menu)) |>
          ggplot() +
          geom_line(aes(x = t, y = quantile_t, 
                    color = Menu,
                    linetype = Menu,
                    group = interaction(Menu))) +
          labs(group = "Menu", x = TeX(r"(Induced Altruism Parameter)"), 
               y = "Percentile") +
          theme_minimal() +
          #ggtitle("Empirical Induced Parameter CDF for each Menu") +
          facet_wrap(~n_obs, nrow = 1) +
          theme(legend.position = "bottom",
                legend.title = element_text(size = 16),
                legend.text = element_text(size = 16),
                strip.text = element_text(size = 16),
                axis.title.y = element_text(size = 16),
                axis.title.x = element_text(size = 16, 
                                      margin = margin(t = 10, r = 0, b = 0, l = 0))) +
          guides(color = guide_legend(nrow = 1))
  
  plot(plot_cdf_induced_param)
  ggsave("output_final/figure3b_altruism.png", plot_cdf_induced_param, 
         height = 5 , width = 10)
  
}

# Step 6. Equality of distributions #
#-----------------------------------#

# Non-parametric testing AD (adjusted for J samples)
# Documentation: https://cran.r-project.org/web/packages/kSamples/kSamples.pdf

similarity_testing_alt = function(df, recover_params,
                                  menu_selection = "all", n_menus_random = 10, seed = 12345,
                                  quant_selection = "percentile",
                                  median_alpha){
  
  # df_list: a list of one or more objects created via choice_data_gen()
  # recover_params: a function which outputs an induced parameter t (like risk_aversion_params)
  # menu_selection: one of c("all", "random", "num_vector")
  #   - all: selects all menus in unique(df$menu)
  #   - random: selects n_menus_random at random from unique(df$menu)
  #             if menu_selection != random, n_menus_random is ignored
  #   - num_vector: a vector of menus to be selected i.e. c(1, 7, 14, 15)
  # seed for replicability (only relevant for menu_selection = random)
  
  # quant_selection: one of c("percentile", "decile", "num_vector")
  #   - percentile: selects quantiles 1:100. quantiles need to be matched exactly
  #   - decile: selects quantiles 10, 20, ..., 100. quantiles need to be matched exactly
  #   - num_vector: a vector of quantiles to be selected i.e. c(19, 31, 42, 89)
  
  # median_alpha: median alpha
  
  # Check and Select relevant menus according to user specification
  df = menu_selection_verif(df, menu_selection, n_menus_random)
  
  # Select quantiles according to user specification
  df = quant_selection_verif(df, quant_selection)
  
  # Recover t-parameter via recover_params
  df = recover_params(df, median_alpha = median_alpha)
  
  # Set seed
  set.seed(seed)
  
  # Select menus according to user specification
  df = menu_selection_verif(df, menu_selection, n_menus_random)
  
  # Testing #
  #---------#
  
  # Drop param = Inf values to compute integral
  if(max(df$t) == Inf){
    print("One or more induced parameters have been inputted to be Infinity")
    df = df |> filter(t < Inf)
  }
  
  test = ad.test(t ~ menu, data = df)
  return(test$ad[1,]) # Version 1 of AD test given continuous data
}

# TABLE 6
similarity_testing_alt_wrap = function(df_list, recover_params,
                                       menu_selection = "all", n_menus_random = 10, seed = 12345,
                                       quant_selection = "percentile",
                                       median_alpha){
  
  # df_list: a list of one or more objects created via choice_data_gen()
  # recover_params: a function which outputs an induced parameter t (like risk_aversion_params)
  # menu_selection: one of c("all", "random", "num_vector")
  #   - all: selects all menus in unique(df$menu)
  #   - random: selects n_menus_random at random from unique(df$menu)
  #             if menu_selection != random, n_menus_random is ignored
  #   - num_vector: a vector of menus to be selected i.e. c(1, 7, 14, 15)
  # seed for replicability (only relevant for menu_selection = random)
  
  # quant_selection: one of c("percentile", "decile", "num_vector")
  #   - percentile: selects quantiles 1:100. quantiles need to be matched exactly
  #   - decile: selects quantiles 10, 20, ..., 100. quantiles need to be matched exactly
  #   - num_vector: a vector of quantiles to be selected i.e. c(19, 31, 42, 89)
  
  # median_alpha: median alpha
  
  test_table =
    sapply(X = df_list, FUN = similarity_testing_alt,
           recover_params = recover_params,
           menu_selection = menu_selection,
           n_menus_random = n_menus_random,
           seed = seed,
           quant_selection = quant_selection,
           median_alpha = median_alpha) |> t() |> as_tibble()
  
  cbind(c("20 obs", "80 obs", "160 obs", "480 obs"), test_table) |>
    kbl(format="latex",
        col.names = c("N obs", "AD Statistic","Stand. AD Stat.", "Asymp p-value"),
        align="c",
        digits = 3) |>
    kable_minimal(full_width = F,  html_font = "Source Sans Pro")
}

#--------------------------------------#
# PARAMETRIC ANALAYSIS
#--------------------------------------#

# Step 7. Mass analysis
#----------------------#
mass_analysis_alt = function(df_list){
  
  # df_list: a list of one or more objects created via choice_data_gen_alt()
  
  # Bind datasets together
  df = do.call(rbind, df_list)
  
  # Mass analysis
  mass_analysis_df =
    df |> mutate(larger_than_0 = ifelse(c_2 == 0, 1, 0)) |>
          group_by(menu, n_obs) |>
          summarize(mass_at_0 = mean(larger_than_0) * 100,
                    largest_v = max(v_j),
                    pi_j = mean(pi_2 / pi_1)) |>
          pivot_wider(names_from = n_obs, 
                      values_from = c("largest_v", "mass_at_0"))
  
  # Output as list  
  list("largest_v_df" = mass_analysis_df |> select(menu, pi_j, starts_with("largest_v")),
       "mass_at_0_df" = mass_analysis_df |> select(menu, pi_j, starts_with("mass_at")))
}

# LaTeX output for mass analysis table
# TABLE 9 and 10
mass_output_alt = function(df){
  
  # df: a tibble created via mass_analysis_alt()
  
  df |> select(menu, pi_j, ends_with("0")) |>
        kbl(format="latex",
            col.names = c("Menu","\\pi_j", "20 obs", "80 obs", "160 obs", "480 obs"),
            align="c",
            digits = 3) |>
        kable_minimal(full_width = F,  html_font = "Source Sans Pro")
}

# Step 8. CLA Property
#---------------------#

# Recover CLA property relevant data (sum of t and sum of log odds)
cla_data_alt = function(df, recover_params, 
                        choices_per_menu = 5, seed = 12345,
                        median_alpha){
  
  # df: an object created via choice_data_gen()
  # choices_per_menu: number of random choices to be selected by menu
  # recover_params: a function which outputs an induced parameter t (like altruism_params)
  # median_alpha: median alpha
  
  # seed for replicability
  
  # Recover t-parameter via recover_params
  df = recover_params(df, median_alpha = median_alpha)
  
  # Initialize df
  cla_df = tibble(menu = NA, t = NA, t_prime = NA, 
                  log_odds_t = NA, log_odds_t_prime = NA,
                  sum_t = NA, sum_log_odds = NA)
  
  # Set seed
  set.seed(seed)
  
  # Recover n observations per menu and merge together
  for (j in 1:length(unique(df$menu))){
    
    # Select menu
    df_menu = df |> filter(menu == j)
    
    # Compute ECDF
    ecdf_t = ecdf(df_menu |> pull(t)) # Including 0s
    
    # Sample interior observations
    index = sample(which(df_menu$t > 0 & df_menu$t < max(df_menu$t)), choices_per_menu)
    
    # Compute log_odds by using ecdf_t
    subset_menu_t = df_menu[index,] |> mutate(ecdf = ecdf_t(t),
                                              log_odds_t = log(ecdf / (1 - ecdf))) |>
      select(t, log_odds_t)
    
    # Create a mirror subset to sum across
    subset_menu_t_prime = subset_menu_t |> rename(t_prime = t,
                                                  log_odds_t_prime = log_odds_t)
    
    # Merge both data sets and sum across them 
    grid_t = expand_grid("t" = subset_menu_t$t, "t_prime" = subset_menu_t$t)
    df_t = grid_t |> left_join(subset_menu_t, by = "t") |>
                     left_join(subset_menu_t_prime, by = "t_prime") |>
                     mutate(sum_t = t + t_prime,
                            sum_log_odds = log_odds_t + log_odds_t_prime,
                            menu = j) |>
                     select(menu, t, t_prime, log_odds_t, 
                            log_odds_t_prime, sum_t, sum_log_odds)
    
    # Bind data
    cla_df = rbind(cla_df, df_t)
    
  }
  
  cla_df[-1,] # Exclude initial empty row
  cla_df$n_obs = unique(df$n_obs)
  cla_df
  
}

# FIGURE 4B
# Plot CLA property
plot_cla_alt = function(df_list, true_params){
  
  # df_list: a list of objects created via cla_data_alt()
  # true_params: a vector of true parameters to be compared to
  
  # Initialize true_params
  tau = true_params[1]
  sigma = true_params[2]
  
  # Bind datasets together
  df = do.call(rbind, df_list)
  
  # Create segment according to true values
  # y = 1 / sigma * x - 2 * tau / sigma
  
  x_start = 0
  x_end = 2.5
  y_start = 1/sigma *(x_start - 2 * tau)
  y_end = 1/sigma *(x_end - 2 * tau)
  
  # Auxiliary vector with true parameter information
  true_log_vector = paste0("Log(", tau,", ", sigma, ")")
  
  # Plot
  cla_plot =
    df |> ggplot() +
          geom_point(aes(x = sum_t, y = sum_log_odds)) +
          geom_segment(aes(x = x_start, y = y_start, xend = x_end, yend = y_end,
                       color = true_log_vector), linewidth = 1) +
          scale_color_manual(values = "red") +
          theme_minimal() +
          labs(x = "Sum of Altruism Parameter t",
               y = TeX(r"(Sum $\log$ $\frac{F_j(t)}{1-F_j(t)}$)"),
               color = "") +
          #ggtitle("CLA Property") +
          facet_wrap(~n_obs, nrow = 1) +
          theme(legend.position = "bottom",           
                legend.title = element_text(size = 16),
                legend.text = element_text(size = 16),
                axis.title.x = element_text(size = 16, 
                                            margin = margin(t = 10, r = 0, b = 0, l = 0)),
                axis.title.y = element_text(size = 16),
                strip.text = element_text(size = 16))
  
  plot(cla_plot)
  ggsave("output_final/figure4b_altruism.png", cla_plot, height = 5 , width = 10)
  
}


# Step 9. Censored logit parametric MLE
#--------------------------------------#

# TABLE 11 A
# Parametric estimation
parametric_estimation_alt = function(df, recover_params, 
                                     median_alpha){
  
  # df: an object created via choice_data_gen_alt()
  # recover_params: a function which outputs an induced parameter t (like altruism_params)
  # media_alpha: median alpha
  
  # Recover t-parameter via recover_params
  if (class(recover_params) == "function"){
    df = recover_params(df, median_alpha = median_alpha)
  }
  
  # Auxiliary variables (see paper for reference and notation)
  M = nrow(df)
  M_0 = nrow(df[df$t == 0,])
  df_geq0 = df |> filter(t > 0) # Strictly positive observations
  
  # Start values in optimization
  start_values = c(mean(df_geq0$t), sd(df_geq0$t))
  
  # Optimization
  optim_results = optim(par = start_values, # Uses plug-in coefficients as starting values
                        fn = likelihood_function_logit,
                        gr = gradient_function_logit,
                        df_geq0 = df_geq0, # Optimize only c_2 > 0 observations
                        M_0 = M_0,
                        M = M,
                        method = 'Nelder-Mead',
                        control=list(fnscale = -1)) # To maximize (not minimize)  
  
  c("tau" = round(optim_results$par[1], 3),
    "sigma" = round(optim_results$par[2], 3),
    "avg_ll" = round(optim_results$value / unique(df$n_obs), 3))
}

# GoF Test via Bootrstrap
#-------------------------#

# Compute AD statistic
ad_statistic_alt = function(df, param_results) {
  
  # df: an object created via choice_data_gen()
  # param_results: a vector created via parametric_estimation_alt()
  
  # Recover parameters given parametric_estimation_alt() output
  tau = param_results["tau"]
  sigma = param_results["sigma"]
  
  # Compute AD test in t > 0 data
  df_geq0 = df |> filter(t > 0)
  
  # True logistic at x = df_geq0$t
  F_logis = plogis((df_geq0$t - tau) / sigma)
  
  # Empirical CDF
  F_emp = ecdf(df_geq0$t)(df_geq0$t)
  
  # AD statistic
  sum((F_emp - F_logis)^2 / (F_logis * (1 - F_logis)))
}

# Simulate data (Logistic) under the null according to parametric bootstrap
simulate_censored_logistic_alt = function(n_draws, param_results) {
  
  # n_draws: number of draws (= nrow(df))
  # param_results: a vector created via parametric_estimation_alt()
  
  tau = param_results["tau"]
  sigma = param_results["sigma"]
  
  t_sim = rlogis(n_draws, location = tau, scale = sigma)
  t_sim[t_sim < 0] = 0  # Apply censoring
  
  # Export data
  tibble(t = t_sim, n_obs = n_draws)
}

# Bootstrap wrap

bootstrap_wrap_alt = function(n_draws, param_results, recover_params, median_alpha){
  
  # n_draws: number of draws (= nrow(df))
  # param_results: a vector created via parametric_estimation_alt()
  # recover_params: a function which outputs an induced parameter t (like altruism_params)
  
  # Simulate data
  df_boot = simulate_censored_logistic_alt(n_draws, param_results)
  
  # Parametric estimation of simulated data
  params_boot = parametric_estimation_alt(df_boot, "no", median_alpha)
  
  # Compute AD statistic
  ad_statistic_alt(df_boot, params_boot)
  
}

# TABLE 11 B
# Bootstrap AD test function
bootstrap_ad_test_alt = function(df, recover_params, B = 1000, seed = 12345, median_alpha) {
  
  # df: an object created via choice_data_gen()
  # recover_params: a function which outputs an induced parameter t (like altruism_params)
  # B: number of bootrstrap draws
  
  # Set seed for replicability
  set.seed(seed)
  
  # Step 1: Estimate parameters via parametric estimation
  params = parametric_estimation_alt(df, recover_params, median_alpha)
  
  # Step 2: Compute AD test on original data
  df = recover_params(df, median_alpha = median_alpha)
  true_ad = ad_statistic_alt(df, params)
  
  # Step 3: Bootstrap
  n_draws = nrow(df)
  
  ad_bootrstrap = future_replicate(B, bootstrap_wrap_alt(
    n_draws = n_draws,
    param_results = params,
    recover_params = recover_params,
    median_alpha = median_alpha))
  
  # Step 4: Compute p-value
  p_val = mean(ad_bootrstrap >= true_ad)
  
  c("true_ad" = true_ad, 
    "p_value" = p_val)
}