# Compute statistics quoted in the paper and write them all to a single CSV.
# One row per quoted statistic, with the page in the manuscript, a short
# description, and the computed value.

# Bring in MSA names
msanames = read.csv('raw_data/msas/msa_names_2020.csv')
msanames$msa = as.character(msanames$CBSA.Code)
msanames = select(msanames, -CBSA.Code, msaname = CBSA.Title)
msanames$shortname = as.character(msanames$shortname)
msanames$msaname = as.character(msanames$msaname)

# Bring in outcomes and control variables
outv = read_csv('intermediate_data/regression_data/msa_reg_outcomes.csv')
controlv = read_csv('intermediate_data/regression_data/msa_reg_predictors.csv')

# Merge
fulldta = full_join(outv %>% select(-region,-st), controlv, by = c('year','racedef','msa'))
stopifnot(nrow(fulldta) == nrow(outv))
stopifnot(!is.na(fulldta$pop_all_full))

# Make south variable
fulldta$south = fulldta$region == 'South'
fulldta$region = factor(fulldta$region, levels = c('Northeast','Midwest','South','West'))

controls = fulldta %>% filter(racedef == 'expansive', incsamp == 'all_faminc_sqrt')

# Container for paper statistics
paper_stats = tibble(page = integer(), name = character(), value = character())

add_stat = function(page, name, value) {
  if (length(value) == 0) value = NA
  paper_stats <<- bind_rows(paper_stats,
                            tibble(page = as.integer(page),
                                   name = name,
                                   value = as.character(value[1])))
}

# Pull a single MSA x year x variable value.
# Some shortnames are ambiguous across states (Cleveland TN vs OH; Columbia MO vs SC;
# Jackson MI/MS/TN), so use CBSA codes for those.
metro_msa = c(
  'Cleveland'  = '17460',  # Cleveland-Elyria, OH
  'Columbia'   = '17900',  # Columbia, SC
  'Jackson'    = '27140'   # Jackson, MS
)
mv = function(metro, yr, var) {
  if (metro %in% names(metro_msa)) {
    v = controls %>% filter(msa == metro_msa[[metro]], year == yr) %>% pull(!!sym(var))
  } else {
    v = controls %>% filter(shortname == metro, year == yr) %>% pull(!!sym(var))
  }
  if (length(v)) v[1] else NA_real_
}


# Pages 3-4: Huntsville, Baton Rouge, Detroit, San Francisco median ratios
add_stat(3, 'Huntsville 1969 median income ratio',     mv('Huntsville',     1969, 'med_rat'))
add_stat(3, 'Baton Rouge 1969 median income ratio',    mv('Baton Rouge',    1969, 'med_rat'))
add_stat(3, 'Huntsville 2017 median income ratio',     mv('Huntsville',     2017, 'med_rat'))
add_stat(3, 'Baton Rouge 2017 median income ratio',    mv('Baton Rouge',    2017, 'med_rat'))
add_stat(4, 'Detroit 1969 median income ratio',        mv('Detroit',        1969, 'med_rat'))
add_stat(4, 'Detroit 2017 median income ratio',        mv('Detroit',        2017, 'med_rat'))
add_stat(4, 'San Francisco 1969 median income ratio',  mv('San Francisco',  1969, 'med_rat'))
add_stat(4, 'San Francisco 2017 median income ratio',  mv('San Francisco',  2017, 'med_rat'))


# Page 16: BA-attainment examples
add_stat(16, 'San Francisco 2017 % white with BA',     mv('San Francisco',   2017, 'pct_w_col_ba'))
add_stat(16, 'Lake Havasu City 2017 % white with BA',  mv('Lake Havasu City',2017, 'pct_w_col_ba'))
add_stat(16, 'Washington 2017 % Black with BA',        mv('Washington',      2017, 'pct_b_col_ba'))
add_stat(16, 'Scranton 2017 % Black with BA',          mv('Scranton',        2017, 'pct_b_col_ba'))


