# Bring in raw IPUMS data
# Clean and merge to MSAs

dir.create('intermediate_data/ipums/')

# Bring in geographic crosswalks
geogs = read.csv('output/crosswalks/puma_msa_allyrs_byrace.csv')
geogs$stgeogs = str_pad(geogs$st,2,pad = '0')

# Drop afactor = 0
geogs = filter(geogs, afactor != 0) %>% select(year,race,stpuma,stgeogs,msa,met2,afactor)



# Bring in full raw data from .Rdata format
load(paste0('raw_data/ipums/',ipumsname,'.Rdata'))

data = zap_labels(data)

# For each one, first go through and make clean version at individual level. Save that as .Rdata. 

for(yr in c(1970,1980,1990,2000,2009,2014,2019)){

# Filter data
rreg = filter(data, YEAR == yr)
  
rreg = select(rreg, -INCSUPP)

# Geographies for merging
rreg$st = str_pad(rreg$STATEFIP,width = 2, pad = '0')

if(yr == 1970){
  rreg$stpuma = str_pad(rreg$CNTYGP97, width = 5, pad = '0')
  
}
if(yr == 1980){
  rreg$cg = str_pad(rreg$CNTYGP98, width = 2, pad = '0')
  rreg$stpuma = paste(rreg$st,rreg$cg,sep = '_')
}
if(yr %in% c(1990,2000,2009,2014,2019)){
  rreg$puma = str_pad(rreg$PUMA,width = 5, pad = '0')
  if(yr %in% c(2009,2014)){ 
    rreg = mutate(rreg, 
                  puma = case_when(
                    puma != '77777' ~ puma,
                    puma == '77777' ~ '01801' # Per https://usa.ipums.org/usa//volii/2000pumas.shtml, these PUMAs were combined after Katrina bc of population loss. All are part of New Orleans metro

  ))
  }
  
  
  rreg$stpuma = paste(rreg$st,rreg$puma,sep = '_')
}

stopifnot(!is.na(rreg$stpuma))

# Clean up income variables 
rreg[rreg$INCTOT == 9999999, 'INCTOT'] = NA
rreg[rreg$FTOTINC == 9999999, 'FTOTINC'] = NA
rreg[rreg$INCSS == 99999, 'INCSS'] = NA
rreg[rreg$INCWELFR == 99999, 'INCWELFR'] = NA


# Make HH and family incomes
rreg$all = 1
rreg$adult = rreg$AGE >= 18
rreg = rreg %>% group_by(SERIAL) %>% mutate(hhinc = sum(INCTOT,na.rm = T),
                                            # hhss = sum(incss,na.rm = T),
                                            # hhwelf = sum(incwelfr,na.rm = T),
                                               hhadults = sum(adult),
                                               hhper = sum(all)) %>% ungroup()
rreg = rreg %>% group_by(SERIAL,FAMUNIT) %>% mutate(faminc = sum(INCTOT,na.rm = T),
                                                    # famss = sum(incss,na.rm = T),
                                                    # famwelf = sum(incwelfr,na.rm = T),
                                                    famadults = sum(adult),
                                            famper = sum(all)) %>% ungroup()

# Normalized incomes
rreg = rreg %>% mutate(
  
  hhinc_nadult = case_when(
    hhadults > 0 ~ hhinc / hhadults,
    hhadults == 0 ~ NA_real_ ),
  hhinc_sqrt = hhinc / sqrt(hhper),
  hhinc_norm = hhinc / hhper,
  hhinc_tot = hhinc,
  
  faminc_nadult = case_when(
    famadults > 0 ~ faminc / famadults,
    famadults == 0 ~ NA_real_ ),  
  faminc_sqrt = faminc / sqrt(famper),
  faminc_norm = faminc / famper,
  faminc_tot = faminc

)

# Make race variables
# Two definitions: expansive is anyone who identifies as Black in any way
# Narrow is only people who identify as only NH Black

rreg = rreg %>% mutate(
  race_white = RACE == 1 & HISPAN == 0,
  # race_whitealt = race == 1 & hispan == 0,
  race_black_n = RACE == 2 & HISPAN == 0,
  race_black_e = (RACE == 2 | RACED %in% c(801, 830:845,901:904,917,930:935,950:955,970:973,980:983,985:986,990)),
  race_other_n = race_white == F & race_black_n == F,
  race_other_e = race_white == F & race_black_e == F,
    race_all = T,
  
  race_narrow = case_when(
    race_white == T ~ 'w',
    race_black_n == T ~ 'b',
    race_other_n == T ~ 'o'
    ),
  race_expansive = case_when(
    race_white == T ~ 'w',
    race_black_e == T ~ 'b',
    race_other_e == T ~ 'o'
  )
  )


## Make citizenship variables
rreg = rreg %>% mutate(
  cit_usborn = CITIZEN %in% c(0,1),
  cit_all = T
)

## Make education, employment, group quarters status
rreg = rreg %>% mutate(
  edu_all = T,
  edu_hs = EDUC %in% 0:7,
  edu_somecol = EDUC %in% 7:9,
  edu_ba = EDUC %in% 10:11,
  edu_over25 = AGE >= 25,
  
  emp_all = T,
  emp_emp = EMPSTAT == 1,
  emp_lf = EMPSTAT %in% c(1,2),
  emp_nolf = EMPSTAT == 3,
  
  gq_hh = GQ %in% c(1,2,5,"Households under 1970 definition","Additional households under 1990 definition","Additional households under 2000 definition"),
  gq_inst = GQ %in% c(3,'Group quarters--Institutions'),
  gq_inst_male = GQ == 3 & SEX == 1,
  gq_incarc = GQTYPE == 2,
  gq_incarc_male = GQTYPE == 2 & SEX == 1,
  gq_all = T,
  gq_all_male = SEX == 1,
  
)


# Add in SPM for 2019
if(yr %in% c(2014,2019)){
  
  # Bring in SPM 
  if(yr == 2014) {
    load('intermediate_data/spm_acs/spm_2010_2014.Rdata') 
    spm = spm_2010_2014 %>% select(CBSERIAL,sporder,MULTYEAR,starts_with('spm_'),SPM_Totval,SPM_ID,dispinc,posttax)
    rm(spm_2010_2014)
  }
  if(yr == 2019){
    load('intermediate_data/spm_acs/spm_2015_2019.Rdata') 
    spm = spm_2015_2019 %>% select(CBSERIAL,sporder,MULTYEAR,starts_with('spm_'),SPM_Totval,SPM_ID,dispinc,posttax)
    rm(spm_2015_2019)
  }
  
  # Merge by CBSERIAL, pernum, and multyear
  sdta = left_join(rreg, spm, by = c('CBSERIAL','PERNUM' = 'sporder','MULTYEAR'))
  
  # Check if merge worked
  stopifnot(sdta$SEX == sdta$spm_sex | is.na(sdta$spm_sex))
  stopifnot(sdta$AGE == sdta$spm_age | is.na(sdta$spm_age))
  
  sdta$testwt = sdta$spm_wt * 5
  cor(sdta$testwt, sdta$PERWT,use = 'c') # 0.96 
  
  # Make SPM unit incomes and persons
  sdta = sdta %>% group_by(MULTYEAR, SPM_ID) %>% mutate(
    spmuinc = case_when(
      is.na(SPM_ID) ~ NA,
      !is.na(SPM_ID) ~ sum(INCTOT, na.rm = T)
    ),
    spmuper = case_when(
      is.na(SPM_ID) ~ NA,
      !is.na(SPM_ID) ~ sum(all)
    )
  ) %>% ungroup()
  
  # Adjust for inflation and normalize by square root of SPM unit size
  sdta = sdta %>% mutate(
    cpi_factor = spmuinc / SPM_Totval,
    spmuinc_disp = dispinc * cpi_factor,
    spmuinc_sqrt = spmuinc / sqrt(spmuper),
    spmuinc_posttx = posttax * cpi_factor,
    spmuinc_disp_sqrt = spmuinc_disp / sqrt(spmuper),
    spmuinc_posttx_sqrt = spmuinc_posttx / sqrt(spmuper)
    
  )

  # Note that it's not quite perfect - lots of variation within each year
  # This matches other papers with the same data
  sdta %>% group_by(MULTYEAR) %>% summarise(
    cpi_mean = mean(cpi_factor,na.rm = T),
    cpi_sd = sd(cpi_factor, na.rm = T)
  )
  
  rreg = sdta
  rm(sdta)

}



## Match to CBSAs ##
# Make puma year variable
rreg$pumayear = as.numeric(as.character(rreg$YEAR))

if(yr == 2019) rreg$pumayear = 2010

if(yr == 2014) {
  rreg = rreg %>% mutate(pumayear = case_when(
    MULTYEAR %in% c(2010,2011) ~ 2000,
    MULTYEAR %in% c(2012,2013,2014) ~ 2010
  ))
}

if(yr == 2009) rreg$pumayear = 2000

# Merge on geography data

# In 1970 have to merge by geography + race, in other years don't have to
if(yr == 1970) alldta = left_join(rreg, geogs, by = c('pumayear' = 'year','stpuma','race_expansive' = 'race'))
if(yr != 1970) alldta = left_join(rreg, geogs %>% select(-race), by = c('pumayear' = 'year','stpuma'),relationship = 'many-to-many') # It's many to many because each PUMA has multiple rows in geogs, and also each row in geogs has many individual obs in rreg to merge to

# alldta = left_join(rreg, geogs, by = c('pumayear' = 'year','stpuma'))

# Replace state variable for those missing 
# I think these are cases where a county group crosses state boundaries? In that case a state wouldn't be reported

if(yr == 1970){
  alldta = mutate(alldta,
                  st = case_when(
                    st != '99' ~ st,
                    st == '99' ~ stgeogs
                  ))
}


#No mismatches
stopifnot(!is.na(alldta$msa))
stopifnot(!is.na(alldta$AGE))

#Make new weights
alldta$pcwt = alldta$PERWT * alldta$afactor
alldta$hcwt = alldta$HHWT * alldta$afactor

#Weights almost sum, off by <1000 out of 500,000,000 
# In 2009 off by 3000
sum(alldta$pcwt) - sum(rreg$PERWT)
sum(alldta$hcwt) - sum(rreg$HHWT)

stopifnot(abs(sum(alldta$pcwt) - sum(rreg$PERWT)) < 10000)

# Make internal migrant variable
alldta = mutate(alldta, 
                bplst = str_pad(BPL,3,pad = '0'),
                stbpl = str_pad(st,3,pad = '0'),
                mig_all = T,
                mig_samest = bplst == stbpl,
                mig_difst = bplst != stbpl & BPL < 130, # This includes US territories
                mig_imm = BPL >= 130)



# Select only variables we want to use in the analysis
dtaexp = select(alldta,
                year = YEAR,st,stpuma,msa,met2,afactor,serial = SERIAL,famunit = FAMUNIT,pcwt,hcwt,all,adult,age = AGE,sex = SEX,
                contains('race_'),hhadults,hhper,famadults,famper,contains('spmuper'),
                inctot = INCTOT,incss = INCSS,incwelfr = INCWELFR,contains('hhinc_'),contains('faminc_'),contains('spmuinc'),
                contains('edu_'),contains('emp_'),contains('gq_'),
                bpl = BPL,contains('mig_'),contains('cit_'),contains('mig_'))

# Replace missing spmuincs with 0 - most of the cases where disp is missing and gq_hh is true are dividing by 0
# But really are missing when GQ != T

dtaexp = dtaexp %>% mutate(
  across(contains('spmuinc'), ~ replace_na(.x, 0))
)

save(dtaexp, file = paste('intermediate_data/ipums/ipums_clean_',yr,'.Rdata',sep = ''))

remove(rreg)
remove(alldta)
remove(dtaexp)
}
