r - 如何根据分组变量对单个序列进行多个序列的回归？

Question

我觉得我的基本问题是如何在单个系列上回归多个系列。尽管我的系列在时间上不相等，但即使我对股票和基准使用相等的时间长度系列(如果需要，我可以提供手动相等的数据)，我也会收到错误。我想估计一个市场模型(即，将所有股票每天的股票返回率与基准返回率进行回归)，并根据长格式的回归制作一个 beta 值的数据框架。因此，对于提供的样本，beta 值数据框中将有 4 个 beta 值(2 个用于 ABC，2 个用于 XYZ)。这是两个股票价格的样本

idf <- structure(list(Firm = c("ABC", "ABC", "ABC", "ABC", "ABC", "ABC", "ABC",
  "ABC", "ABC", "ABC", "ABC", "ABC", "ABC", "ABC", "ABC", "XYZ", "XYZ", "XYZ",
  "XYZ", "XYZ", "XYZ", "XYZ", "XYZ", "XYZ", "XYZ", "XYZ", "XYZ", "XYZ", "XYZ",
  "XYZ"), Date = structure(c(NA, 1451642400, 1451646000, 1451649600, 1451653200,
  1451656800, 1451660400, 1451664000, 1451898000, 1451901600, 1451905200,
  1451908800, 1451912400, 1451916000, 1451919600, NA, 1451642400, 1451646000,
  1451649600, 1451653200, 1451656800, 1451660400, 1451664000, 1451898000,
  1451901600, 1451905200, 1451908800, 1451912400, 1451916000, 1451919600),
  tzone = "UTC", class = c("POSIXct", "POSIXt")), Price = c(1270.9, 1277,
  1273.25, 1273.85, 1273.75, 1272, 1265.35, 1265.35, 1253.1, 1248.1, 1242,
  1248.15, 1241.1, 1246.5, 1242.5, 225.75, 225.7, 225.5, 225.45, 228.6, 227.7,
  227.8, 227.8, 226, 225.1, 222.35, 222.25, 221.1, 221.2, 220.7), rt = c(NA,
  0.00478826614451489, -0.0029408902678254, 0.000471124032113579,
  -7.8505259892836e-05, -0.00137484063686699, -0.00524170116535849, 0,
  -0.00972828256347036, -0.00399808624683118, -0.00489941143098349,
  0.00493947152112462, -0.00566437187907365, 0.00434154082813709, 
  -0.00321414499282824, NA, -0.000221508473604359, -0.000886524880757023, 
  -0.000221754075639957, 0.0138753464936165, -0.00394477829101625, 
  0.000439077943387822, 0, -0.00793305174079517, -0.00399025135959974, 
  -0.0122920309559813, -0.000449842562691316, -0.00518778653061336, 
  0.000452181784779349, -0.00226295638548901), day = structure(c(NA, 16801,
  16801, 16801, 16801, 16801, 16801, 16801, 16804, 16804, 16804, 16804, 16804,
  16804, 16804, NA, 16801, 16801, 16801, 16801, 16801, 16801, 16801, 16804,
  16804, 16804, 16804, 16804, 16804, 16804), class = "Date")),
  .Names = c("Firm", "Date", "Price", "rt", "day"),
  row.names = c(NA, -30L), class = c("tbl_df", "data.frame"))

以及基准指数样本

imdf <- structure(list(Date = structure(c(NA, 1451642400, 1451646000, 1451649600,
  1451653200, 1451656800, 1451660400, 1451664000, 1451898000, 1451901600,
  1451905200, 1451908800, 1451912400, 1451916000, 1451919600, 1451923200),
  class = c("POSIXct", "POSIXt"), tzone = "UTC"), Price = c(3443.1, 3450.5,
  3453.85, 3447.9, 3456.9, 3468.45, 3472.1, 3472.1, 3484.75, 3466.45, 3417.5,
  3416.05, 3401, 3425.75, 3425.25, 3425.25), rt = c(NA, 0.0021469197059254, 
  0.000970402793278424, -0.00172420081111291, 0.00260688364525663, 
  0.00333557458000655, 0.00105178994070698, 0, 0.0036367074012027, 
  -0.00526529010184795, -0.0142217259100423, -0.000424376794422088, 
  -0.00441540679648789, 0.00725091981911596, -0.000145964093093198, 
  0), day = structure(c(NA, 16801, 16801, 16801, 16801, 16801, 16801, 16801,
  16804, 16804, 16804, 16804, 16804, 16804, 16804, 16804), class = "Date")),
  .Names = c("Date", "Price", "rt", "day"), row.names = c(NA, -16L),
  class = c("tbl_df", "tbl", "data.frame"))

以下是我基于 dplyr 的代码，用于通过基准返回率回归日内返回率来估计股票 ABC 和 XYZ 两天的当前贝塔值。

require(dplyr)
q3 <- idf %>%
  group_by(Firm, day) %>%
  summarise(Beta = PerformanceAnalytics::CAPM.beta(idf$rt, imdf$rt)) %>%
  na.omit()

但它给出了以下错误:

Error: The data cannot be converted into a time series. If you are trying to pass in names from a data object with one column, you should use the form 'data[rows, columns, drop = FALSE]'. Rownames should have standard date formats, such as '1985-03-15'.

对于

q3 <- idf %>%
 group_by(Firm, day) %>%
 summarise(Beta = coef(summary(lm(idf$rt~imdf$rt)))[2,1]) %>%
 na.omit()

