35行程式碼搞定事件研究法（下）

作者簡介：

祝小宇，個人公眾號：大貓的R語言課堂

前文推送：

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 <- 1m1 <- 10m2 <- 5
# do cardo_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萬個左右的事件。一般的事件研究法的資料量極少超過這個量級。

在這裡附上生成模擬資料集的程式碼：

n.day <- 1000n.stk <- 2500p <- 0.02stk.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語言課堂有用，請多多支援關注哦！

公眾號後臺回覆關鍵字即可學習

回覆 爬蟲            爬蟲三大案例實戰
回覆 Python       1小時破冰入門
回覆 資料探勘     R語言入門及資料探勘
回覆 人工智慧     三個月入門人工智慧
回覆 資料分析師  資料分析師成長之路 
回覆 機器學習     機器學習的商業應用
回覆 資料科學     資料科學實戰
回覆 常用演算法     常用資料探勘演算法