##### WORKSPACE & PACKAGES #####
rm(list = ls())

setwd("...")

lapply(c("tidyverse", "lubridate", "fixest", "data.table"), library, character.only = TRUE)

summarize <- dplyr::summarize
select <- dplyr::select
rename <- dplyr::rename

.

##### FUNCTIONS #####
winsorize <- function(object, lower_percentile, upper_percentile, exclude = NULL){
  
  if (is.data.frame(object) == TRUE){
    
    reference_columns <- names(object)
    
    numeric_columns <- sapply(object, is.numeric)
    exclude_columns <- if(!is.null(exclude)) exclude %>% as.vector() %>% as.character() else vector() %>% as.character()
    
    numeric_subset <- object[ , numeric_columns]
    nonnumeric_subset <- cbind(object[ , !numeric_columns], object[ , exclude_columns])
    
    winsorized_variables <- numeric_subset %>% select(- one_of(exclude_columns))
    
    columns <- ncol(winsorized_variables)
    rows <- nrow(winsorized_variables)
    
    for (i in 1:columns){
      L <- (lower_percentile / 100)
      U <- (upper_percentile / 100)
      percentile <- as.matrix(quantile(winsorized_variables[[i]], c(L, U) , na.rm = T))
      min <- percentile[1,1]
      max <- percentile[2,1]
      
      for (j in 1:rows){
        winsorized_variables[j, i] <- if (is.na(winsorized_variables[j, i]) | is.infinite(winsorized_variables[j, i])) NA else
          if (winsorized_variables[j, i] < min) min else if (winsorized_variables[j, i] > max) max else winsorized_variables[j, i]
      }
    }
    
    final_output <- cbind(nonnumeric_subset, winsorized_variables)
    final_output <- final_output[, c(reference_columns, setdiff(names(final_output), reference_columns))]
    
    cat("\n", "In total,", columns, "variables were winsorized:", colnames(winsorized_variables))
    
    return(final_output)
  }
  
  else if (is.vector(object) == TRUE){
    L <- (lower_percentile / 100)
    U <- (upper_percentile / 100)
    percentile <- as.matrix(quantile(object, c(L, U) , na.rm = T))
    min <- percentile[1,1]
    max <- percentile[2,1]
    for (i in 1:length(object)){
      object[i] <- if (is.na(object[i]) | is.infinite(object[i])) NA else
        if (object[i] < min) min else if (object[i] > max) max else object[i]
    }
    return(object)
  }
}

##### MERGE IN ALL DATA #####
load("Environments/AGPT - Environment.RData"); gc()

##### TABLE 1: SUMMARY STATISTICS #####
### NUMBER OF OBSERVATIONS ###
Data %>% filter(!is.na(AvB2)) %>% nrow
Data$inventor_id[!is.na(Data$AvB2)] %>% unique %>% length
Data$patent_number[!is.na(Data$AvB2)] %>% unique %>% length

### PATENT-LEVEL SUMMARY STATISTICS ###
Descriptives <- Data %>% 
  filter(!is.na(AvB2)) %>% 
  group_by(patent_number) %>% 
  mutate(patent_competitor = max(patent_competitor, na.rm = T)) %>% 
  ungroup() %>% 
  transmute(patent_number,
            individual_patent = individual_patent / 100, 
            patent_competitor = patent_competitor / 100,
            early_publication_request = early_publication_request / 100,
            early_publication_request_time = ifelse(early_publication_request == 1, early_publication_request_time, NA),
            patent_originality = patent_originality / 100,
            patent_generality = patent_generality / 100,
            n_inventors) %>% 
  distinct() %>% 
  select(-patent_number)

psych::describe(Descriptives, quant = c(.25, .75)) %>%
  as.data.frame %>%
  rownames_to_column() %>% 
  select(rowname, n, mean, sd, Q0.25, median, Q0.75) %>% 
  print(., digits = 3)

