名师讲堂｜使用 R 语言测算上市公司人力资本流动、创新产出及流动方向

3283 字 约 12 分钟 名师讲堂 R 语言上市公司创新与知识

今天给大家分享使用 R 语言测算上市公司人力资本流动、创新产出及流动方向的方法。该指标来源于「空气污染、人力资本流动与创新活力 —— 基于个体专利发明的经验证据」，附件中也提供了该文献的 pdf 文件。

在下面的代码中，上市公司专利数据使用的是之前分享的：

1985～2024 年上市公司与专利数据匹配结果（版本3， 含申请、授权信息）：https://rstata.duanshu.com/#/brief/course/04100321f88b411f90429be934934bff

论文中提到使用 AQI 作为空气质量的测度，不过这个指标通常只有 2014 年之后的，平台上也有分享：

2014～2024 年中国各省市区县分年、分月、逐日 AQI 面板数据：https://rstata.duanshu.com/#/brief/course/aa737515f4234c89bcc593a85ad50bc2

如果想进行更长年度的研究，也可以考虑使用 PM2.5 代替，例如这里我使用的就是：

1980 年 1 月～2024 年 12 月各省市区县地表 PM2.5 质量浓度：https://rstata.duanshu.com/#/brief/course/76698bc838ef49c5ab57633aa64b4e06

首先我们读取专利数据进行初步处理：

  1. 仅保留授权专利；
  2. 根据公开公告号和申请号去除重复专利（专利数据里面同时包含了很多专利的申请和授权公告）；
  3. 去除外观设计专利；
  4. 只保留需要的变量。

实际上原论文也保留了外观设计专利。如果想把外观设计专利纳入分析中，可以保留该部分专利，然后外观设计专利分类号采用 LOC 分类，结构为大类（两位数字）-小类（两位数字）。可以使用前两位数字作为下面的 class 变量。另外也有一些论文认为发明专利的创新性更强，仅选择发明专利进行研究，这样也是不错的，更有说服力，也可以节省一些计算量。

library(tidyverse)
# 读取专利数据
read_csv("2010~2014年上市公司与专利数据匹配结果.csv") %>%
  filter(!is.na(授权公告号), 授权公告号 != "") %>%
  select(newipzlid, 股票代码,
         公开公告号, 申请号, IPC, 年份, 专利类型, 市, 发明人) %>%
  filter(str_detect(专利类型, "发明") | str_detect(专利类型, "实用新型")) %>%
  mutate(公开公告号 = str_remove_all(公开公告号, "[A-Z]$")) %>%
  distinct(股票代码, 公开公告号, .keep_all = T) %>%
  distinct(申请号, .keep_all = T) %>%
  select(-公开公告号, -申请号, -专利类型) -> dfall

dfall
#> # A tibble: 302,390 × 6
#>     newipzlid 股票代码 IPC                                     年份 市    发明人
#>         <dbl> <chr>    <chr>                                  <dbl> <chr> <chr>
#>  1 2010000002 002705   A47J37/08                               2010 佛山市…… 郭建刚; …
#>  2 2010000004 000333   A47J37/12                               2010 佛山市…… 杨云; 冯…
#>  3 2010000005 000333   A47J37/12; A21B5/08                     2010 佛山市…… 杨云; 冯…
#>  4 2010000006 000333   A47J37/12                               2010 佛山市…… 杨云; 冯…
#>  5 2010000091 603630   A61K8/97; A61K8/92; A61K8/66; A61K8/6…  2010 汕头市…… 李世忠
#>  6 2010000188 002728   A61K36/9068; A61P11/14; A61K33/02       2010 江门市…… 许丹青
#>  7 2010000262 301372   B01D24/08; B01D24/46                    2010 北京市…… 葛敬
#>  8 2010000281 000912   B01D53/78; B01D53/50; C05G1/00; C05C1…  2010 泸州市…… 宁忠培; …
#>  9 2010000350 002001   B01J27/232; C07C33/03; C07C29/17        2010 绍兴市…… 朱伟东; …
#> 10 2010000353 600938   B01J29/12; B01J37/02; C10G49/08; C10G…  2010 天津市…… 于海斌; …
#> # ℹ 302,380 more rows

