[R] Tidyquant 的練習－利用200天平均線的測試

最初大約是3-4月，當時參加了一個 R User group的 Machine Learning Beginner 的課堂，接觸到 "tidyverse" library。Tidyverse" 中的 "dplyr" 在數據整理上十分方便，不過課堂中處理的都是一些非時間序列的數據，好奇在時間序列的股票方面有什麼好工具可以使用呢？那時候開始留意到有一個叫 "tidyquant" 的 R Library，就一直想找時候來學習一下。例如，好不好把之前的股價研究換成用 Tidyquant再試一次？終於，上星期在 Feedly 看到一篇評論信報某專欄的blog post，就正好用這研究來作為練習，有興趣的話建議先看以下兩篇原文：

分析數據還是以數據「作」分析呢？
https://htmichael.blogspot.hk/2017/08/blog-post.html

【EJFQ信析】港股未脫超買 不離三種下場
http://www2.hkej.com/instantnews/market/article/1636469

今次關心的他如何用200天移動平均線去作為後市預測，甚至作為買賣指標。想到的第一個問題，先是「偏離200天移動平均線的程度，與未來20天的平均回報是否負相關？」所以整個步驟分為：
１. 提取股價和求取200天移動平均數值
２. 製作紀錄200天移動平均線偏離度、和未來20天的平均回報的數據列
３. 製作Scatter plot，Correlation、r^2 等用作觀察
４. 製作迴歸線Regression Line、預測值區間Prediction Intervals.


今次數據整理的過程中經常使用Tidyverse中的 mutate()方法－（見：R User Group 的 Machine Learning for Beginners）。R的金融套件以往是quantmod, TTR, 等，但套件一多而各自對資料有不同格式風格的要求，就會混亂和不便，Tidyquant 是一個為了讓 Quants 們也可以融入到 Tidyverse 的格式環境中，今次就用到當中的 tq_get()、tq_mutate()、tq_performance()。tq_get()是用來提取股價的，好比getsymbols()、get.hist.quote()；tq_mutate() 整作指標等的新數據列；tq_performance() 用作分析股票/組合的表現。後面兩個都是調用 TTR、quantmod、tseries 等底層套件的功能，兩者可用的功能可以用 tq_mutate_fun_options()、tq_performance_fun_options() 查詢。

tq_mutate_fun_options() 查看可以進行的數據操作

tq_performance_fun_options() 查看可以進行的表現分析

先看2010年至今的 "偏離200天移動平均線" 和 "隨後20個交易日的平均表現" 的比較。[dSMA] 的是顯著不等於零的( 當 Pr ( > | t | ) 足夠細小)。每比SMA.200 高 1% ，回報就低0.0884% （雖然還是太小了吧？但你不會認為這麼容易可以明顯地賺大錢吧？）。以 31Aug 高約14%為例，預測回報為 -0.884 %，95%區間在大約 +9%至-11% 之間。 區間之闊就是一個簡單線性模型所未考慮的一切因素和隨機誤差。這是否值搏就見仁見智，和你如何控制注碼了。

2010年1月1日 － 2017年8月31日

但另一篇提出一個數據的自我閹割問題，如果用更多數據，加入2010年之前的股價資料，又是否成立呢？改一改開始日期變成2000年，就可以看到下圖，基本上就是一條平線，再沒有斜度的顯著性。Why ?

2000年1月1日 － 2017年8月31日

當再觀看第一張圖的恆指和200天線的走勢，會看到2010年前後的港股表現有所不同：2000-2010年大致上是大升大跌，2010年後就是圍繞200天線起跌上落。所以猜測這是否因為大勢的不同，所以這套方法的效果不同？  回到這個預測方法的原意，它是在相信 Mean Reversion 的觀念；對比另一派是看突破、趁勢。那一個想法才正確，就只有當下的市場才是真理。

如果當偏離200天線後，就會回歸200天線。這就會表現為Scatter plot上的負數斜率；如果愈升(跌)愈高(低)的話，這應該會表現為Scatter plot上的正數斜率，如下圖。

2000年1月1日 － 2010年1月1日

到這裡要問一個重要的問題， 現在（甚至未來）到底是「上落市」還是「單邊市」？我就不知道了，都希望有人可以指教一下。




回到 Tidyquant 的試驗，希望最終是可以建立到買賣紀錄和測試投資組合表現的方式。要以投資組合的角度做比較，而不單是大量個別交易的回報統計，但這個就跟注碼控制有關，每個交易機會時應該把多少資金從Cash轉到Stock？ 雖然這方面要另外測試，但我猜想之下應該會是需要記錄以下資料。

新增說明文字

以下的設定絕對有待優化，簡單地假設期間為2010年1月1日－2017年8月31日，因為要會200天數據作指標，實際在2010年10月22日開始，開始時有23517元現金（當時的恆指水平），當 dSMA < 15% 買入 1個單位，設為點數的十分之一，約二千多元，10天後放出（dSMA > 15% 賣入 1個單位，10天後買回），如果認真操作的話應該有更仔細的退出策略。表現如下圖藍線。能夠以Portfolio的⻆度考慮的話，甚至可以用 Tidyquant 中的 tq_performance() 計算：Alpha、Beta、Sharpe Ratio等：

另一方面，只統計大量單筆交易的結果還算容易，只要看每次從觸發交易到退出期間的回報做個總結。以下可見平均是有約3.6%的回報：

R Code：

## Example of 200 Day Testing ##
library(tidyquant)   # 引入 tidyquant library

dateStart <- as_date("2010-01-01")  # 設定觀察期初
dateEnd <- as_date("2017-08-31")  # 設定觀察期尾

