《巨量資料的第一步－R語言與商業應用》
《實戰Java－9個別具特色的實作經驗》
《Arduino錦囊妙計》
《Raspberry Pi 機器人自造專案》
這本學習R的書，算是我看過的R中文書當中很好的一本，特別是案例的部份的學習價值就很高。基本內容都有提及一些我不熟悉的。例如時間和時間序列的類型、資料連接SQL資料庫、處理遺留數據等。
上月初，周末加班時就帶著電腦和書，在坐車時跟著書去嘗試。誰知把書放在椅背休息一下，落車時就忘了帶走。好在一星期後發現有人已經代為把書還了圖書館，感謝這個好心人。
首先跟著這本書試的是這一章「股票市場的預測」。 它是用quandmod套件中的getSymbols()功能獲得xts時間序列類型的股價資訊。建立一個自行定義的回報觀察值去衡量價格變動，再用常見的技術分析指標作變數，建立決策樹Decision Tree的模型。最後評價模型的預測誤準確度。
按書中的方式和參數所建立的算是一個短線投資的指標（T.Index），因為觀察的是日子t當日後的十天內的回報情況。如果T值大於某個臨界點a1則該買入，小於某個臨界點a2則該賣出。
$$P_t = \frac{C_t + H_t + L_t}{3}; V_{(t,i)} = \frac{P_{t+i} - C_t}{C_t}$$
$$T.Index_t = \sum_{i=1}^{k=10} (V_{t+i}:|V_{t+i}| > p ) $$
p值的設立，是要只考慮顯著變動的回報。
而 T.index（觀察值, 或預測值）又透過所建立的模型而與一堆指標（自變數）所相關。所用的自變數包括10天內的價格變動、及其移動平均數和標準差，收市價的ADX, SMI, MACD,....等10個可以用TTR套件得到的技術指標，參數都只是用default的設定，可以參考文件(https://cran.r-project.org/web/packages/TTR/TTR.pdf)。
- ATR (Average True Range)
- SMI (Stochastic Momentum Index)
- ADX (Average Directional Index)
- Arron
- Blooinger Bands
- Chaikin Volatility
- CLV (Close Location Value)
- MACD
- MFI (Money Flow Index)
- SAR (Parabolic-Stop and-Reverse)
書中用的是上證指數“^SSEC”在 1990-01-01 至 2013-06-30期間的價格。2013年之前的是 Training Set，2013年開始的6個月是Testing Set。準確度有86.36%。－ $ = \frac{(2)+(11)}{(8+2+0)+(0+1+11)}$
$$ Accuracy = \frac{ (true.Buy | pred.Buy) + (true.Sell | pred.Sell) }{( \ldots | pred.Buy) + ( \ldots | pred.Sell)}$$
可惜，這個高準確度源自一個coding上的錯誤，令應該是Testing Set的數據也包括了在Training Set 當中。一旦修改過後，準確度下跌至25%。
Testing Set: 01Jan2013 - 30Jun2013:
| Training Set: 01Jan1990 - 30Jun2013 | Training Set: 01Jan1990 - 31Dec2012 |
|---|---|
 |
 |
加入一些改動後，改用港股的：恆指、匯控、新地、港鐵、大家樂、澳博、騰訊。結果如下:
Training Set: 1990-01-01 至 2013-12-31；
Testing Set: 2014-01-01 至 2016-09-30：
| myStock: | ^HSI | 0005.HK | 0016.HK | 0066.HK | 0341.HK | 0880.HK | 0700.HK |
|---|---|---|---|---|---|---|---|
| Accuracy (%): | 32.94 | 34.93 | 31.50 | 27.33 | 28.52 | 36.90 | 36.54 |
以恆指為例，實際分類結果：
| ^HSI |
|---|
 |
不過，這個準確度的定義也可以斟酌，因為所謂的失誤大多是預計上升(下跌)後，實際只是持平，而非相反方向。可能在調試參數後會有較好的分辨能力。還有是否該引入考慮回報值而不單單是估中升跌的次數？短線投資的benchmark應該是什麼？而且，這本書的角度也是電腦軟件R 和 建立模型過程 上的指導，有了這個框架後，實際的回報還是自行研究吧。
源碼的內容如下，當中略作了修改以方便在港股上使用。
1) 例如"2800.HK" 的格式：GetSymbols() 在命名變數時，港股以數字開頭和當中的點號，都令到調用變數時做成障礙。所以要指定變數和資料列的名字
2) 可以改用Yahoo下載後的檔案作匯入。
3) 加入levels 確保在單邊市下原本的準確度的函數仍能正常運作
4) 加入了as.Date()，去修正錯誤的Training Set範圍
5) 刪去了p=0.025的條件
另外幾篇包括1)「鳶尾花的分類」：用多個不同的方法去建立分類模型，再比較一不同模型的準確度。2)「關聯分析」處理遺留的數據，使用Apriori演算法。3)「推薦系統」使用RODBC套件連接mySQL，衡量消費者之間／物品之間的相似度，並實現推薦。這個mySQL的連接是另一個很期待的實例，可惜原來Mac上使用RODBC有點麻煩，要摸索另一個RJDBC的套件。
R Code:
# Import Libraries
library(quantmod); library(PerformanceAnalytics);library(tseries); library(e1071)
#% Get Price from Yahoo (package:quantmod return xts))
myStock <- getSymbols(Symbols = "^SSEC", from="1990-01-01", to="2013-6-30", src="yahoo", auto.assign=FALSE)
myStock <- getSymbols(Symbols = "^HSI", from="2000-01-01", to="2016-9-30", src="yahoo", auto.assign=FALSE)
colnames(myStock) = paste("myStock", c("Open","High","Low","Close","Volume","Adjusted"), sep=".")
#% Get 2800.HK Price from downloaded Yahoo file
# ColClasses = c("character","numeric","numeric","numeric","numeric","numeric","numeric")
# ColNames = c("Date", paste("myStock", c("Open","High","Low","Close","Volume","Adjusted"), sep="."))
# myStock = read.csv("HK2800.csv", skip=0, sep=",", colClasses=ColClasses, col.names=ColNames, header=TRUE)
# rm(ColClasses, ColNames)
# myStock <-transform(transform(myStock, Date=as.Date(Date, format="%Y-%m-%d"))[order(myStock $Date),] )
# myStock <- as.xts(myStock[,-c(1)], order.by= myStock $Date)
p <- 0.025
k <- 10
#% Define function as for measuring performance
T.index <- function(data,p,k){
require(quantmod)
hlc = HLC(data)
P <- rowMeans(hlc)
V <- matrix(NA, ncol=k, nrow=NROW(P))
for (i in 1:k) {
V[,i] <- Next(Delt(P,k=i),k=i)
}
# T <- apply(V,1,function(x) sum(x[abs(x)>p]))
T <- apply(V,1,function(x) sum(x,na.rm=TRUE))
T <- xts(x=T, order.by=time(data))
return(T)
}
# Define function as indicators for prediction
myTTR <- function(data){
require(TTR)
require(quantmod)
names(data) <- sapply(X=names(data), FUN=function(x) strsplit(x, split=".", fixed=TRUE)[[1]][2]) # change [2]
myATR <- ATR(HLC(data))$atr
mySMI <- SMI(HLC(data))$SMI
myADX <- ADX(HLC(data))$ADX
myAroon <- aroon(HLC(data)[,-3])$oscillator
myBBands <- BBands(HLC(data))$pctB
myChaikin <- Delt(chaikinVolatility(HLC(data)[,-3]))[,1]
myCLV <- EMA(CLV(HLC(data)))[,1]
myMACD <- MACD(data[,"Close"])[,2]
myMFI <- MFI(HLC(data), data[,"Volume"])
mySAR <- SAR(data[,c("High", "Close")])[,1]
result <- cbind(myATR, mySMI, myADX, myAroon, myBBands, myChaikin, myCLV, myMACD, myMFI, mySAR)
colnames(result) <- cbind("myATR", "mySMI", "myADX", "myAroon", "myBBands", "myChaikin", "myCLV", "myMACD", "myMFI", "mySAR")
return(result)
}
% Variable Selection 1 - randomForest
require(randomForest)
set.seed(42)
model <- specifyModel(formula=T.index(myStock, p=0.025, k=10) ~ Delt(Cl(myStock), k=1:10) + myTTR(myStock) + runMean(Cl(myStock)) + runSD(Cl(myStock)) ,na.rm=TRUE)
rf <- buildModel(x=model, method="randomForest", training.per=c(start(myStock), index(myStock["2015-01-05"])))
varImpPlot(rf@fitted.model)
x <- importance(rf@fitted.model)
rownames(x)[x>12]
# Variable Selection 2 - refined model setup
rm(myTTR.data, model.data, train.data, test.data)
myTTR.data <- myTTR(myStock)[,c("myATR", "myADX", "myMACD", "mySAR", "mySMI")]
model.data <- specifyModel(formula=T.index(myStock, p=0.025, k=10) ~ myTTR.data + runMean(Cl(myStock)) )
train.data <- as.data.frame(modelData(model.data,
data.window=c(start(myStock), as.Date("2013-12-31")) )) #Books Error was here
test.data <- as.data.frame(modelData(model.data,
data.window=c(as.Date("2014-01-01"), end(myStock)) ))
colnames(train.data) <- c("T", "myATR", "myADX", "myMACD", "mySAR", "mySMI", "runMean")
colnames(test.data) <- c("T", "myATR", "myADX", "myMACD", "mySAR", "mySMI", "runMean")
form <- as.formula("T~.")
# Fitting into the model
svm.model <- svm(form, train.data, cost=100)
svm.predict <- predict(svm.model, na.omit(test.data))
# Transform from T.Index to Buy/Hold/Sell Signal
T2Signal <- function(x, a1=-0.01, a2=-a1){
result <- ifelse(x<a1, "Sell", ifelse(x>a2, "Buy", "Hold"))
result <- factor(result, levels=c("Buy", "Hold", "Sell"))
return(result)
}
accuracy <- function(prediction, true){
t <- table(prediction, true)
result <- (t["Sell", "Sell"] + t["Buy", "Buy"]) / (sum(t["Buy",]) + sum(t["Sell",]))
return(result)
}
accuracy2 <- function(prediction, true){
t <- table(prediction, true)
result <- (t["Sell", "Sell"] + t["Buy", "Buy"]) / (t["Sell", "Sell"] + t["Buy", "Buy"] + t["Sell", "Buy"] + t["Buy", "Sell"])
return(result)
}
# Result & Accuracy
signal.pred <- T2Signal(x=svm.predict, a1=-0.095, a2=0.095)
signal.true <- T2Signal(x=na.omit(test.data)$T, a1=-0.095, a2=0.095)
table(signal.pred, signal.true)
accuracy(signal.pred, signal.true)
accuracy2(signal.pred, signal.true)
# For Parameter Testing
a1 <- -seq(0.01, 0.2, length.out=50)
a2 <- -a1
ac <- function(a1,a2){
signal.pred <- T2Signal(x=svm.predict, a1=a1, a2=a2)
signal.true <- T2Signal(x=na.omit(test.data)$T, a1=a1, a2=a2)
accuracy(signal.pred, signal.true)
}
result <- outer(X=a1, Y=a2, FUN=Vectorize(ac, vectorize.args=c("a1", "a2")))
ind <- which(max(result)==result, arr.ind=TRUE)
cbind(a1[ind[,1]], a2[ind[,2]])
filled.contour(x=-a1*100, y=a2*100,z=result, color=terrain.colors,
plot.title=title(main="The Accuracy",
xlab=expression(paste(alpha[1], " ",10^-2)),
ylab=expression(paste(alpha[2], " ",10^-2))
),
key.title=title(main="Accuracy"),
key.axes=axis(4,seq(0,1,by=0.05)))
####### Simple version after function loaded ##########
myStock <- getSymbols(Symbols = "^HSI", from="2000-01-01", to="2016-9-30", src="yahoo", auto.assign=FALSE)
colnames(myStock) = paste("myStock", c("Open","High","Low","Close","Volume","Adjusted"), sep=".")
myTTR.data <- myTTR(myStock)[,c("myATR", "myADX", "myMACD", "mySAR", "mySMI")]
model.data <- specifyModel(formula=T.index(myStock, p=0.025, k=10) ~ myTTR.data + runMean(Cl(myStock)) )
train.data <- as.data.frame(modelData(model.data, data.window=c(start(myStock), as.Date("2013-12-31"))))
test.data <- as.data.frame(modelData(model.data, data.window=c(as.Date("2014-01-01"), end(myStock))))
colnames(train.data) <- c("T", "myATR", "myADX", "myMACD", "mySAR", "mySMI", "runMean")
colnames(test.data) <- c("T", "myATR", "myADX", "myMACD", "mySAR", "mySMI", "runMean")
form <- as.formula("T~.")
rm(svm.model, svm.predict, signal.pred, signal.true)
svm.model <- svm(form, train.data, cost=100)
svm.predict <- predict(svm.model, na.omit(test.data))
signal.pred <- T2Signal(x=svm.predict, a1=-0.095, a2=0.095)
signal.true <- T2Signal(x=na.omit(test.data)$T, a1=-0.095, a2=0.095)
table(signal.pred, signal.true)
accuracy(signal.pred, signal.true)
accuracy2(signal.pred, signal.true)
###########################################
沒有留言:
張貼留言