每个专利都可能包含多个 IPC 以及多个申请人，下面使用 tidytext::unnest_tokens() 进行拆分：

# 拆分 IPC 和申请人、保留需要的变量
dfall %>%
  tidytext::unnest_tokens("IPC", "IPC", token = stringr::str_split,
                          pattern = "; ", to_lower = F) %>%
  mutate(IPC = str_remove_all(IPC, "[\\s.]"),
         class = str_sub(IPC, 1, 3)) %>%
  select(-IPC) %>%
  tidytext::unnest_tokens("发明人", "发明人", token = stringr::str_split,
                          pattern = "; ", to_lower = F) %>%
  mutate(发明人 = str_remove_all(发明人, "[\\s.]")) %>%
  distinct() -> patent_data

patent_data
#> # A tibble: 1,408,696 × 6
#>     newipzlid 股票代码  年份 市     发明人 class
#>         <dbl> <chr>    <dbl> <chr>  <chr>  <chr>
#>  1 2010000002 002705    2010 佛山市 郭建刚 A47
#>  2 2010000002 002705    2010 佛山市 曾雄   A47
#>  3 2010000002 002705    2010 佛山市 敬启文 A47
#>  4 2010000004 000333    2010 佛山市 杨云   A47
#>  5 2010000004 000333    2010 佛山市 冯敏强 A47
#>  6 2010000004 000333    2010 佛山市 高天兵 A47
#>  7 2010000005 000333    2010 佛山市 杨云   A47
#>  8 2010000005 000333    2010 佛山市 冯敏强 A47
#>  9 2010000005 000333    2010 佛山市 高天兵 A47
#> 10 2010000005 000333    2010 佛山市 杨云   A21
#> # ℹ 1,408,686 more rows
# 保存
patent_data %>%
  write_rds("patent_data.rds")

再计算创新产出前，我们需要先去除发明人重名的。

根据原文逻辑：同一年、同一发明人、不同城市、不同IPC大类 -> 视为重名：

# 1. 剔除发明人重名问题
patent_data %>%
  group_by(年份, 发明人) %>%
  mutate(
    city_count = n_distinct(市),
    class_count = n_distinct(class)
  ) %>%
  ungroup() %>%
  # 重名条件：同一年同一发明人，且城市数 > 1，且大类数 > 1
  filter(!(city_count > 1 & class_count > 1)) %>%
  select(-city_count, -class_count) -> patent_clean

patent_clean
#> # A tibble: 1,038,381 × 6
#>     newipzlid 股票代码  年份 市     发明人 class
#>         <dbl> <chr>    <dbl> <chr>  <chr>  <chr>
#>  1 2010000002 002705    2010 佛山市 郭建刚 A47
#>  2 2010000002 002705    2010 佛山市 曾雄   A47
#>  3 2010000002 002705    2010 佛山市 敬启文 A47
#>  4 2010000004 000333    2010 佛山市 冯敏强 A47
#>  5 2010000004 000333    2010 佛山市 高天兵 A47
#>  6 2010000005 000333    2010 佛山市 冯敏强 A47
#>  7 2010000005 000333    2010 佛山市 高天兵 A47
#>  8 2010000005 000333    2010 佛山市 冯敏强 A21
#>  9 2010000005 000333    2010 佛山市 高天兵 A21
#> 10 2010000006 000333    2010 佛山市 冯敏强 A47
#> # ℹ 1,038,371 more rows
patent_clean %>%
  write_rds("patent_clean.rds")

# 查看剔除情况
cat("原始数据量：", nrow(patent_data), "\n")
cat("剔除重名后数据量：", nrow(patent_clean), "\n")

# 总发明人数
length(unique(patent_data$发明人))
# 去除重名的：
length(unique(patent_clean$发明人))

# 重名的人占：
(length(unique(patent_data$发明人)) - length(unique(patent_clean$发明人))) / length(unique(patent_data$发明人))