# Page 22: Pretax vs. disposable income national median ratios (2012 and 2017).
# Text claim: "disparities in posttax income are about 7 percentage points smaller
# than those in pretax income". Equivalently, the median income ratio is ~7pp higher.
natdta_long = read_csv('intermediate_data/sumstats/sumstats_nat_raceex.csv', show_col_types = FALSE) %>%
  filter(race %in% c('b','w'),
         (subpop == 'gq_hh'   & incdef == 'faminc_sqrt') |
         (subpop == 'gq_hh' & incdef == 'spmuinc_disp_sqrt'),
         year %in% c(2014, 2019)) %>%
  select(year, subpop, incdef, race, med) %>%
  pivot_wider(names_from = race, values_from = med) %>%
  mutate(med_rat = b / w,
         spec = ifelse(incdef == 'spmuinc_disp_sqrt', 'disposable', 'pretax baseline'),
         yr_paper = ifelse(year == 2014, 2012, 2017))

for (yr in c(2012, 2017)) {
  row_pre  = natdta_long %>% filter(yr_paper == yr, spec == 'pretax baseline')
  row_disp = natdta_long %>% filter(yr_paper == yr, spec == 'disposable')
  add_stat(22, paste0('National pretax median income ratio, ', yr), row_pre$med_rat[1])
  add_stat(22, paste0('National disposable income median ratio, ', yr), row_disp$med_rat[1])
  add_stat(22, paste0('National disposable - pretax median ratio difference, ', yr),
           row_disp$med_rat[1] - row_pre$med_rat[1])
}
add_stat(22, 'Average disposable - pretax median ratio difference (2012 and 2017)',
         mean(c(natdta_long$med_rat[natdta_long$spec == 'disposable'] -
                natdta_long$med_rat[natdta_long$spec == 'pretax baseline'])))


# Page 31: National median Black-white ratio over time
natdta = read_csv('intermediate_data/sumstats/sumstats_nat_raceex.csv') %>%
  pivot_wider(id_cols = c('year','msa','racedef','subpop','incdef'),
              names_from = 'race', values_from = 'med') %>%
  mutate(med_rat = b / w)

nat_subset = natdta %>%
  filter(incdef == 'faminc_sqrt', subpop == 'all',
         year %in% c(1970, 1980, 1990, 2000, 2009, 2014, 2019))

nat_paper_year = c('1970' = 1969, '1980' = 1979, '1990' = 1989, '2000' = 1999,
                   '2009' = 2007, '2014' = 2012, '2019' = 2017)
for (i in seq_len(nrow(nat_subset))) {
  yr_data = as.character(nat_subset$year[i])
  yr_paper = nat_paper_year[yr_data]
  add_stat(31, paste0('National median income ratio, ', yr_paper),
           nat_subset$med_rat[i])
}


# Pages 31-32: 1969 metros - South high-disparity examples
for (m in c('Huntsville','Montgomery','Baton Rouge','Jackson','Memphis')) {
  add_stat(31, paste0(m, ' 1969 median income ratio'), mv(m, 1969, 'med_rat'))
}


# Page 32: 1969 metros outside the South - small-disparity examples
for (m in c('Cleveland','Seattle','San Jose','Rockford','Dayton','Flint')) {
  add_stat(32, paste0(m, ' 1969 median income ratio'), mv(m, 1969, 'med_rat'))
}


# Page 32: Southern 2017 metros with smallest disparities
for (m in c('San Antonio','Tampa','Nashville')) {
  add_stat(32, paste0(m, ' 2017 median income ratio'), mv(m, 2017, 'med_rat'))
}


# Page 32: Southern metros where inequality lessened the most
for (m in c('Memphis','Orlando','Columbia')) {
  add_stat(32, paste0(m, ' 1969 median income ratio'), mv(m, 1969, 'med_rat'))
  add_stat(32, paste0(m, ' 2017 median income ratio'), mv(m, 2017, 'med_rat'))
}