### INVENTOR-LEVEL SUMMARY STATISTICS ###
Descriptives <- Data %>% 
  filter(!is.na(AvB2)) %>% 
  left_join(Moves.all,
            by = c("inventor_id", "priority_year")) %>% 
  left_join(Moves.to_control,
            by = c("inventor_id", "priority_year")) %>% 
  left_join(Inventors %>%
              filter(inventor_country == "US") %>% 
              left_join(Combined %>% 
                          select(patent_number, first_priority_date),
                        by = "patent_number") %>% 
              filter(!is.na(first_priority_date)) %>%
              mutate(priority_year = lubridate::year(first_priority_date)) %>%
              mutate(inventor_state = inventor_state) %>% 
              distinct() %>% 
              na.omit %>% 
              group_by(inventor_state, priority_year) %>%
              summarize(inventors_in_state = n()) %>%
              ungroup(),
            by = c("inventor_state", "priority_year")) %>% 
  left_join(Inventors %>%
              filter(inventor_country == "US") %>% 
              left_join(Combined %>% 
                          select(patent_number, first_priority_date),
                        by = "patent_number") %>% 
              filter(!is.na(first_priority_date)) %>%
              mutate(priority_year = lubridate::year(first_priority_date) + 1) %>%
              mutate(inventor_state = inventor_state) %>%
              distinct() %>% 
              na.omit %>% 
              group_by(inventor_state, priority_year) %>%
              summarize(lag_inventors_in_state = n()) %>%
              ungroup(),
            by = c("inventor_state", "priority_year")) %>% 
  left_join(Assignees %>%
              filter(assignee_country == "US") %>% 
              left_join(Combined %>% 
                          select(patent_number, first_priority_date),
                        by = "patent_number") %>% 
              filter(!is.na(first_priority_date)) %>%
              mutate(priority_year = lubridate::year(first_priority_date)) %>%
              mutate(inventor_state = assignee_state) %>% 
              distinct() %>% 
              na.omit %>% 
              group_by(inventor_state, priority_year) %>%
              summarize(assignees_in_state = n()) %>%
              ungroup(),
            by = c("inventor_state", "priority_year")) %>% 
  left_join(Assignees %>%
              filter(assignee_country == "US") %>% 
              left_join(Combined %>% 
                          select(patent_number, first_priority_date),
                        by = "patent_number") %>% 
              filter(!is.na(first_priority_date)) %>%
              mutate(priority_year = lubridate::year(first_priority_date) + 1) %>%
              mutate(inventor_state = assignee_state) %>% 
              distinct() %>% 
              na.omit %>% 
              group_by(inventor_state, priority_year) %>%
              summarize(lag_assignees_in_state = n()) %>%
              ungroup(),
            by = c("inventor_state", "priority_year")) %>%
  mutate(inventors_in_state = inventors_in_state %>% replace_na(., 0) %>% winsorize(., 1, 99),
         lag_inventors_in_state = lag_inventors_in_state %>% replace_na(., 0) %>% winsorize(., 1, 99),
         assignees_in_state = assignees_in_state %>% replace_na(., 0) %>% winsorize(., 1, 99),
         lag_assignees_in_state = lag_assignees_in_state %>% replace_na(., 0) %>% winsorize(., 1, 99),
         agglomeration_state.inventors = (inventors_in_state - lag_inventors_in_state) %>% winsorize(., 1, 99),
         agglomeration_state.assignees = (assignees_in_state - lag_assignees_in_state) %>% winsorize(., 1, 99),
         agglomeration_state.inventors = agglomeration_state.inventors / 1000,
         agglomeration_state.assignees = agglomeration_state.assignees / 1000) %>% 
  transmute(inventor_id,
            AvB2,
            changes_employer,
            changes_to_control,
            agglomeration_state.inventors,
            agglomeration_state.assignees,
            n_patents,
            same_state) %>% 
  distinct() %>% 
  group_by(inventor_id) %>% 
  summarize(AvB2 = mean(AvB2, na.rm = T),
            changes_employer = mean(changes_employer, na.rm = T),
            changes_to_control = mean(changes_to_control, na.rm = T),
            agglomeration_state.inventors = mean(agglomeration_state.inventors, na.rm = T),
            agglomeration_state.assignees = mean(agglomeration_state.assignees, na.rm = T),
            n_patents = mean(n_patents, na.rm = T),
            same_state = mean(same_state, na.rm = T)) %>% 
  ungroup() %>% 
  select(-inventor_id)

psych::describe(Descriptives, quant = c(.25, .75)) %>%
  as.data.frame %>%
  rownames_to_column() %>% 
  select(rowname, n, mean, sd, Q0.25, median, Q0.75) %>% 
  print(., digits = 3)

