35行代码搞定事件研究法（下）

Hello亲爱的小伙伴们，上期已经讲到如何对单一事件日计算超额收益，本期将会教大家如何针对多个股票多个事件日计算超额收益，Let's go!

注意 I，本代码主要使用data.table包完成，关于data.table包的相应知识会在涉及的时候进行讲解。在以后的课堂中，我们会重点介绍data.table这个包。

注意 II， 本代码还使用了partial()函数，它来自于pryr这个包

用data.table包处理多个事件日

本期课堂的核心代码只有下面5行（应用了data.table包的语法）：

> car <- event[, {

> ns <- which(event.flg == 1);

> lapply(ns, partial(do_car, r = r, rm = rm, date = date)) %>% rbindlist()

> },

> by = stk.id]

虽然看起来似乎有些难懂，但如果我们将他分解为三部分，理解起来就容易多了。

首先，这5行代码可以抽象为如下形式：

> event[,

> {...},

> by = stk.id]

其中，event数据集就是我们在上节课讲到的包含有股票代码、日期、股票收益率、市场收益率、事件日标识的数据集（什么你忘了？快去看上节课教程！就是那个黑色的图）。请观察在上面这个抽象后的代码，大家应该可以看出我们对event数据集做了三件事情，具体分别为：

选取event中所有的行（第一行代码）。此处，我们没有添加任何条件，因此默认选中event的所有行。

对选中的变量进行操作（第二行代码）。此处，所有的操作都用大括号{}包裹了起来。

对event按照stk.id进行分组（第三行代码）。加了这一行代码后，第二行代码中所有的操作都会对每个stk.id分组运行一遍（这一步很关键！）。

讲到这，大家一定会发现，上述代码的关键部分就在大括号{...}所括起来的内容。我们一行一行来看：

ns <- which(event.flg == 1);

这一行代码的作用找到每个股票的所有事件日的序号 ns。大家应该还记得在上一讲中我们用 n 来表示单一事件日的序号吧？在这里，which(event.flg == 1)的意思是返回所有event.flg变量等于1的那些行的序号，很自然的，在这里 ns 应该是一个向量。

在上一讲中，我们已经给出了函数 do_car() 用来求单个事件日的超额收益，因此很自然的，我们希望对于事件日向量 ns 中的每个元素，都应用一遍 do_car()这个函数。为了做到这一点，我们运用了lapply() 函数。因此代码就变成了

lapply(ns, do_car)

那么，在最初给的那段代码中，partial()函数是用来干什么的呢？在这里我们不妨先回忆一下上一讲中的do_car() 函数有哪些参数：

do_car <- function(n, r, rm, date) {

....

}

看到了没有？do_car() 要求我们提供n, r, rm, date 四个参数，但是向量 ns 只能提供 n 这一个参数的值，因此我们需要用pryr包中的partial() 函数把剩下的几个变量补充完整（感谢pryr的作者Hadley！如果不是你，我们需要写许多非常冗长的代码）。

最后，将处理的结果赋值给car，我们的任务就完成了！下图是最终的输出结果（部分）：

其中，stk.id是股票代码，date是事件日（非事件日不输出结果），coef是该事件日对应的收益率模型的系数（alpha、beta），ars是对应的超额收益率。在我们的例子中，我们只计算T日前后各一日的收益，因而ars一共有三个元素。

完整的代码

到此为止，求超额收益的计算就完成了，现在大猫给出完整的代码：

c1 <- 1

c2 <- 1

m1 <- 10

m2 <- 5

# do car

do_car <- function(n, r, rm, date) {

stopifnot(m1 > m2)

if (n - m1 < 0) {

cat("n =", n, "is too small \n")

} else if (n + c2 > length(r)) {

cat("n =", n, "is too large \n")

} else {

i1 <- max(1, n - m1)

i2 <- n - m2

i3 <- n - c1

i4 <- n + c2

r.model <- r[i1:i2]

rm.model <- rm[i1:i2]

r.car <- r[i3:i4]

rm.car <- rm[i3:i4]

model <- lm(r.model ~ I(r.model - rm.model))

coef <- coef(model)

ars <- r.car - predict(model, list(r.model = r.car, rm.model = rm.car))

list(date = date[n], coef = list(coef), ars = list(ars))

}

}

car <- event[, {

ns <- which(event.flg == 1);

lapply(ns, partial(do_car, r = r, rm = rm, date = date)) %>% rbindlist()

},

by = stk.id]

最上面三行注释用来描述数据结构，如果去掉的话，所有代码加起来35行都不到，是不是很神奇！

性能测试

大猫在这里给出的代码已经经过高度优化，是在尝试众多可行方法后给出的计算速度最快的版本。小伙伴大可不必担心自己的数据太多计算机跑不起来。但是口说无凭，大猫在这里给出用模拟数据得到的测试结果。

在测试中，大猫设置了一个极端条件：模拟2500个股票（差不多是A股股票数），每个股票拥有1000个交易日的记录（差不多有4年的时间），平均50个交易日出现一个事件（模拟盈利公告这类事件的出现频率）。因此在整个数据集中，一共有250万条观测，5万个左右的事件。一般的事件研究法的数据量极少超过这个量级。

在这里附上生成模拟数据集的代码：

n.day <- 1000

n.stk <- 2500

p <- 0.02

stk.id <- rep(1:n.stk, each = n.day)

event.flg <- rbinom(n.day * n.stk, 1, p)

date <- rep(seq(from = as.Date("2000-01-01"), by = "day", length.out = n.day))

event <- data.table(stk.id, date = rep(date, n.stk), r = runif(n.stk * n.day), rm = runif(n.stk * n.day),

event.flg = event.flg)

我们接着设定 c1 = c2 = 1, m1 = 90, m2 = 5，也即求 [T - 1, T + 1] 区间的超额收益，并用 [T - 90, T - 5] 这个区间估计收益率模型。这也是一个比较常见的设定。

大猫用这个数据集在自己的surface pro 4 i5版上连续跑了三遍，每一次的耗时分别为：

79s 81s 82s

三次平均耗时在80秒左右。可以说，这是一个非常优秀的成绩了。况且我们平时遇到的数据集应该远远小于模拟数据集，小伙伴还担心什么嗯？

对CAR进行 T 检验

既然已经算出了超额收益AR，那么下面我们自然希望把AR加起来得到累计超额收益CAR并进行T检验。例如，我们想知道每个事件日对应的累计超额收益，那么代码就为：

car[,

car := vapply(ars, sum, double(1))

]

其中，car数据集是上面计算得到的所有事件日对应的超额收益率。语句“car :=” 表示在原数据集中新建一个名为 car 的变量，vapply(ars, sum)的含义是把超额收益率向量ars中的元素相加，double(1)指定输出的必须是一个标量（因为对于每个事件日，CAR是唯一的）

再比如，如果我们想计算逐日的累计超额收益率，那么代码就为：

car[,

cumcar := lapply(ars, cumsum)

]

cumsum() 是累计求和函数。注意，此时最终得到的cunsum应该是一个和ars长度相等的向量。

如果我们希望对每个股票的CAR进行T检验，那么代码就为：

ttest <- car[,

.(t.test = sapply(ars, function(x) t.test(x)$statistic),

p.ttest = sapply(ars, function(x) t.test(x)$p.value)),

by = .(stk.id)

]

最终的结果为：

其中，t.test给出了 t 值，p.ttest 给出了对应的p值。

其实，还有很多别的后续工作可以扩展，大猫就不一一介绍啦，小伙伴们可以自行实验。最后，如果觉得大猫的R语言课堂有用，请多多支持关注哦！