我收到以下错误:

Error: error in evaluating the argument 'object' in selecting a method for function 'summary': Error in model.frame.default(formula = idf$rt ~ imdf$rt, drop.unused.levels = TRUE) : variable lengths differ (found for 'imdf$rt')`

我还尝试了以下基于聚合函数的操作:

beta.est= function(x) coef(summary(lm(x~imdf$rt)))[2,1]
betas = aggregate(cbind(rt) ~   Firm + day , idf, FUN = beta.est)

但它也给出了同样的错误:

Error in summary(lm(x ~ imdf$rt)) : error in evaluating the argument 'object' in selecting a method for function 'summary': Error in model.frame.default(formula = x ~ imdf$rt, drop.unused.levels = TRUE) :    variable lengths differ (found for 'imdf$rt')

然后我尝试了“for 循环”:

betas=list()
for(i in idf$Firm ){
 for(j in  idf$day){
  betas$day= coef(summary(lm(idf[i,4]|j~imdf$rt|j)))[2,1]}}

我再次明白了:

Error in summary(lm(idf[i, 4] | j ~ imdf$rt | j)) : error in evaluating the argument 'object' in selecting a method for function 'summary': Error in model.frame.default(formula = idf[i, 4] | j ~ imdf$rt | j, drop.unused.levels = TRUE) : variable lengths differ (found for 'imdf$rt | j')

当我对股票和基准使用相同时间长度系列时，出现以下错误:

Error in summary(lm(idf[i, 4] | j ~ imdf$rt | j)) : error in evaluating the argument 'object' in selecting a method for function 'summary': Error in model.frame.default(formula = idf[i, 4] | j ~ imdf$rt | j, drop.unused.levels = TRUE) : variable lengths differ (found for 'imdf$rt | j')

虽然上述代码用于估计当前贝塔值，但如果我也可以回归股票基准指数的滞后性，那就太好了。请帮助我。

Answer 1

我希望你不介意我使用一些额外的软件包。我认为 dplyr、tidyr、purrr 和 broom 的组合确实简化了多重线性回归的工作:

为了确保我们不会将数据列名称与基准列名称混淆，我们只需添加基准名称:

names(imdf) <- paste("bm1", names(imdf), sep ="_")

您还可以添加滞后日期变量和/或加入另一个基准

set.seed(123)
imdf$bm1_Date_lag <- imdf$bm1_Date - 3600     # shift Date by 1 hour
imdf$bm2_rt <- rnorm(nrow(imdf), sd = 0.005)  # add a dummy benchmark

首先，您必须按公司和日对数据进行分组，将其与日期上的基准数据连接并嵌套，这将创建对于每个组，都有一个名为 data 的 data.frame，其中包含所有其他变量。

idf_grouped <- idf %>% group_by(Firm, day) %>%
     left_join(imdf, by = c("Date" = "bm1_Date")) %>% 
     nest()

删除没有 day 值的组，以便您能够在每个组上拟合线性模型，并使用 purrr 中的 map 将其添加为新列。您可以一步对多个模型执行此操作。

idf_grouped <- idf_grouped %>% filter(!is.na(day)) %>%
     mutate(model = map(data,~lm(rt~bm1_rt, data = .)),
            model2 = map(data,~lm(rt~bm2_rt, data = .)))

为了提取线性模型的系数，我更喜欢使用 broom 中的 tidy，因为它返回一个 data.frame，它是 - 在dplyr - 更容易处理。 unnest() 将您的数据转换回单个 data.frame。因为您对拦截不感兴趣，所以我们过滤 bm1_rt

idf_grouped %>% 
    unnest(bm1 = model %>% map(tidy, quick = TRUE),
             bm2 = model2 %>% map(tidy, quick = TRUE),.sep = "_") %>%
    filter(bm1_term == "bm1_rt")
##    Firm        day bm1_term bm1_estimate bm2_term bm2_estimate
##   <chr>     <date>   <fctr>        <dbl>   <fctr>        <dbl>
## 1   ABC 2016-01-01   bm1_rt   0.08879514   bm2_rt   -0.2781275
## 2   ABC 2016-01-04   bm1_rt   0.25888489   bm2_rt    0.3845765
## 3   XYZ 2016-01-01   bm1_rt   0.66791986   bm2_rt   -0.2891714
## 4   XYZ 2016-01-04   bm1_rt   0.45735812   bm2_rt   -0.5014824

要在滞后基准上计算这些 Beta，只需在第一步中加入 "Date"= "bm1_Date_lag" 即可。

编辑:为了实现行业返回，您必须了解每个公司所属的行业。为了便于说明，我只是添加了另一家公司 "DEF" 作为 "XYZ" 的副本，并将其映射到与 "ABC" 相同的行业

idf_2 <- idf %>% filter(Firm == "XYZ") %>% mutate(Firm = "DEF") %>% bind_rows(idf)
firm_map <- data.frame(Firm = c("ABC", "DEF", "XYZ"), Industry = c(1,1,2), stringsAsFactors = FALSE)

只需将其加入idf_2

idf_map <- idf_2 %>% left_join(firm_map, by = "Firm")

并计算例如。各行业平均返回

idf_map %>% left_join(idf_map %>% group_by(Industry, Date) %>% 
    summarise(ind_rt = mean(rt, na.rm = TRUE)), by = c("Industry", "Date"))

现在ind_rt可以用作回归中的解释变量。