### FIRM-LEVEL SUMMARY STATISTICS ###
Descriptives <- Data.firm %>% 
  filter(!is.na(AvB2.ex_ante)) %>% filter(!is.infinite(AvB2.ex_ante)) %>% # drop Texas
  transmute(AvB2.ex_ante,
            n_patents,
            n_inventors,
            same_state,
            perc_treat_state_patents = perc_treat_state_patents / 100,
            perc_cross_treat_control_patent = perc_cross_treat_control_patent / 100,
            avg_perc_disclosure.early = avg_perc_disclosure.early / 100,
            hhi_class = hhi_class / 100)

psych::describe(Descriptives, quant = c(.25, .75)) %>%
  as.data.frame %>%
  rownames_to_column() %>% 
  select(rowname, n, mean, sd, Q0.25, median, Q0.75) %>% 
  print(., digits = 3)

.

##### TABLE 2: ALCATEL VS. BROWN AND PATENT ASSIGNMENT #####
### PANEL A: INDIVIDUAL PATENTS ###
MODEL <- feols(individual_patent ~ AvB2 + sw0(asinh(n_patents) + asinh(n_inventors) + same_state) | firm_year + inventor_id + inventor_city + hqstate_year + class_year, Data); summary(MODEL, cluster = c("firm", "inventor_id", "class")) %>% lapply(., print)

### PANEL B: FUTURE COMPETITOR PATENTS ###
MODEL <- feols(patent_competitor ~ AvB2 + sw0(asinh(n_patents) + asinh(n_inventors) + same_state) | firm_year + inventor_id + inventor_city + hqstate_year + class_year, Data); summary(MODEL, cluster = c("firm", "inventor_id", "class")) %>% lapply(., print)

##### TABLE 3: ALCATEL VS. BROWN AND INVENTOR MOBILITY - MOVES ANYWHERE #####
treatment_date <- ymd("2005-01-01")

Data2 <- Data %>% 
  left_join(Moves.all,
            by = c("inventor_id", "priority_year")) %>% 
  mutate(post = ifelse(filing_date >= (treatment_date %m-% years(1)), 1, 0), # set post period one year earlier to avoid overlapping DVs
         treat = ifelse(!inventor_state %in% c("CA", "DE", "IL", "KS", "MN", "NC", "NV", "UT", "WA"), 1, 0),
         treat2 = ifelse(inventor_state == "TX", NA, treat),
         AvB = post * treat,
         AvB2 = post * treat2,
         changes_employer = changes_employer * 100) %>%
  filter(filing_date >= (treatment_date %m-% years(1))) %>% # use only "2004" for static design (i.e., year prior to treatment)
  filter(filing_date <= (treatment_date - 1)) %>% # use only "2004" for static design (i.e., year prior to treatment)
  distinct(inventor_id, .keep_all = T)

MODEL <- feols(changes_employer ~ AvB2 + sw0(asinh(n_patents) + asinh(n_inventors) + same_state) | firm_year + hqstate_year + class_year, Data2); summary(MODEL, cluster = "firm") %>% lapply(., print)

##### TABLE 4: ALCATEL VS. BROWN AND INVENTOR MOBILITY - MOVES TO CONTROL STATES #####
treatment_date <- ymd("2005-01-01")

Data2 <- Data %>% 
  left_join(Moves.to_control,
            by = c("inventor_id", "priority_year")) %>% 
  mutate(post = ifelse(filing_date >= (treatment_date %m-% years(1)), 1, 0), # set post period one year earlier to avoid overlapping DVs
         treat = ifelse(!inventor_state %in% c("CA", "DE", "IL", "KS", "MN", "NC", "NV", "UT", "WA"), 1, 0),
         treat2 = ifelse(inventor_state == "TX", NA, treat),
         AvB = post * treat,
         AvB2 = post * treat2,
         changes_to_control = changes_to_control %>% replace_na(., 0),
         changes_to_control = changes_to_control * 100) %>% 
  filter(filing_date >= (treatment_date %m-% years(1))) %>% # use only "2004" for static design (i.e., year prior to treatment)
  filter(filing_date <= (treatment_date - 1)) %>% # use only "2004" for static design (i.e., year prior to treatment)
  distinct(inventor_id, .keep_all = T)

MODEL <- feols(changes_to_control ~ AvB2 + sw0(asinh(n_patents) + asinh(n_inventors) + same_state) | firm_year + hqstate_year + class_year, Data2); summary(MODEL, cluster = c("firm", "inventor_id")) %>% lapply(., print)