# Page 33: Northern metros where inequality worsened the most
for (m in c('San Jose','Minneapolis','Milwaukee','Buffalo')) {
  add_stat(33, paste0(m, ' 1969 median income ratio'), mv(m, 1969, 'med_rat'))
  add_stat(33, paste0(m, ' 2017 median income ratio'), mv(m, 2017, 'med_rat'))
}


# Page 33: Aggregate regional ratios, South vs. rest of country
regdta = read_csv('intermediate_data/sumstats/sumstats_region_raceex.csv') %>%
  pivot_wider(id_cols = c('year','racedef','subpop','incdef','south'),
              names_from = 'race', values_from = 'med') %>%
  mutate(med_rat = b / w)

reg_subset = regdta %>%
  filter(incdef == 'faminc_sqrt', subpop == 'all', year %in% c(1970, 2019))

for (i in seq_len(nrow(reg_subset))) {
  yr_paper = ifelse(reg_subset$year[i] == 1970, 1969, 2017)
  reg_label = ifelse(reg_subset$south[i] == 'South', 'South', 'Rest of country')
  add_stat(33, paste0(reg_label, ' ', yr_paper, ' median income ratio'),
           reg_subset$med_rat[i])
}


# Page 34: Raleigh and Austin Black/white income growth vs national
mch = read_csv('output/trends/mchange.csv') %>% filter(racedef == 'expansive')

for (m in c('Raleigh','Austin')) {
  r = mch %>% filter(shortname == m)
  if (nrow(r) > 0) {
    add_stat(34, paste0(m, ' Black income growth as % of national'), r$growp_b[1])
    add_stat(34, paste0(m, ' white income growth as % of national'), r$growp_w[1])
  }
}


# Page 34: Shares of Black population in improving and worsening metros.
# Uses changep_b/changep_w (ratio of MSA's position relative to national in 2017 vs
# 1969) rather than growp_b/growp_w. The two should be equivalent, but the +100 buffer
# in grow_b in graphs_3 distorts growp_b for some low-Black-income MSAs; changep_b is
# what the bbetter/wbetter/bothbetter splits in graphs_3 use and what the article quotes.
imp = mch %>% filter(change_rat >= 0.05)            # local ratio improved 5+ pp
wors = mch %>% filter(change_rat <= -0.05)          # local ratio worsened 5+ pp

if (nrow(imp) > 0) {
  pct_imp_b_grew_faster = sum(imp$pop_b_2017[imp$changep_b > 1], na.rm = TRUE) /
                          sum(imp$pop_b_2017, na.rm = TRUE)
  pct_imp_both_grew_faster = sum(imp$pop_b_2017[imp$changep_b > 1 & imp$changep_w > 1], na.rm = TRUE) /
                             sum(imp$pop_b_2017, na.rm = TRUE)
  add_stat(34, 'Share of Black pop in improving metros where Black income grew faster than national',
           pct_imp_b_grew_faster)
  add_stat(34, 'Share of Black pop in improving metros where both Black and white incomes grew faster than national',
           pct_imp_both_grew_faster)
}

if (nrow(wors) > 0) {
  pct_wors_w_also_fell = sum(wors$pop_b_2017[wors$changep_w < 1], na.rm = TRUE) /
                         sum(wors$pop_b_2017, na.rm = TRUE)
  pct_wors_w_grew_faster = sum(wors$pop_b_2017[wors$changep_w > 1], na.rm = TRUE) /
                           sum(wors$pop_b_2017, na.rm = TRUE)
  add_stat(34, 'Share of Black pop in worsening metros where white income also grew slower than national',
           pct_wors_w_also_fell)
  add_stat(34, 'Share of Black pop in worsening metros where white income grew faster than national',
           pct_wors_w_grew_faster)
}


# Page 37: % Black coefficient illustration (Memphis vs Seattle, 1969)
coefs = read_csv('output/regressions/fulloutput/all_faminc_sqrt/coefs_r_all_faminc_sqrt_med_ldiff_expansive_Baseline.csv')

