微信号:gh_a9d9c3cb1ca5

介绍:市面上的教材太水?学了半天还是不知道如何处理工作中的数据?大猫的R语言课堂每期都和你分享一个工作研究中的实际问题,期期干货,就是这么实在!

35行代码搞定事件研究法(下)

2016-04-06 11:21 大猫
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语言课堂有用,请多多支持关注哦!

下期预告
第一讲就告一段落,我是大猫,最后放一个次回预告! 如何优雅地求移动平均!


 
大猫的R语言课堂 更多文章 大猫的R语言课堂开课啦+次回预告 35行代码搞定事件研究法(上) 如何用R进行中文分词? 经验总结 | 最有效的R学习路径(二) R Tricks: 如何处理Gaps &amp; Islands问题?
猜您喜欢 R语言股指波动分析——股票系列(一) [ECUG专题回顾]《深入理解容器技术》-田琪(京东资深架构师) 怎样交付业余项目 首档大数据对话网络节目《职话大数据》在华信智原(北京)录制