# 用tidyquant的 tq_get 取得 tibble 格式的股價資料
StockHSI <- c("^HSI") %>% 
    tq_get(get = "stock.price",
    from = dateStart, 
    to = dateEnd)
myData <- StockHSI %>% filter(!is.na(close))  # 觀察過資料後，留意到幾筆 NA 資料 (例如23/08的颱風停市)，直接刪除

# 用 geom_ma 製作 一個 20天 和 200天 移動平均線 的圖
p1 <- myData %>%
    ggplot(aes(x = date, y = close)) +
    geom_barchart(aes(open = open, high = high, low = low, close = close),
        color_up = "darkgreen", color_down = "darkred") +
    geom_ma(ma_fun = SMA, n = 200, col = "darkblue", size = 0.5) +
    geom_ma(ma_fun = SMA, n = 20, col = "blue", size = 0.5) 
p1

# 自定一個函數，計算未來20天的平均回報
avgReturn <- function(x, p=1, ... ){
   V <- matrix(NA, ncol=p, nrow=NROW(x))
   for (i in 1:p) { V[,i] <- lead(x, n=p)/x }
   avgReturn <- apply(V,1,function(x) sum(x)) / p
   return(avgReturn)
}

# 用 tq_mutate、mutate 的製作紀錄200天移動平均線偏離度、和20天內平均回報 的數據列
myData <- myData %>% tq_mutate(select=close, mutate_fun = SMA, n = 200) %>% rename(SMA200 = SMA)
myData <- myData %>% mutate(dSMA = close / SMA200 -1, 
    avgRet = avgReturn(close,20) - 1)

# 記下當日的偏離度，之後去除NA資料
todayData <- myData %>% filter(date==dateEnd) 
myData <- myData %>% filter(!is.na(avgRet) & !is.na(dSMA)) 

# 製作 200天移動平均線偏離度 的圖
p2 <- myData %>%
    ggplot(aes(x = date, y = dSMA)) +
    geom_line()
p2

# 將數據用作線性回歸模型
lm1 <- lm(avgRet ~ dSMA, data = myData)
summary(lm1)
pAvgRet <- predict(lm1, interval="prediction")
myData <- cbind(myData, pAvgRet)

# 製作 Scatter Plot，加上regression line、 預測區間、最近的200天偏離度。
ggplot(myData, aes(x = dSMA, y = avgRet)) + 
  geom_point(size=0.5) +
  geom_vline(xintercept=todayData$dSMA, col="grey", size=0.5, linetype="dashed") +
  stat_smooth(method = "lm", col = "red", size=0.75) +
  geom_line(aes(y=lwr), color = "red", size=0.5, linetype = "dashed") +
  geom_line(aes(y=upr), color = "red", size=0.5, linetype = "dashed") +
  labs(x = "Deviation from SMA.200", y = "Average Return.20") 


# 紀錄觸發點
myData <- myData %>% mutate(trigger = ifelse(dSMA>0.15, -1, ifelse(dSMA < (-0.15), 1, 0) ))

# 每個機會出現後20天的回報總結:
myData %>% filter(trigger == 1) %>% select(avgRet)  %>% summary() # 買入信號
myData %>% filter(trigger == -1) %>% select(avgRet)  %>% summary() # 賣入信號

# 自定一個函數，將交易紀錄轉變成投資組合的價值
trade2Hold <- function(price, trade, initCash=100, initStock=0, ... ){
   V <- matrix(NA, nrow=length(trade), ncol=6)
   V[1,1] = price[1]; # price
   V[1,2] = trade[1]; # trade
   V[1,3] = initStock + trade[1]; # Holding Stock Number
   V[1,4] = V[1,1] * V[1,3]; # Holding Stock Value
   V[1,5] = initCash - V[1,1] * V[1,2]; # Holding Cash Value
   V[1,6] = V[1,4] + V[1,5]   
   for (i in 2:length(trade)) {
       V[i,1] = price[i]; # price
       V[i,2] = trade[i]; # trade
       V[i,3] = V[i-1,3] + trade[i]; # Holding Stock Number
       V[i,4] = V[i,1] * V[i,3]; # Holding Stock Value
       V[i,5] = V[i-1,5] - V[i,1] * V[i,2]; # Holding Cash Value
       V[i,6] = V[i,4] + V[i,5]      
   }
   return(V[,6])
}

# 製作買賣紀錄 和 投資組合價值 的數據列
initialCash = myData$close[1] # 將開始的現金定為恆指的期初點數，方便繪圖比較myData <- myData %>% mutate(trade = 0.1*(trigger - lag(myData$trigger, n=10, default=0))) myData <- myData %>% mutate(pfValue = trade2Hold(close, trade, initCash= initialCash))   # 組合價值
myData <- myData %>% tq_mutate(select=close, mutate_fun=dailyReturn, col_rename="bmRet") %>% 
          tq_mutate(select=pfValue, mutate_fun=dailyReturn, col_rename="pfRet") #  組合和基準的回報

# 製作 組合 和 基準 的圖
p3 <- myData %>% ggplot(aes(x = date)) +
    geom_line(aes(y = pfValue), col = "blue", size = 0.5) +
    geom_line(aes(y = close), col = "black", size = 0.5) +
    geom_line(aes(y = SMA200), col = "black", size = 0.1, linetype = "dotdash")
p3

# CAPM's Alpha, Beta;  Sharpe Ratio
myData %>% tq_performance(Ra = pfRet, 
                   Rb = bmRet, 
                   performance_fun = table.CAPM)

myData %>% tq_performance(Ra = pfRet, 
                          Rb = NULL, 
                          performance_fun = SharpeRatio,
                          Rf=0.03/365,
                          p=0.99)