memphis_pct_b_1969 = mv('Memphis', 1969, 'pct_b_pop')
seattle_pct_b_1969 = mv('Seattle', 1969, 'pct_b_pop')
add_stat(37, 'Memphis 1969 % Black population', memphis_pct_b_1969)
add_stat(37, 'Seattle 1969 % Black population',  seattle_pct_b_1969)

coef_pctb_1969 = coefs %>% filter(year == 1969, var == 'pct_b_pop') %>% pull(coef)
add_stat(37, '1969 coefficient on % Black', coef_pctb_1969)

predicted_pctb_effect = exp(coef_pctb_1969 * (memphis_pct_b_1969 - seattle_pct_b_1969)) - 1
add_stat(37, 'Predicted Memphis-vs-Seattle ratio difference from % Black (proportion)',
         predicted_pctb_effect)

huntsville_medrat_1969 = mv('Huntsville', 1969, 'med_rat')
sanjose_medrat_1969    = mv('San Jose',   1969, 'med_rat')
add_stat(37, 'Huntsville 1969 median income ratio (range endpoint)', huntsville_medrat_1969)
add_stat(37, 'San Jose 1969 median income ratio (range endpoint)',   sanjose_medrat_1969)
total_range_1969 = huntsville_medrat_1969 / sanjose_medrat_1969 - 1
add_stat(37, 'Huntsville/San Jose 1969 ratio difference as fraction', total_range_1969)
add_stat(37, '% Black effect as share of full 1969 range',
         predicted_pctb_effect / total_range_1969)

# Correlation between proportion Black and HRR factor, 1969 and 2017
for (yr in c(1969, 2017)) {
  d = controls %>% filter(year == yr, !is.na(pct_b_pop), !is.na(hrr), !is.na(pop_b))
  c = weights::wtd.cor(d$pct_b_pop, d$hrr, weight = d$pop_b)
  add_stat(37, paste0('Black pop-weighted correlation of % Black with HRR, ', yr),
           c[1, 'correlation'])
}


# Page 39: Manufacturing coefficient illustration (Detroit vs Miami, 1969)
detroit_manuf_1969 = mv('Detroit', 1969, 'cemp_pct_manuf')
miami_manuf_1969   = mv('Miami',   1969, 'cemp_pct_manuf')
add_stat(39, 'Detroit 1969 % employment in manufacturing',  detroit_manuf_1969)
add_stat(39, 'Miami 1969 % employment in manufacturing',    miami_manuf_1969)

coef_manuf_1969 = coefs %>% filter(year == 1969, var == 'cemp_pct_manuf') %>% pull(coef)
add_stat(39, '1969 coefficient on % manufacturing', coef_manuf_1969)
add_stat(39, 'Predicted Detroit-vs-Miami ratio difference from manufacturing (proportion)',
         exp((detroit_manuf_1969 - miami_manuf_1969) * coef_manuf_1969) - 1)


# Page 40: Public administration coefficient illustration (DC vs Dallas, 1969)
dc_pub_1969     = mv('Washington', 1969, 'cemp_pct_public')
dallas_pub_1969 = mv('Dallas',     1969, 'cemp_pct_public')
add_stat(40, 'Washington DC 1969 % employment in public administration', dc_pub_1969)
add_stat(40, 'Dallas 1969 % employment in public administration',         dallas_pub_1969)

coef_pub_1969 = coefs %>% filter(year == 1969, var == 'cemp_pct_public') %>% pull(coef)
add_stat(40, '1969 coefficient on % public administration', coef_pub_1969)
add_stat(40, 'Predicted DC-vs-Dallas ratio multiplier from public administration',
         exp(coef_pub_1969 * (dc_pub_1969 - dallas_pub_1969)))


# Page 41: Segregation coefficient illustration (Chicago vs Riverside, 2017)
chicago_seg_2017   = mv('Chicago',   2017, 'dissimilarity_index')
riverside_seg_2017 = mv('Riverside', 2017, 'dissimilarity_index')
add_stat(41, 'Chicago 2017 dissimilarity index',   chicago_seg_2017)
add_stat(41, 'Riverside 2017 dissimilarity index', riverside_seg_2017)

