35行代码搞定事件研究法（下）
作者簡介：
祝小宇，個人公眾號：大貓的R語言課堂
前文推送：
35行代碼搞定事件研究法（上）
Hello親愛的小伙伴們，上期已經講到如何對單一事件日計算超額收益，本期將會教大家如何針對多個股票多個事件日計算超額收益，Let's go!
注意 I，本代碼主要使用data.table包完成，關于data.table包的相應知識會在涉及的時候進行講解。在以后的課堂中，我們會重點介紹data.table這個包。
注意 II， 本代碼還使用了partial()函數，它來自于pryr這個包
用data.table包處理多個事件日?
本期課堂的核心代碼只有下面5行（應用了data.table包的語法）：
雖然看起來似乎有些難懂，但如果我們將他分解為三部分，理解起來就容易多了。
首先，這5行代碼可以抽象為如下形式：
其中，event數據集就是我們在上節課講到的包含有股票代碼、日期、股票收益率、市場收益率、事件日標識的數據集（什么你忘了？快去看上節課教程！就是那個黑色的圖）。請觀察在上面這個抽象后的代碼，大家應該可以看出我們對event數據集做了三件事情，具體分別為：
選取event中所有的行（第一行代碼）。此處，我們沒有添加任何條件，因此默認選中event的所有行。
對選中的變量進行操作（第二行代碼）。此處，所有的操作都用大括號{}包裹了起來。
對event按照stk.id進行分組（第三行代碼）。加了這一行代碼后，第二行代碼中所有的操作都會對每個stk.id分組運行一遍（這一步很關鍵！）。
講到這，大家一定會發現，上述代碼的關鍵部分就在大括號{...}所括起來的內容。我們一行一行來看：
這一行代碼的作用找到每個股票的所有事件日的序號 ns。大家應該還記得在上一講中我們用 n 來表示單一事件日的序號吧？在這里，which(event.flg == 1)的意思是返回所有event.flg變量等于1的那些行的序號，很自然的，在這里 ns 應該是一個向量。
在上一講中，我們已經給出了函數 do_car() 用來求單個事件日的超額收益，因此很自然的，我們希望對于事件日向量 ns 中的每個元素，都應用一遍 do_car()這個函數。為了做到這一點，我們運用了lapply() 函數。因此代碼就變成了
那么，在最初給的那段代碼中，partial()函數是用來干什么的呢？在這里我們不妨先回憶一下上一講中的do_car() 函數有哪些參數：
}
看到了沒有？do_car() 要求我們提供n, r, rm, date 四個參數，但是向量 ns 只能提供 n 這一個參數的值，因此我們需要用pryr包中的partial() 函數把剩下的幾個變量補充完整（感謝pryr的作者Hadley！如果不是你，我們需要寫許多非常冗長的代碼）。
最后，將處理的結果賦值給car，我們的任務就完成了！下圖是最終的輸出結果（部分）：
其中，stk.id是股票代碼，date是事件日（非事件日不輸出結果），coef是該事件日對應的收益率模型的系數（alpha、beta），ars是對應的超額收益率。在我們的例子中，我們只計算T日前后各一日的收益，因而ars一共有三個元素。
到此為止，求超額收益的計算就完成了，現在大貓給出完整的代碼：
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 ") ? ? } else if (n + c2 > length(r)) { ? ? ? ? cat("n =", n, "is too large ") ? ? } 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萬個左右的事件。一般的事件研究法的數據量極少超過這個量級。
在這里附上生成模擬數據集的代碼：
我們接著設定 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 :=” 表示在原數據集中新建一個名為 car 的變量，vapply(ars, sum)的含義是把超額收益率向量ars中的元素相加，double(1)指定輸出的必須是一個標量（因為對于每個事件日，CAR是唯一的）
再比如，如果我們想計算逐日的累計超額收益率，那么代碼就為：
cumsum() 是累計求和函數。注意，此時最終得到的cunsum應該是一個和ars長度相等的向量。
如果我們希望對每個股票的CAR進行T檢驗，那么代碼就為：
最終的結果為：
其中，t.test給出了 t 值，p.ttest 給出了對應的p值。
其實，還有很多別的后續工作可以擴展，大貓就不一一介紹啦，小伙伴們可以自行實驗。最后，如果覺得大貓的R語言課堂有用，請多多支持關注哦！
公眾號后臺回復關鍵字即可學習
回復?爬蟲????????????爬蟲三大案例實戰
回復?Python???????1小時破冰入門
回復?數據挖掘?????R語言入門及數據挖掘
回復?人工智能?????三個月入門人工智能
回復?數據分析師??數據分析師成長之路?
回復?機器學習?????機器學習的商業應用
回復?數據科學?????數據科學實戰
回復?常用算法? ? ?常用數據挖掘算法
總結
以上是生活随笔為你收集整理的35行代码搞定事件研究法（下）的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 计算机电源选平衡,电脑里选择电源计划哪个
- 下一篇： 你的每行代码值多少钱？