多重滚动回归（Fama-Macbeth 过程）

Question

问题介绍

大家好，

我目前正在使用著名的 Fama-Macbeth 滚动窗口程序来完成我的论文，以估计具有许多自变量和不平衡面板数据的模型。将与不同的多因素模型进行比较。在这里我要求完整的一个（8 个因素），因为完整的代码将用作较小模型的参考。

问题就出现在这里，因为网络上的许多教程都处理滚动回归，但很少有人以多元回归的方式展示它。如果没有足够的对前面示例的参考，我的编码知识目前还不足以解决这样的问题。 此外，如果有人更了解 Python，我也非常感谢此代码的 Python 版本。

我研究了许多类似的问题，但除了一个，我找不到适合我的数据的东西。似乎最快的方法之一是使用矩阵代数或“rollRegres”。好吧，“rollRegres”版本可以正常运行，直到我的数据显示缺失值（空白单元格）自动转换为“NA”。但是，这些缺失的数据显示不同证券的发行日期不同，这意味着不应对这些缺失的单元格采取任何行动。

因此，如果有人知道如何加速此代码（在较大数据集的情况下远非高效），提出其他建议，或正确地向我展示如何以矩阵形式执行，我将非常感激。

总结：

• 滚动回归（8 ind. + 1300 dep. var's with window 60 days） - 需要测试版

  • 大型数据集“Y”（1.300 个深度变量，每个变量包含 1000 个条目）
  • 数据集“X”（8 个独立变量 - 因子 - 全部包含 1000 个条目）

• 使其健壮的任何可能性（将 lm() 替换为库 MASS 的 rlm() 没有成功）
• 由于发布日期不同而导致数据集“Y”中的 NA（数据特殊性）
• 效率（速度问题）
• 输出（导出问题）

代码

我的数据集的公平再现如下。请记住，原始数据集要大得多，这样的代码效率远非如此。

###LIBRARIES
library(zoo)


###DATASET

set.seed(10000000)

df <- data.frame(

##Date
Date = seq(as.Date("2000/01/01"), by = "day", length.out = 5000),  

##Factors (indep var)
MKTRF = sample(1:100, 5000, replace=TRUE),
SMB = sample(1:100, 5000, replace=TRUE),
HML = sample(1:100, 5000, replace=TRUE),
RMW = sample(1:100, 5000, replace=TRUE),
CMA = sample(1:100, 5000, replace=TRUE),
MOM = sample(1:100, 5000, replace=TRUE),
TERM = sample(1:100, 5000, replace=TRUE),
DEF = sample(1:100, 5000, replace=TRUE),

##Securities (dep var)
#type 0
B1 = c(rep(NA, 1000), sample(1:100, 3000, replace=TRUE), rep(NA, 1000)),
B2 = c(rep(NA, 1500), sample(1:100, 3000, replace=TRUE), rep(NA, 500)),
B3 = c(sample(1:100, 3000, replace=TRUE), rep(NA, 2000)),
B4 = c(rep(NA, 1000), sample(1:100, 4000, replace=TRUE)),
B5 = c(sample(1:100, 5000, replace=TRUE)),

#type 1
GB1 = c(rep(NA, 1500), sample(1:100, 2000, replace=TRUE), rep(NA, 1500)),
GB2 = c(rep(NA, 1500), sample(1:100, 3000, replace=TRUE), rep(NA, 500)),
GB3 = c(sample(1:100, 3000, replace=TRUE), rep(NA, 2000)),
GB4 = c(rep(NA, 1000), sample(1:100, 4000, replace=TRUE)),
GB5 = c(sample(1:100, 5000, replace=TRUE)))


###REGRESSION

#Setting the problem
z <- read.zoo(df, FUN = as.yearmon, format = "%Y/%m/%d")
View(z)

getCoef <- function(z, lhs, rhs) {
  if (all(is.na(z[, lhs]))) "Empty"         #"Empty" is for me to spot empty cells due to different issue dates
  else coef(lm(paste(lhs, "~", rhs), z))
}

roll <- function(z, lhs, rhs = "MKTRF + SMB + HML + RMW + CMA + MOM + TERM + DEF") {
  rollapplyr(z, 60, getCoef, by.column = FALSE, coredata = FALSE, lhs = lhs, rhs = rhs)
}

ynames <- c("B1", "B2", "B3", "B4", "B5", "GB1", "GB2", "GB3", "GB4", "GB5")

#Regress
L <- lapply(ynames, roll, z = z)

#Save outputs
output <- Map(fortify.zoo, L)

此代码来自 Stack-Overflow 上的一个类似问题，可以运行，但在输出中显示了一个小问题（屏幕截图将变得清晰 - 金字塔形结果显示在第一行数据中）。此外，它对于长数据集绝对没有效率。最后，我想以简洁明了的方式导出每种证券的所有测试版（在本例中为 B1:GB5）；现在这些保留在“输出”列表中；最好是在一个带有注释的表格中打开。

如果有人愿意提出任何类型的改进建议，我们将不胜感激。此外，我认为这对于其他处理此类回归的人来说可能仍然是一个可靠的参考。

Image - Pyramidal_outputs

提前致谢！

Answer 1

我们可以做一些改进：

• 使用矩阵 X 和 Y
• 对于每个子集，我们可以一次回归所有 Y 列，从而将回归函数的调用次数减少一个数量级
• 处理 NA 错误
• 滚动应用索引而不是数据
• 使用 .lm.fit 代替 lm
• 将结果重塑为 3d 数组

这会产生以下代码。

library(zoo)

nr <- nrow(df)
X <- cbind(1, as.matrix(df[2:9]))
Y <- as.matrix(df[10:19])

coef_fun <- function(ix) {
  co <- try(coef(.lm.fit(X[ix, ], Y[ix, ])))
  if (inherits(co, "try-error")) NA else c(co)
}
out <- rollapplyr(1:nr, 60, coef_fun, fill = NA)

out <- array(out, c(nr, ncol(X), ncol(Y)), 
  dimnames = list(NULL, colnames(X), colnames(Y)))