利用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])







## MACD(end)

转载自http://blog.sina.com.cn/s/blog_6789c0a80102vugv.html

转载自http://blog.csdn.net/newusb/article/details/49445125

标签