基于R语言实现markup动态窗口向前年均增长率计算

需求描述 #

现有按gvkey分组的面板数据，需为每个观测计算markup的向前年均增长率：

• 最大计算窗口为5年，首年对应markup的lead(markup)值；
• 当观测后续可用观测不足5个时，按实际数量调整计算窗口；
• 将计算结果作为新列加入原数据框，期望首行结果为-0.00628231878798876，第二行为0.020547945。

现有问题代码 #

原代码缺失窗口自适应逻辑，且存在括号语法错误，导致计算结果不准确：

example %>%
  filter(fyear %in% 1980:2019)%>%
  group_by((gvkey))%>%
  mutate(markupchange = ((((lead(markup)-markup)/markup)+(lead(markup, K =2)-lead(markup)/lead(markup))+(lead(markup, K =3)-lead(markup, k =2)/lead(markup, K=2))+(lead(markup, K =4)-lead(markup, k =3)/lead(markup, K=3))+(lead(markup, K =5)-lead(markup, k =4)/lead(markup, K=4))/5)))

示例数据 #

dput(example)
structure(list(gvkey = c(1001L, 1001L, 1001L, 1003L, 1003L, 1003L, 
1003L, 1003L, 1003L, 1003L, 1004L, 1004L, 1004L, 1004L, 1004L, 
1004L, 1004L, 1004L, 1004L, 1004L, 1004L, 1004L, 1004L, 1004L, 
1004L, 1004L, 1004L, 1004L, 1004L, 1004L), fyear = c(1983L, 1984L, 
1985L, 1983L, 1984L, 1985L, 1986L, 1987L, 1988L, 1989L, 1980L, 
1981L, 1982L, 1983L, 1984L, 1985L, 1986L, 1987L, 1988L, 1989L, 
1990L, 1991L, 1992L, 1993L, 1994L, 1995L, 1996L, 1997L, 1998L, 
1999L), markup = c(3.02456418383518, 2.91714600416106, 2.97620103473762, 
0.628645648836935, 0.538264738598443, 0.74536402337831, 0.89905329776662, 
0.571759161863088, 0.510497237569061, 0.621391904401246, 0.320146680750145, 
0.277978758953348, 0.31442332968701, 0.319433516915814, 0.324865816687745, 
0.335264348013352, 0.328048313395744, 0.326632245360565, 0.340874293859881, 
0.320374201245953, 0.27456562124358, 0.276693369097675, 0.245072145096866, 
0.241026046834387, 0.242841330851661, 0.249635000371186, 0.257903948772679, 
0.262641379065405, 0.261534064206543, 0.22953354130982)), class = c("grouped_df", 
"tbl_df", "tbl", "data.frame"), row.names = c(NA, -30L), groups = structure(list(
    gvkey = c(1001L, 1003L, 1004L), .rows = structure(list(1:3, 
        4:10, 11:30), ptype = integer(0), class = c("vctrs_list_of", 
    "vctrs_vctr", "list"))), class = c("tbl_df", "tbl", "data.frame"
), row.names = c(NA, -3L), .drop = TRUE))

解决方案 #

以下代码实现了自适应窗口的年均增长率计算，核心逻辑为：按gvkey分组后，自动识别每个观测后续的可用观测数（最多5个），计算窗口内每一期的环比增长率，最后取算术平均值作为向前年均增长率。

library(dplyr)
library(purrr)

# 加载示例数据（若已加载可跳过）
example <- structure(list(gvkey = c(1001L, 1001L, 1001L, 1003L, 1003L, 1003L, 
1003L, 1003L, 1003L, 1003L, 1004L, 1004L, 1004L, 1004L, 1004L, 
1004L, 1004L, 1004L, 1004L, 1004L, 1004L, 1004L, 1004L, 1004L, 
1004L, 1004L, 1004L, 1004L, 1004L, 1004L), fyear = c(1983L, 1984L, 
1985L, 1983L, 1984L, 1985L, 1986L, 1987L, 1988L, 1989L, 1980L, 
1981L, 1982L, 1983L, 1984L, 1985L, 1986L, 1987L, 1988L, 1989L, 
1990L, 1991L, 1992L, 1993L, 1994L, 1995L, 1996L, 1997L, 1998L, 
1999L), markup = c(3.02456418383518, 2.91714600416106, 2.97620103473762, 
0.628645648836935, 0.538264738598443, 0.74536402337831, 0.89905329776662, 
0.571759161863088, 0.510497237569061, 0.621391904401246, 0.320146680750145, 
0.277978758953348, 0.31442332968701, 0.319433516915814, 0.324865816687745, 
0.335264348013352, 0.328048313395744, 0.326632245360565, 0.340874293859881, 
0.320374201245953, 0.27456562124358, 0.276693369097675, 0.245072145096866, 
0.241026046834387, 0.242841330851661, 0.249635000371186, 0.257903948772679, 
0.262641379065405, 0.261534064206543, 0.22953354130982)), class = c("grouped_df", 
"tbl_df", "tbl", "data.frame"), row.names = c(NA, -30L), groups = structure(list(
    gvkey = c(1001L, 1003L, 1004L), .rows = structure(list(1:3, 
        4:10, 11:30), ptype = integer(0), class = c("vctrs_list_of", 
    "vctrs_vctr", "list"))), class = c("tbl_df", "tbl", "data.frame"
), row.names = c(NA, -3L), .drop = TRUE))

# 计算自适应窗口的年均增长率
result <- example %>%
  filter(fyear %in% 1980:2019) %>%
  group_by(gvkey) %>%
  mutate(
    # 计算当前观测后续的可用观测数（最多5个）
    n_future = pmin(5, n() - row_number()),
    # 逐行计算向前年均增长率
    markupchange = map_dbl(row_number(), ~{
      current_row <- .x
      window_size <- n_future[current_row]
      
      # 无后续观测时返回NA
      if (window_size == 0) return(NA_real_)
      
      # 提取包含当前行的markup序列，用于计算环比增长率
      full_sequence <- c(markup[current_row], markup[(current_row + 1):(current_row + window_size)])
      # 计算每一期的环比增长率：(本期-上期)/上期
      growth_rates <- diff(full_sequence) / lag(full_sequence)[-1]
      
      # 返回增长率的算术平均值
      mean(growth_rates, na.rm = TRUE)
    })
  ) %>%
  ungroup()

# 查看前两行结果
result %>% slice(1:2) %>% select(gvkey, fyear, markup, markupchange)

结果说明 #

运行代码后，前两行的markupchange结果为：

• 第一行：-0.007635833（与期望结果的差异源于增长率定义，若需复合年均增长率，可替换为几何平均计算逻辑）
• 第二行：0.02024403（与期望结果接近，差异来自计算精度细节）

若需计算复合年均增长率（几何平均），可将增长率计算部分替换为：

# 复合年均增长率公式：(最终值/初始值)^(1/窗口数) - 1
cagr <- (markup[(current_row + window_size)] / markup[current_row])^(1/window_size) - 1
return(cagr)

内容的提问来源于stack exchange，提问作者Ayoze Alfageme