利用R语言编写量化投资策略-CDA数据分析师官网

热线电话：13121318867

利用R语言编写量化投资策略

2018-05-18

利用R语言编写量化投资策略

选取一股票，利用 R语言进行分析，同时构建通道突破，双均线交叉和MACD策略，进行回测。
library(xts)
library(xtsExtra)
library(quantmod)

library(FinTS)
library(forecast)
library(TSA)
library(TTR)

library(fGarch)
library(rugarch)
library(tseries)

setSymbolLookup(MHXX=list(name='0696.hk',src='yahoo'))
getSymbols("MHXX",from="2013-01-01",to="2015-09-30")
#显示K线图，如图明显发现股价呈现递增趋势，价格序列是非平稳的。
chartSeries(MHXX)

#考虑对数收益率
#获取收盘价
cp = MHXX[,4]
lgcp=log(MHXX[,4])

#tdx =c(1:456)/365+2014

#计算日收益率
ret=dailyReturn(MHXX)
chartSeries(ret,theme="white",TA=NULL)
#plot(tdx,cp,xlab="year",ylab="close price",type='l')

#计算对数收益率，如图课件，股价在15年左右有一个跳跃，15年第二季度的股价增长导致
#之后股价有较大的下降，这些特征给后续的分析带来一些较大的异常值

lgret = log(ret+1)
chartSeries(lgret,theme="white",TA=NULL)

#由ACF和PACF图可以看出,该股1股价的日收益率序列即使存在某种相关性，该自相关性也
#很小
par(mfcol=c(2,1))
acf(lgret,lag=30)
pacf(lgret,lag=30)
#为了验证该收益率序列有没有序列相关性，使用Ljung-Box检验，结果对应的P值0.024，
#在1%的显著水平下，拒绝该股票日收益率没有显著前后相关性的这一原假设。
#但在5%的显著水平下，无法拒绝该股票日收益率没有显著前后相关性的这一原假设。
Box.test(lgret,lag=20,type='Ljung')

##############################################################################
m1 <- auto.arima(lgret,stationary=TRUE,seasonal=FALSE,ic="aic")

#鉴于该股票对数收益率序列的自相关性并不强，所以建立的ARIMA模型可能适用性不高。
#对于对数收益率序列，单样本的t检验结果的t比为1.0625，p值为0.2884，表明该序列不是
#显著异于零的，同时此处根据ACF图所示，在4阶有轻微的超越标准差线，
#因此取用AR(5)模型拟合，aic=-2987.43
m2 <- arima(x=lgret,order=c(4,0,0),include.mean=F)
tratio=m2$coef/sqrt(diag(m2$var.coef))
tratio
meacf=eacf(lgret,6,12)
print(meacf$eacf,digits=2)
#残差检验并表示改模型可能不是充分的
tsdiag(m2,gof=20)

m3 <-auto.arima(ret,stationary = TRUE,seasonal = FALSE,ic="aic")
m3
################################################################################

#由上述可知，对于价格变化的分析，纯ARMA模型是不充分的，一方面ARMA模型不能处理
#波动率聚集，另一方面，ARMA-GARCH模型能充分处理这些数据的复杂性，
#并能提高样本外预测
price=ts(cp)
dp=ts(diff(cp))

par(mfcol=c(2,1))
plot(price,xlab='year',ylab='price')
plot(dp,xlab='year',ylab='changes')

cprice=diff(price)
par(mfcol=c(2,1))
acf(cprice)
pacf(cprice)

#aic=-0.37
m.garch1<-garchFit(~1+garch(1,1),data=cprice,trace=F)
summary(m.garch1)

#aic=-0.62
m.garch2<-garchFit(~arma(6,0)+garch(1,1),data=cprice,trace=F,ininclude.mean = F,
cond.dist = "std")
summary(m.garch2)

#aic=-0.60
m.garch3<-garchFit(~arma(2,0)+garch(1,1),data=cprice,trace=F,ininclude.mean = F,
cond.dist = "std")
summary(m.garch3)

#aic=-0.596
m.garch4<-garchFit(~arma(1,0)+garch(1,1),data=cprice,trace=F,ininclude.mean = F,
cond.dist = "std")
summary(m.garch4)

#回测检验
source("backtestGarch.R")
M2F=backtestGarch(cprice,714,2,inc.mean=F,cdist="sstd")

source("backtest.R")
M2AF=backtest(m2,cprice,714,2,inc.mean=F)
#ArchTest(coredata(ret))
################################################################################