##### TABLE 5: ALCATEL VS. BROWN AND NEW HIRES IN AGGLOMERATIONS #####
MODEL <- feols(firmchoice ~ AvB2 * agglomeration_state.inventors + same_state | firm_year + inventor_city + hqstate_year, Data.agglomeration); summary(MODEL, cluster = c("firm", "inventor_id"))

MODEL <- feols(firmchoice ~ AvB2 * agglomeration_state.assignees + same_state | firm_year + inventor_city + hqstate_year, Data.agglomeration); summary(MODEL, cluster = c("firm", "inventor_id"))

##### TABLE 6: ALCATEL VS. BROWN AND EARLY PUBLICATION REQUESTS #####
### EARLY PUBLICATION REQUESTS - EXTENSIVE MARGINS ###
# NO INVENTOR FIXED EFFECTS #
MODEL <- feols(early_publication_request ~ AvB2 + sw0(asinh(n_patents) + asinh(n_inventors) + same_state) | firm_year + inventor_city + hqstate_year + class_year, Data %>% filter(disclosure_days < 546)); summary(MODEL, cluster = c("firm", "inventor_id", "class")) %>% lapply(., print)

# INVENTOR FIXED EFFECTS #
MODEL <- feols(early_publication_request ~ AvB2 + sw0(asinh(n_patents) + asinh(n_inventors) + same_state) | firm_year + inventor_id + inventor_city + hqstate_year + class_year, Data %>% filter(disclosure_days < 546)); summary(MODEL, cluster = c("firm", "inventor_id", "class")) %>% lapply(., print)

### EARLY PUBLICATION REQUESTS - INTENSIVE MARGINS ###
# NO INVENTOR FIXED EFFECTS #
MODEL <- feols(asinh(early_publication_request_time) ~ AvB2 + sw0(asinh(n_patents) + asinh(n_inventors) + same_state) | firm_year + inventor_city + hqstate_year + class_year, Data %>% filter(early_publication_request == 100)); summary(MODEL, cluster = c("firm", "inventor_id", "class")) %>% lapply(., print)

# INVENTOR FIXED EFFECTS #
MODEL <- feols(asinh(early_publication_request_time) ~ AvB2 + sw0(asinh(n_patents) + asinh(n_inventors) + same_state) | firm_year + inventor_id + inventor_city + hqstate_year + class_year, Data %>% filter(early_publication_request == 100)); summary(MODEL, cluster = c("firm", "inventor_id", "class")) %>% lapply(., print)

##### TABLE 7: ALCATEL VS. BROWN AND NATURE OF INNOVATION #####
### ORIGINALITY - HHI OF BACKWWARD CITES BASED ON ENTITIES ###
MODEL <- feols(patent_originality ~ AvB2 + sw0(asinh(n_patents) + asinh(n_inventors) + same_state) | firm_year + inventor_id + inventor_city + hqstate_year + class_year, Data); summary(MODEL, cluster = c("firm", "inventor_id", "class")) %>% lapply(., print)

### GENERALITY - HHI OF FORWARD CITES BASED ON ENTITIES ###
MODEL <- feols(patent_generality ~ AvB2 + sw0(asinh(n_patents) + asinh(n_inventors) + same_state) | firm_year + inventor_id + inventor_city + hqstate_year + class_year, Data); summary(MODEL, cluster = c("firm", "inventor_id", "class")) %>% lapply(., print)

##### TABLE 8: ALCATEL VS. BROWN AND AGGREGATE FIRM-LEVEL EFFECTS #####
### ACROSS STATES ###
MODEL <- feols(c(perc_treat_state_patents, # treat vs control patents
                 perc_cross_treat_control_patent) # cross control-treat patents
               ~ AvB2.ex_ante + asinh(n_patents) + asinh(n_inventors) + same_state | filing_year, Data.firm); summary(MODEL, cluster = "firm") %>% lapply(., print)

### OVERALL PORTFOLIO ###
MODEL <- feols(c(avg_perc_disclosure.early, # average patent disclosure timeliness
                 hhi_class) # diversity among technology classes
               ~ AvB2.ex_ante + asinh(n_patents) + asinh(n_inventors) + same_state | filing_year, Data.firm); summary(MODEL, cluster = "firm") %>% lapply(., print)