然后分组统计每年每位发明人申请的专利总数即为创新产出：

# 2. 计算创新产出（Patent）
patent_clean %>%
  select(发明人, 年份, newipzlid) %>%
  distinct() %>%
  group_by(发明人, 年份) %>%
  summarise(
    Patent = n(),  # 当年申请专利总数
    .groups = 'drop'
  ) -> innovation_output

innovation_output
#> # A tibble: 334,883 × 3
#>    发明人        年份 Patent
#>    <chr>        <dbl>  <int>
#>  1 A·R·纳戈帕尔  2012      1
#>  2 A·仁加萨米    2010      1
#>  3 A·德特默斯    2013      1
#>  4 A·德特默斯    2014      2
#>  5 A·查克拉      2014      2
#>  6 A·米勒        2013      1
#>  7 A·维斯戈尔    2013      1
#>  8 A·维斯戈尔    2014      2
#>  9 A·马赛拉      2011      2
#> 10 BA彼得斯      2014      1
#> # ℹ 334,873 more rows
innovation_output %>%
  haven::write_dta("innovation_output.dta")

如果同一个发明人在不同的年份（如 2005 年和 2006 年）出现在两个不同的上市公司（A 和 B），如果 2006 年后该发明人不再出现 在 A 企业，同时 2006 年前 B 企业也没有该发明人，则认为发明人在 A 和 B 之间发生了流动：

# 3. 识别人力资本流动（Flow）
# 首先构建发明人-年份-企业面板
patent_clean %>%
  filter(!str_detect(发明人, "不公开"), !str_detect(发明人, "不公布")) %>%
  distinct(发明人, 年份, 股票代码, 市) %>%
  arrange(发明人, 年份) -> inventor_panel

inventor_panel
#> # A tibble: 341,909 × 4
#>    发明人        年份 股票代码 市
#>    <chr>        <dbl> <chr>    <chr>
#>  1 A·R·纳戈帕尔  2012 002415   杭州市
#>  2 A·仁加萨米    2010 688553   内江市
#>  3 A·德特默斯    2013 300195   天津市
#>  4 A·德特默斯    2014 300195   天津市
#>  5 A·查克拉      2014 002415   杭州市
#>  6 A·米勒        2013 300195   天津市
#>  7 A·维斯戈尔    2013 300195   天津市
#>  8 A·维斯戈尔    2014 300195   天津市
#>  9 A·马赛拉      2011 000157   长沙市
#> 10 BA彼得斯      2014 688114   深圳市
#> # ℹ 341,899 more rows
# 标记流动：如果发明人在相邻年份出现在不同企业，且满足流动条件
inventor_panel %>%
  group_by(发明人) %>%
  arrange(年份) %>%
  mutate(
    prev_firm = dplyr::lag(股票代码),
    next_firm = dplyr::lead(股票代码),
    prev_city = dplyr::lag(市),
    next_city = dplyr::lead(市),
    Flow = ifelse(股票代码 != prev_firm &
                    股票代码 != next_firm &
                    !is.na(prev_firm) & !is.na(next_firm), 1, 0)
  ) %>%
  ungroup() -> flow_data

# 发生流动的
flow_data %>%
  filter(Flow > 0)
#> # A tibble: 7,081 × 9
#>    发明人  年份 股票代码 市     prev_firm next_firm prev_city next_city  Flow
#>    <chr>  <dbl> <chr>    <chr>  <chr>     <chr>     <chr>     <chr>     <dbl>
#>  1 丁建新  2010 601360   苏州市 601313    601313    苏州市    苏州市        1
#>  2 丁鹏    2010 002594   深圳市 000063    600021    深圳市    上海市        1
#>  3 万修根  2010 000550   南昌市 200550    200550    南昌市    南昌市        1
#>  4 严军    2010 600021   上海市 600019    000063    上海市    深圳市        1
#>  5 严钦山  2010 000625   重庆市 200625    200625    重庆市    重庆市        1
#>  6 乔艳军  2010 000625   重庆市 200625    200625    重庆市    重庆市        1
#>  7 于忠    2010 600848   上海市 900928    000825    上海市    太原市        1
#>  8 于杰    2010 600019   上海市 688073    600066    上海市    郑州市        1
#>  9 于海龙  2010 601600   北京市 600938    600938    北京市    北京市        1
#> 10 井振伟  2010 000625   重庆市 200625    200625    重庆市    重庆市        1
#> # ℹ 7,071 more rows

