我自己想出来的。代码如下：

#Split df into list of df's by permno
df.list <- split(Merged_df , f = Merged_df$permno )

# Function that loops through list of dfs and calculates yearly 
# covariance between return and narratives
cov_df_wide_list[1:2] <- lapply(df.list, function(df) {
  # create a new column with year information
  df <- df %>%
    mutate(year = lubridate::year(date))

  # group the dataframe by year and each parameter column
  # and calculate the covariance with the return column
  cov_df <- df %>%
    select(year, permno, Natural.disasters:Revenue.growth, ret) %>%
    gather(key = "param_col", value = "value", -year, -permno, -ret) %>%
    group_by(year, permno, param_col) %>%
    summarize(cov = cov(value, ret))

  # spread the dataframe to get each year's covariances
  cov_df_wide <- cov_df %>%
    spread(key = "param_col", value = "cov")

  return(cov_df_wide)
})

基本上，该函数循环遍历一个 Dataframe 列表，并计算我的返回列和每个风险参数之间的协方差。