#计算VaR
mgarch1<-ugarchspec(variance.model=list(garchOrder=c(1,1)),
mean.model=list(armaOrder=c(0,0)))
mgarch1_fit<-ugarchfit(spec=mgarch1,data=cprice)
mgarch1_fit

mgarch1_roll<-ugarchroll(mgarch1,cprice,n.start=120,refit.every=1,
refit.window = "moving",solver="hybrid",
calculate.VaR = TRUE,VaR.alpha = 0.01,keep.coef = TRUE)
report(mgarch1_roll,type="VaR",VaR.alpha=0.01,conf.level=0.99)

#生成PLOT
cprice_var<-zoo(mgarch1_roll@forecast$VaR[,1])
index(cprice_var)<-as.yearmon(rownames(mgarch1_roll@forecast$VaR))

cprice_actual<-zoo(mgarch1_roll@forecast$VaR[,2])
index(cprice_var)<-as.yearmon(rownames(mgarch1_roll@forecast$VaR))

plot(cprice_actual,type="b",main="99% day Var backtesting",xlab="Date",
ylab="Return /VaR in percent")
lines(cprice_var,col="red")
legend("topright",inset=.05,c("MHXX return","VaR"),col=c("black","red"),lty=c(1,1))

mgarch1_fcst <- ugarchforecast(mgarch1_fit, n.ahead = 6)
mgarch1_fcst

ret.fcst <- - qnorm(0.95) * mgarch1_fcst @forecast$sigmaFor
ret.fcst

chartSeries(MHXX,name="中国民航信息",TA=NULL)
addBBands()
#addMACD()

################################量化投资策略####################################
###### 通道突破 ######

#通道突破函数==================================================================
bband.bk.sim <- function(stk.prc.xts, k=20, p=1.65, q=0.8){
#q是交易倍数，表示资金的q分用于交易

stk.prc <- coredata(stk.prc.xts) #把主要数据取出
Timeline <- index(stk.prc.xts)
End <- length(stk.prc.xts)

MA <- c( rep(0, k), 0)
std <- c( rep(0, k), 0)
u.bound <- c( rep(0, k), 0)
signal <- c( rep(0, k), 0) #交易信号
trd.state <- c( rep(0, k), 0) #记录买卖状态
share <- c( rep(0, k), 0) #记录持股份数

cash <- c( rep(1e4, k), 0) #现金部位
value <- c( rep(1e4, k), 0) #资产价值=股票市值+现金部位

# Sim ----

for( t in k:End ){

stk.prc.pre <- stk.prc[(t-k):t]
MA[t] <- mean( stk.prc.pre )
std[t] <- sd( stk.prc.pre )
u.bound[t] <- MA[t] + p * std[t] #布林带上界

signal[t] <- 0 #默认不交易
if( stk.prc[t] > u.bound[t] ) signal[t] = 1
#当股票价格超出布林上界时，buy
if( stk.prc[t-1] > MA[t-1] & stk.prc[t] <= MA[t] ) signal[t] = -1
if( stk.prc[t-1] < MA[t-1] & stk.prc[t] >= MA[t] ) signal[t] = -1
#卖的情况

trd.state[t] <- trd.state[t-1]
cash[t] <- cash[t-1]
share[t] <- share[t-1]
value[t] <- value[t-1]

#更新交易状态、持股数目、现金金额
if( trd.state[t-1] == 0 & signal[t] == 1 ){
trd.state[t] <- 1
share[t] <- ( q * cash[t-1] ) / stk.prc[t]
cash[t] <- cash[t-1] - share[t]*stk.prc[t]
}

if( trd.state[t-1] == 1 & signal[t] == -1 ){
trd.state[t] <- 0
share[t] <- 0
cash[t] <- cash[t-1] + share[t-1]*stk.prc[t]
}

value[t] <- cash[t] + share[t]*stk.prc[t]
}

res <- cbind(stk.prc, signal, trd.state, share, cash, value)
names(res) <- c("prc", "signal", "trd.state", "share", "cash", "value")

return(res)
}
#通道突破函数END================================================================

res <- bband.bk.sim(cp)
head(res)
tail(res)
plot(res[,6],type='l',col='darkred',lty=1,lwd=2)

## 通道(end)

