# NOTE: To load data, you must download both the extract's data and the DDI
# and also set the working directory to the folder with these files (or change the path below).

dir.create('intermediate_data/predictors')

if (!require("ipumsr")) stop("Reading IPUMS data into R requires the ipumsr package. It can be installed using the following command: install.packages('ipumsr')")

ddi <- read_ipums_ddi(paste0("raw_data/union_cps/",cpsname,".xml"))
data <- read_ipums_micro(ddi)


# Bring in crosswalk of old MSA to current MSAs
mxw1 = read_csv('raw_data/union_cps/cps_metro_xw_1.csv') %>% filter(!is.na(METAREA)) %>% select(METAREA,met_standard)

cun = left_join(data,mxw1, by = 'METAREA')

# 1990 is first year (not 1989!)
data %>% filter(!is.na(UNION)) %>% count(YEAR)

# Restrict to only 1990+, keep only ones that are not missing met standard in all 5 years 
undat = cun %>% filter(YEAR >= 1990) %>% group_by(met_standard,YEAR) %>% summarise(cnt = length(YEAR))
undat2 = undat %>% group_by(met_standard) %>% summarise(cnt = length(YEAR))
metros = undat2 %>% filter(cnt == 5) %>% select(met_standard) %>% unlist()

xw = filter(cun, met_standard %in%  metros,!is.na(met_standard), YEAR == 2017) %>% group_by(met_standard,METFIPS) %>% summarise(cnt = length(YEAR)) 
write_csv(xw,'intermediate_data/predictors/cps_metro_xw.csv')


# Now get unionization rate 
union = cun %>% filter(UNION %in% c(1,2,3)) 
union = left_join(union, xw, by = 'met_standard',suffix = c('_raw',''))

mun = union %>% group_by(YEAR,METFIPS) %>% summarise(union_pct = wtd.mean(UNION %in% c(2,3),ASECWT),cnt = length(YEAR)) %>% filter(!is.na(METFIPS))
mun[mun$YEAR == 1990,'YEAR'] = 1989

# From merge, Dayton OH code has been updated since 2013
mun[mun$METFIPS == 19380,'METFIPS'] = 19430

unionall = left_join(union, xw, by = 'met_standard',suffix = c('_raw',''))
unionall[unionall$YEAR == 1990,'YEAR'] = 1989



## Now do all year x metros with 50 obs in that year
bigyears = filter(mun, cnt >= 50)

ua2 = left_join(unionall, bigyears %>% select(-union_pct), by = c('METFIPS','YEAR'),suffix = c('','_big'))

um2 = filter(ua2, !is.na(cnt_big))
ump2 = um2 %>% group_by(YEAR,METFIPS) %>% summarise(union_pct = wtd.mean(UNION %in% c(2,3),ASECWT),cnt = length(YEAR)) #%>% filter(!is.na(METFIPS))

write_csv(ump2, 'intermediate_data/predictors/union_50oneyr_cbsa.csv' )

us2 = filter(ua2, is.na(cnt_big))

usp2 = us2 %>% group_by(YEAR,STATEFIP) %>% summarise(union_pct = wtd.mean(UNION %in% c(2,3),ASECWT),cnt = length(YEAR)) #%>% filter(!is.na(METFIPS))
usp2$st = str_pad(usp2$STATEFIP,2,pad = '0')

write_csv(usp2, 'intermediate_data/predictors/union_50oneyr_st.csv' )


## State-level union density from Hirsch-Macpherson series
## Reshape the wide State_Union_Membership_Density_1964-2018.csv into long form
## with columns fips (zero-padded 2-digit), year, unionrate (proportion, not percent)
us_wide = read_csv('raw_data/unionstats/State_Union_Membership_Density_1964-2018.csv')
us_long = us_wide %>%
  select(-`State Name`, -`State ID`) %>%
  pivot_longer(cols = -fips, names_to = 'year', values_to = 'unionrate') %>%
  mutate(fips = str_pad(fips, 2, pad = '0'),
         year = as.integer(year),
         unionrate = unionrate / 100)

write_csv(us_long, 'intermediate_data/predictors/union_statewide.csv')