coef_seg_2017 = coefs %>% filter(year == 2017, var == 'dissimilarity_index') %>% pull(coef)
add_stat(41, '2017 coefficient on dissimilarity index', coef_seg_2017)
add_stat(41, 'Predicted Chicago-vs-Riverside ratio difference from segregation (proportion)',
         exp((chicago_seg_2017 - riverside_seg_2017) * coef_seg_2017) - 1)


# Page 43: BA ratio coefficient illustration (Milwaukee vs Nashville, 2017)
controls$barat = controls$pct_w_col_ba / controls$pct_b_col_ba
milwaukee_barat_2017 = controls %>% filter(shortname == 'Milwaukee', year == 2017) %>% pull(barat)
nashville_barat_2017 = controls %>% filter(shortname == 'Nashville', year == 2017) %>% pull(barat)
milwaukee_lograt_2017 = mv('Milwaukee', 2017, 'pct_ba_lograt_bw')
nashville_lograt_2017 = mv('Nashville', 2017, 'pct_ba_lograt_bw')
add_stat(43, 'Milwaukee 2017 white/Black BA ratio',       milwaukee_barat_2017[1])
add_stat(43, 'Nashville 2017 white/Black BA ratio',       nashville_barat_2017[1])
add_stat(43, 'Milwaukee 2017 logged BA ratio',            milwaukee_lograt_2017)
add_stat(43, 'Nashville 2017 logged BA ratio',            nashville_lograt_2017)

coef_baratio_1969 = coefs %>% filter(year == 1969, var == 'pct_ba_lograt_bw') %>% pull(coef)
coef_baratio_2017 = coefs %>% filter(year == 2017, var == 'pct_ba_lograt_bw') %>% pull(coef)
add_stat(43, '1969 coefficient on logged BA ratio', coef_baratio_1969)
add_stat(43, '2017 coefficient on logged BA ratio', coef_baratio_2017)
add_stat(43, 'Ratio of BA coefficients 2017 / 1969',
         coef_baratio_2017 / coef_baratio_1969)
add_stat(43, 'Predicted Milwaukee-vs-Nashville ratio difference from BA, 2017 coefs',
         exp((milwaukee_lograt_2017 - nashville_lograt_2017) * coef_baratio_2017) - 1)
add_stat(43, 'Predicted Milwaukee-vs-Nashville ratio difference from BA, 1969 coefs',
         exp((milwaukee_lograt_2017 - nashville_lograt_2017) * coef_baratio_1969) - 1)


# Page 44: Incarceration coefficient illustration.
# incar_comb_bwlrat is constructed at the state level; MSAs in the same state share
# the same value. Minneapolis carries the Minnesota value and Atlanta carries the Georgia value.
controls$incarrat = controls$incar_comb_b / controls$incar_comb_w
minn_incar_2017 = controls %>% filter(shortname == 'Minneapolis', year == 2017)
atl_incar_2017  = controls %>% filter(shortname == 'Atlanta',     year == 2017)
add_stat(44, 'Minnesota 2017 logged Black/white incarceration ratio',
         minn_incar_2017$incar_comb_bwlrat[1])
add_stat(44, 'Minnesota 2017 Black/white incarceration rate ratio',
         minn_incar_2017$incarrat[1])
add_stat(44, 'Georgia 2017 logged Black/white incarceration ratio',
         atl_incar_2017$incar_comb_bwlrat[1])
add_stat(44, 'Georgia 2017 Black/white incarceration rate ratio',
         atl_incar_2017$incarrat[1])

coef_incar_2017 = coefs %>% filter(year == 2017, var == 'incar_comb_bwlrat') %>% pull(coef)
add_stat(44, '2017 coefficient on logged incarceration ratio', coef_incar_2017)
add_stat(44, 'Predicted Minnesota-vs-Georgia ratio difference from incarceration (proportion)',
         exp((minn_incar_2017$incar_comb_bwlrat[1] - atl_incar_2017$incar_comb_bwlrat[1]) * coef_incar_2017) - 1)