结合空气质量数据就可以构造流动方向了。Flow_up 表示发明人 流向空气质量较好的城市，此时 Flow_up 取值为 1，否则为 0：

# 此处代码需下载讲义材料查看～

再添加公司信息，方便下面的分公司汇总：

patent_clean %>%
  left_join(innovation_output, by = c("发明人", "年份")) %>%
  left_join(flow_direction %>% select(发明人, 年份, Flow, Flow_up, Flow_down),
            by = c("发明人", "年份"), relationship = "many-to-many") -> final_data

final_data
#> # A tibble: 1,082,837 × 10
#>     newipzlid 股票代码  年份 市     发明人 class Patent  Flow Flow_up Flow_down
#>         <dbl> <chr>    <dbl> <chr>  <chr>  <chr>  <int> <dbl>   <dbl>     <dbl>
#>  1 2010000002 002705    2010 佛山市 郭建刚 A47      323     0       0         0
#>  2 2010000002 002705    2010 佛山市 曾雄   A47        5     0       0         0
#>  3 2010000002 002705    2010 佛山市 敬启文 A47        2     0       0         0
#>  4 2010000004 000333    2010 佛山市 冯敏强 A47       13     0       0         0
#>  5 2010000004 000333    2010 佛山市 高天兵 A47       13     0       0         0
#>  6 2010000005 000333    2010 佛山市 冯敏强 A47       13     0       0         0
#>  7 2010000005 000333    2010 佛山市 高天兵 A47       13     0       0         0
#>  8 2010000005 000333    2010 佛山市 冯敏强 A21       13     0       0         0
#>  9 2010000005 000333    2010 佛山市 高天兵 A21       13     0       0         0
#> 10 2010000006 000333    2010 佛山市 冯敏强 A47       13     0       0         0
#> # ℹ 1,082,827 more rows
final_data %>%
  haven::write_dta("final_data.dta")

分公司、年份统计流入流出人员数量。继续使用前面的 flow_data 或 final_data，这里可以根据需要分三种方法汇总：

  1. 基于流动数据直接统计;
  2. 根据发明人进行精确的流动识别；
  3. 按流动方向进一步细分。

方法一：基于流动数据直接统计

# 统计企业层面的发明人流动情况
# 1. 统计流入企业的发明人数量
flow_data %>%
  filter(Flow == 1) %>%
  group_by(股票代码, 年份) %>%
  summarise(
    inflow_count = n(),  # 流入该企业的发明人数量
    .groups = 'drop'
  ) -> inflow_by_firm

inflow_by_firm
#> # A tibble: 1,849 × 3
#>    股票代码  年份 inflow_count
#>    <chr>    <dbl>        <int>
#>  1 000002    2010            2
#>  2 000012    2010            2
#>  3 000012    2011            2
#>  4 000012    2012            3
#>  5 000012    2013            7
#>  6 000016    2010           11
#>  7 000016    2011           12
#>  8 000016    2012            9
#>  9 000016    2013           10
#> 10 000016    2014            3
#> # ℹ 1,839 more rows
# 2. 统计流出企业的发明人数量
# 需要找出流动前的企业（prev_firm）
flow_data %>%
  filter(Flow == 1) %>%
  group_by(prev_firm, 年份) %>%
  summarise(
    outflow_count = n(),  # 从该企业流出的发明人数量
    .groups = 'drop'
  ) %>%
  rename(股票代码 = prev_firm) -> outflow_by_firm