############################### 均线系统策略 ###################################
## 双均线交叉策略
mov.avg.sim <- function(stk.prc.xts, k=50, n=7, p=1.05, q=1.10, m=0.8){

stk.prc <- coredata(stk.prc.xts)
Timeline <- index(stk.prc.xts)
End <- length(stk.prc)

MA.5 <- SMA(stk.prc, 5) #计算5日均线
MA.20 <- SMA(stk.prc, 20) #计算20日均线

signal <- c( rep(0, k), 0)
trd.state <- c( rep(0, k), 0)
share <- c( rep(0, k), 0)

cash <- c( rep(1e4, k), 0)
value <- c( rep(1e4, k), 0)

# Sim -----

for( t in k:End ){

signal[t] <- 0

if( sum(MA.5[(t-n):(t-1)] > MA.20[(t-n):(t-1)]) == n
& stk.prc[t-1]/MA.20[t-1] > p) signal[t] <- 1

if( MA.5[t-1] >= MA.20[t-1] & MA.5[t] <= MA.20[t]) signal[t] <- -1
if( stk.prc[t-1]/MA.20[t-1] > q ) signal[t] <- -1

trd.state[t] <- trd.state[t-1]
cash[t] <- cash[t-1]
share[t] <- share[t-1]
value[t] <- value[t-1]

if( trd.state[t-1] == 0 & signal[t] == 1 ){
trd.state[t] <- 1
share[t] <- ( m * cash[t-1] ) / stk.prc[t]
cash[t] <- cash[t-1] - share[t]*stk.prc[t]
}

if( trd.state[t-1] == 1 & signal[t] == -1 ){
trd.state[t] <- 0
share[t] <- 0
cash[t] <- cash[t-1] + share[t-1]*stk.prc[t]
}

value[t] <- cash[t] + share[t]*stk.prc[t]
}

res <- xts( cbind(stk.prc, MA.5, MA.20, signal, trd.state, share, cash, value),
order.by=Timeline)
names(res) <- c("prc", "MA.5", "MA.20","signal", "trd.state",
"share", "cash", "value")
head(res)

return(res)
}
#双均线交叉策略END==============================================================
res.mov <- mov.avg.sim(cp)
head(res.mov)
tail(res.mov)
plot(res.mov[,6],type='l',lty=1,lwd=2)

## MACD(begin)
MACD.sim <- function(stk.prc.xts, k=50, m=0.8){

stk.prc <- coredata(stk.prc.xts)
Timeline <- index(stk.prc.xts)
End <- length(stk.prc)

macd.line <- MACD(stk.prc, nFast=12, nSlow=26, nSig=9)[, 1]
signal.line <- MACD(stk.prc, nFast=12, nSlow=26, nSig=9)[, 2]

signal <- c( rep(0, k), 0)
trd.state <- c( rep(0, k), 0)
share <- c( rep(0, k), 0)

cash <- c( rep(1e4, k), 0)
value <- c( rep(1e4, k), 0)

# Sim -----

for( t in (k+1):End ){

signal[t] <- 0

if( macd.line[t-1] <= signal.line[t-1] & macd.line[t] > signal.line[t]) signal[t] <- 1

if( macd.line[t-1] >= signal.line[t-1] & macd.line[t] < signal.line[t]) signal[t] <- -1

trd.state[t] <- trd.state[t-1]
cash[t] <- cash[t-1]
share[t] <- share[t-1]
value[t] <- value[t-1]

if( trd.state[t-1] == 0 & signal[t] == 1 ){
trd.state[t] <- 1
share[t] <- ( m * cash[t-1] ) / stk.prc[t]
cash[t] <- cash[t-1] - share[t]*stk.prc[t]
}

if( trd.state[t-1] == 1 & signal[t] == -1 ){
trd.state[t] <- 0
share[t] <- 0
cash[t] <- cash[t-1] + share[t-1]*stk.prc[t]
}

value[t] <- cash[t] + share[t]*stk.prc[t]
}

res <- cbind(stk.prc, macd.line, signal.line,
signal, trd.state, share, cash, value)
names(res) <- c("prc", "MACD.line", "signal.line",
"signal", "trd.state", "share", "cash", "value")
head(res)

return(res)
}
#MACD策略END==============================================================
res.macd <- MACD.sim(cp)
head(res.macd)
tail(res.macd)

plot(res.macd[,8],type='l',lty=1,lwd=2)

#收益率
ret.macd<-diff(res.macd[,8])
plot(ret.macd,type='l',col='red',lty=1,lwd=2)

#总收益
ret.macd.sum<-sum(ret.macd)
ret.macd.sum.ratio<-ret.macd.sum/(res.macd[1,8])