# Page 45: HRR coefficient illustration (Jackson vs Amarillo, 1969 and 2017)
jackson_hrr_1969  = mv('Jackson',  1969, 'hrr')
amarillo_hrr_1969 = mv('Amarillo', 1969, 'hrr')
add_stat(45, 'Jackson MS HRR factor',     jackson_hrr_1969)
add_stat(45, 'Amarillo TX HRR factor',    amarillo_hrr_1969)

coef_hrr_1969 = coefs %>% filter(year == 1969, var == 'hrrv') %>% pull(coef)
coef_hrr_2017 = coefs %>% filter(year == 2017, var == 'hrrv') %>% pull(coef)
add_stat(45, '1969 coefficient on HRR factor (hrrv)', coef_hrr_1969)
add_stat(45, '2017 coefficient on HRR factor (hrrv)', coef_hrr_2017)
add_stat(45, 'Predicted Jackson-vs-Amarillo ratio shortfall from HRR, 1969 (proportion)',
         1 - exp((jackson_hrr_1969 - amarillo_hrr_1969) * coef_hrr_1969))
add_stat(45, 'Predicted Jackson-vs-Amarillo ratio shortfall from HRR, 2017 (proportion)',
         1 - exp((jackson_hrr_1969 - amarillo_hrr_1969) * coef_hrr_2017))


# Page 47: Decomposition - regional means and component shares
segmsas   = filter(controls, year == 1969,  !is.na(dissimilarity_index)) %>% pull(msa) %>% unique()
incarmsas = filter(controls, year == 2017,  !is.na(incar_comb_bwlrat))   %>% pull(msa) %>% unique()
msas_consist = intersect(segmsas, incarmsas)

for (yr in c(1969, 2017)) {
  d = controls %>% filter(msa %in% msas_consist, year == yr) %>%
    group_by(south) %>%
    summarise(medrat = wtd.mean(med_rat,    pop_b),
              lmed   = wtd.mean(med_ldiff, pop_b), .groups = 'drop')
  for (i in seq_len(nrow(d))) {
    reg = ifelse(d$south[i], 'South', 'Rest of country')
    add_stat(47, paste0(reg, ' ', yr, ' consistent-sample mean median income ratio'),
             d$medrat[i])
    add_stat(47, paste0(reg, ' ', yr, ' consistent-sample mean log-difference'),
             d$lmed[i])
  }
}

# Decomposition shares (page 47): from output/simulations/dta_explode_17_std.csv
dta_explode_17 = read_csv('output/simulations/dta_explode_17_std.csv', show_col_types = FALSE)
dta_explode_69 = read_csv('output/simulations/dta_explode_69_std.csv', show_col_types = FALSE)

# dif_pct = share of the change in the regional gap attributed to each component
for (i in seq_len(nrow(dta_explode_17))) {
  add_stat(47,
           paste0('Decomposition (2017 distribution ordering): ', dta_explode_17$name[i],
                  ' share of regional gap change'),
           dta_explode_17$dif_pct[i])
}
# Footnote on p. 47: alternative ordering applying coefficients first
for (i in seq_len(nrow(dta_explode_69))) {
  add_stat(47,
           paste0('Decomposition (1969 distribution ordering, footnote): ', dta_explode_69$name[i],
                  ' share of regional gap change'),
           dta_explode_69$dif_pct[i])
}


# Page 53: R-squared values for the baseline model
r2s = read_csv('output/regressions/all_r2s.csv', show_col_types = FALSE) %>%
  filter(racedef == 'expansive',
         incsamp == 'all_faminc_sqrt',
         mod == 'Baseline',
         outvar == 'med_ldiff')
for (yr in c(1969, 1979, 1989, 1999, 2007, 2012, 2017)) {
  v = r2s %>% filter(year == yr) %>% pull(r2)
  add_stat(53, paste0('Baseline R-squared, ', yr), v)
}


# Write the assembled CSV
dir.create('output/paper_statistics', showWarnings = FALSE, recursive = TRUE)
write_csv(paper_stats, 'output/paper_statistics/paper_statistics.csv')