outflow_by_firm
#> # A tibble: 1,889 × 3
#>    股票代码  年份 outflow_count
#>    <chr>    <dbl>         <int>
#>  1 000012    2010             3
#>  2 000012    2011             2
#>  3 000012    2012             4
#>  4 000012    2013             1
#>  5 000012    2014             4
#>  6 000016    2010            15
#>  7 000016    2011            12
#>  8 000016    2012            13
#>  9 000016    2013            11
#> 10 000016    2014             5
#> # ℹ 1,879 more rows
# 此处代码需下载讲义材料查看～

方法二：使用更精确的流动识别

# 构建更精确的流动记录
flow_data %>%
  filter(Flow == 1) %>%
  mutate(
    # 流动发生年份（假设在观察到新企业的年份）
    flow_year = 年份,
    # 流出企业
    source_firm = prev_firm,
    # 流入企业
    target_firm = 股票代码
  ) %>%
  select(发明人, flow_year, source_firm, target_firm, 市, prev_city) -> precise_flow_records

precise_flow_records
#> # A tibble: 7,081 × 6
#>    发明人 flow_year source_firm target_firm 市     prev_city
#>    <chr>      <dbl> <chr>       <chr>       <chr>  <chr>
#>  1 丁建新      2010 601313      601360      苏州市 苏州市
#>  2 丁鹏        2010 000063      002594      深圳市 深圳市
#>  3 万修根      2010 200550      000550      南昌市 南昌市
#>  4 严军        2010 600019      600021      上海市 上海市
#>  5 严钦山      2010 200625      000625      重庆市 重庆市
#>  6 乔艳军      2010 200625      000625      重庆市 重庆市
#>  7 于忠        2010 900928      600848      上海市 上海市
#>  8 于杰        2010 688073      600019      上海市 上海市
#>  9 于海龙      2010 600938      601600      北京市 北京市
#> 10 井振伟      2010 200625      000625      重庆市 重庆市
#> # ℹ 7,071 more rows
# 此处代码需下载讲义材料查看～

方法三：按流动方向进一步细分

# 按流动方向统计
flow_direction %>%
  filter(Flow == 1) %>%
  group_by(股票代码, 年份) %>%
  summarise(
    inflow_total = n(),
    inflow_up = sum(Flow_up),    # 因空气质量改善流入
    inflow_down = sum(Flow_down), # 因空气质量恶化流入
    .groups = 'drop'
  ) -> directional_inflow

directional_inflow

flow_direction %>%
  filter(Flow == 1) %>%
  group_by(prev_firm, 年份) %>%
  summarise(
    outflow_total = n(),
    outflow_up = sum(Flow_up),    # 因空气质量改善流出
    outflow_down = sum(Flow_down), # 因空气质量恶化流出
    .groups = 'drop'
  ) %>%
  rename(股票代码 = prev_firm) -> directional_outflow

directional_outflow

# 合并方向性流动数据
# 此处代码需下载讲义材料查看～

最后还可以计算一些统计信息：

cat("企业流动统计摘要：\n")
cat("平均每年流入发明人数量：", mean(firm_mobility_final$inflow_count), "\n")
cat("平均每年流出发明人数量：", mean(firm_mobility_final$outflow_count), "\n")
cat("净流动为正的企业比例：",
    mean(firm_mobility_final$net_flow > 0) * 100, "%\n")

# 按年份查看流动趋势
yearly_trend <- firm_mobility_final %>%
  group_by(年份) %>%
  summarise(
    avg_inflow = mean(inflow_count),
    avg_outflow = mean(outflow_count),
    total_inflow = sum(inflow_count),
    total_outflow = sum(outflow_count)
  )

print(yearly_trend)
#> # A tibble: 5 × 5
#>    年份 avg_inflow avg_outflow total_inflow total_outflow
#>   <dbl>      <dbl>       <dbl>        <int>         <int>
#> 1  2010      0.109       0.109          426           426
#> 2  2011      0.392       0.392         1534          1532
#> 3  2012      0.562       0.562         2199          2197
#> 4  2013      0.550       0.549         2151          2149
#> 5  2014      0.197       0.197          771           771

点击这里跳转到 RStata 短书平台获取附件：名师讲堂｜使用 R 语言测算上市公司人力资本流动、创新产出及流动方向

相关文章

评论