35行程式碼搞定事件研究法(下)

R語言中文社群發表於2018-07-25

640?wx_fmt=gif


作者簡介:

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


前文推送:

35行程式碼搞定事件研究法(上)


Hello親愛的小夥伴們,上期已經講到如何對單一事件日計算超額收益,本期將會教大家如何針對多個股票多個事件日計算超額收益,Let's go!


640?wx_fmt=png

注意 I,本程式碼主要使用data.table包完成,關於data.table包的相應知識會在涉及的時候進行講解。在以後的課堂中,我們會重點介紹data.table這個包。


注意 II, 本程式碼還使用了partial()函式,它來自於pryr這個包

640?wx_fmt=png


用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,我們的任務就完成了!下圖是最終的輸出結果(部分):

640?wx_fmt=jpeg

其中,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)]


最終的結果為:

640?wx_fmt=jpeg

其中,t.test給出了 t 值,p.ttest 給出了對應的p值。


其實,還有很多別的後續工作可以擴充套件,大貓就不一一介紹啦,小夥伴們可以自行實驗。最後,如果覺得大貓的R語言課堂有用,請多多支援關注哦!


640?wx_fmt=jpeg

640?wx_fmt=jpeg

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

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

相關文章