【原创附代码】基于R软件的数据挖掘预测模型的股票交易系统

(3)(3)4基于数据挖掘预测模型的股票交易系统根据上市保险公司的业务分析和财务分析来看，对投资者来投资中国平安的价值最高，由于实验运行时间较长，数据挖掘分析的方法相同，与选择哪家上市保险公司历史交易数据无关。本文选择了中国人寿历史交易数据进行了数据挖掘与分析。4.1数据来源本文所用数据为中国人寿(601628)历史交易数据，数据来源于雅虎财经网站(网址：)。获取方法为如下R代码：library(tseries)CLI_Web_1<-as.xts(get.hist.quote("601628.ss",start="2007-01-09",quote=c("Open","High","Low","Close","Volume","AdjClose")))head(CLI_Web_1)并将所下载数据转换为R中的时间序列对象(xts对象)，本实验的数据开始时间为2007年1月9日，结束时间为2016年6月4日。4.2建模过程4.2.1数据处理用R函数colnames将下载数据整理成如下统一格式：OpenHighLowCloseVolumeAdjusted34.162007-01-0937.0040.2037.0038.9331901890034.162007-01-1039.8040.3038.7239.466861020034.632007-01-1138.8039.6037.0138.294390250033.602007-01-1237.7939.7437.5039.504217740034.662007-01-1539.8243.4538.9543.455613190038.132007-01-1645.2846.8844.0645.054456770039.534.2.2定义数据挖掘任务本模型所要解决的数据挖掘任务为预测任务。该预测的目标变量为Ti(R代码中名为T.ind)2007-01-1039.8040.3038.7239.466861020034.632007-01-1138.8039.6037.0138.294390250033.602007-01-1237.7939.7437.5039.504217740034.662007-01-1539.8243.4538.9543.455613190038.132007-01-1645.2846.8844.0645.054456770039.534.2.2定义数据挖掘任务本模型所要解决的数据挖掘任务为预测任务。该预测的目标变量为Ti(R代码中名为T.ind)，其定义如下:T=工{vgV:v>p%vvp%}i(1)其中P—Ciijij=1戸=C+H+LPiiii3(2)C、H、

iiL分别为交易日i的收盘价、最高价、最低价。V(i=1,...,k)为从交易ii日i开始接下来的1〜k天平均价格相对交易日i收盘价的变动率(称收益率)。Ti第1页共35页第第#页共35页间顺序。重复该过程R次也将保证训练集与测试集的充分变化，这增加了估计的可靠性。此外，如果使用同一个随机选取的时间点集R以评价不同的模型，则可以进行成对比较，并获得平均表现差值的统计置信水平。4.3.2实验结果的比较与模型选择结果本实验从2016年6月4日22:27产生第一个结果文件svmR.RData至2016年6月5日02:04产生最后一个(第5个)结果文件，耗时约2.5小时。本实验共产生360个学习模型，概述如下：>join(subset(svmR,stats=tgtStats),+subset(svmC,stats=tgtStats),+subset(nnetR,stats=tgtStats),+subset(nnetC,stats=tgtStats),+subset(earth,stats=tgtStats),+by='variants')==MonteCarloExperiment==20repetitionsMonteCarloSimulationusing:seed=1234trainsize=0.5xNROW(DataSet)testsize=0.4xNROW(DataSet)360learningsystemstestedon1datasets分别按准确率(prec)、收益率(Ret)、可获益交易的百分比(PercProf)、累计最大亏损(MaxDD)、Sharpe比率(SharpeRatio)(后二者与风险有关)各指标分别为学习模型排序，表现最佳的前五个模型分别为(R代码见附录)：>rankSystems(allSysRes,5,maxs=c(T,T,T,F,T))$SP500$SP500$prec.sbsystemscoregrow.svmC.v50.3387977grow.svmC.v130.3387977grow.svmC.v210.3387977grow.svmC.v60.3284392grow.svmC.v140.3284392$SP500$Retsystemscoregrow.svmR.v1975.5605grow.svmR.v2064.8880grow.earth.v4432.2415single.nnetR.v1230.0270single.earth.v2229.2260$SP500$PercProfsystemscoregrow.svmR.v466.1595grow.svmR.v366.0185grow.svmR.v165.9255grow.svmR.v265.4300slide.svmR.v565.1710$SP500$MaxDDsystemscoregrow.svmR.v750387.93grow.svmR.v453116.55grow.svmR.v154153.04slide.svmR.v554761.33grow.svmR.v354830.22$SP500$SharpeRatiosystemscoregrow.svmR.v190.0560grow.svmR.v200.0475grow.earth.v440.0315grow.svmR.v40.0310slide.svmR.v170.0305从中可见，支持向量机类模型svrR表现最优，时间窗口类型为增加型的模型表现表现突出。如果按同时满足(1)合理的交易次数至少为50次(2)平均收益率至少为0.7%(3)获益的交易百分比至少为60%三个条件选择模型，则结果如下：>intersect(intersect(s1,s2),s3)[1]"single.svmR.v2""slide.svmR.v2""slide.svmR.v5""slide.svmR.v6"[5]"grow.svmR.v1""grow.svmR.v2""grow.svmR.v3""grow.svmR.v4"可见，这些模型都是支持向量机模型，但滑动时间窗口与增长时间窗口几乎各占一半。根据模型建议交易决策20个例子：>head(t1@trading,20)CloseOrderMoneyN.StocksEquity2009-03-2522.1301000000.001000000.02009-03-2622.7201000000.001000000.02009-03-2722.4701000000.001000000.02009-03-3022.0601000000.001000000.0

2009-03-3122.9901000000.001000000.02009-04-0123.1101000000.001000000.02009-04-0223.7601000000.001000000.02009-04-0323.5101000000.001000000.02009-04-0723.7401000000.001000000.02009-04-0822.6401000000.001000000.02009-04-0923.59-11200173.4-8834991779.42009-04-1024.401990166.60990166.62009-04-1325.34-11193305.1-8116987645.62009-04-1425.391985368.20985368.22009-04-1524.90-11180965.6-7762987691.82009-04-1624.6301180965.6-7762989787.52009-04-1724.4401180965.6-7762991262.32009-04-2024.761993737.20993737.22009-04-2124.09-11186380.2-8027993009.82009-04-2223.451997564.10997564.1其中1表示买入指令，-1卖出指令，0为持有指令。附录本文用到的R代码：#####################################################ReadingthedatafromtheCSVfile###################################################library(DMwR)library(xts)CLI<-as.xts(read.zoo('ChinaLifeInsurance.csv',header=T,sep=","))######################################################GettingthedatafromtheWeb#########################################################################################################Thefirstway:Usingget.hist.quoteintseries###################################################library(tseries)CLI_Web_1<-as.xts(get.hist.quote("601628.ss",start="2007-01-09",quote=c("Open","High","Low","Close","Volume","AdjClose")))head(CLI_Web_1)######################################################Thesecondway:UsinggetSymbolsinquantmod###################################################library(quantmod)getSymbols('601628.ss')CLI_Web_2<-get("601628.SS")head(CLI_Web_2)colnames(CLI_Web_2)<-c("Open","High","Low","Close","Volume","Adjusted")head(CLI_Web_2)######################################################Thesecondway:getseveraldatasets###################################################setSymbolLookup(CLI_Web_3=list(name='601628.ss',src='yahoo'),NCI=list(name='601336.ss',src='yahoo'),CPA=list(name='601318.ss',src='yahoo'),CPac=list(name='601601.ss',src='yahoo'))getSymbols(c('CLI_Web_3','NCI','CPA','CPac'))head(CLI_WEB_3)head(NCI)head(CPA)head(CPAC)colnames(CLI_WEB_3)<-c("Open","High","Low","Close","Volume","Adjusted")

colnames(NCI)<-c("Open","High","Low","Close","Volume","Adjusted")colnames(CPA)<-c("Open","High","Low","Close","Volume","Adjusted")colnames(CPAC)<-c("Open","High","Low","Close","Volume","Adjusted")head(CLI_WEB_3)head(NCI)head(CPA)head(CPAC)此前加载数据######################################################DefiningthePredictionTasks###################################################T.ind<-function(quotes,tgt.margin=0.025,n.days=10){v<-apply(HLC(quotes),1,mean)r<-matrix(NA,ncol=n.days,nrow=NROW(quotes))##Thefollowingstatmentiswronginthebook(page109)!for(xin1:n.days)r[,x]<-Next(Delt(Cl(quotes),v,k=x),x)x<-apply(r,1,function(x)sum(x[x>tgt.margin|x<-tgt.margin]))if(is.xts(quotes))xts(x,time(quotes))elsex}此前计算文中所述目标变量CLI<-GSPC1<-CLI_WEB_3xts_date<-index(GSPC1)[coredata(GSPC1)[,5]!=0]xts_data<-coredata(GSPC1)[coredata(GSPC1)[,5]!=0,]GSPC<-xts(xts_data,order.by=xts_date)candleChart(last(CLI,'3months'),theme='white',TA=NULL)avgPrice<-function(p)apply(HLC(p),1,mean)addAvgPrice<-newTA(FUN=avgPrice,col=1,legend='AvgPrice')addT.ind<-newTA(FUN=T.ind,col='red',legend='tgtRet')addAvgPrice(on=1)addT.ind()library(DMwR)library(xts)CLI<-as.xts(read.zoo('002024.csv',header=T,sep=","))colnames(CLI)<-c("Open","High","Low","Close","Volume","Adjusted")library(quantmod)T.ind<-function(quotes,tgt.margin=0.025,n.days=10){v<-apply(HLC(quotes),1,mean)r<-matrix(NA,ncol=n.days,nrow=NROW(quotes))##Thefollowingstatmentiswronginthebook(page109)!for(xin1:n.days)r[,x]<-Next(Delt(Cl(quotes),v,k=x),x)x<-apply(r,1,function(x)sum(x[x>tgt.margin|x<-tgt.margin]))if(is.xts(quotes))xts(x,time(quotes))elsex}GSPC<-GSPC1<-CLIcandleChart(last(GSPC1,'3months'),theme='white',TA=NULL)avgPrice<-function(p)apply(HLC(p),1,mean)addAvgPrice<-newTA(FUN=avgPrice,col=1,legend='AvgPrice')addT.ind<-newTA(FUN=T.ind,col='red',legend='tgtRet')addAvgPrice(on=1)addT.ind()【已调试代码】library(TTR)myATR<-function(x)ATR(HLC(x))[,'atr']mySMI<-function(x)SMI(HLC(x))[,'SMI']myADX<-function(x)ADX(HLC(x))[,'ADX']myAroon<-function(x)aroon(x[,c('High','Low')])$oscillatormyBB<-function(x)BBands(HLC(x))[,'pctB']myChaikinVol<-function(x)Delt(chaikinVolatility(x[,c("High","Low")]))[,1]myCLV<-function(x)EMA(CLV(HLC(x)))[,1]myEMV<-function(x)EMV(x[,c('High','Low')],x[,'Volume'])[,2]myMACD<-function(x)MACD(Cl(x))[,2]myMFI<-function(x)MFI(x[,c("High","Low","Close")],x[,"Volume"])mySAR<-function(x)SAR(x[,c('High','Close')])[,1]myVolat<-function(x)volatility(OHLC(x),calc="garman")[,1]此前，计算12个技术指标，作为模型解释变量library(randomForest)data.model<-specifyModel(T.ind(GSPC)~Delt(Cl(GSPC),k=1:10)+myATR(GSPC)+mySMI(GSPC)+myADX(GSPC)+myAroon(GSPC)+myBB(GSPC)+myChaikinVol(GSPC)+myCLV(GSPC)+CMO(Cl(GSPC))+EMA(Delt(Cl(GSPC)))+myEMV(GSPC)+myVolat(GSPC)+myMACD(GSPC)+myMFI(GSPC)+RSI(Cl(GSPC))+mySAR(GSPC)+runMean(Cl(GSPC))+runSD(Cl(GSPC)),na.rm=TRUE)set.seed(1234)rf<-buildModel(data.model,method='randomForest',training.per=c(start(GSPC),index(GSPC["2012-01-09"])),ntree=50,importance=T)varImpPlot(,type=1)imp<-importance(,type=1)rownames(imp)[which(imp>6.5)]data.model<-specifyModel(T.ind(GSPC)~#Delt(Cl(GSPC),k=1)#+Delt(Cl(GSPC),k=2)#+Delt(Cl(GSPC),k=3)#+Delt(Cl(GSPC),k=4)#+Delt(Cl(GSPC),k=5)#+Delt(Cl(GSPC),k=6)#+Delt(Cl(GSPC),k=7)#+Delt(Cl(GSPC),k=8)#+Delt(Cl(GSPC),k=9)#+Delt(Cl(GSPC),k=10)+myATR(GSPC)+mySMI(GSPC)+myADX(GSPC)#+myAroon(GSPC)#+myBB(GSPC)#+myChaikinVol(GSPC)+myCLV(GSPC)#+CMO(Cl(GSPC))#+EMA(Delt(Cl(GSPC)))#+myEMV(GSPC)+myVolat(GSPC)+myMACD(GSPC)#+myMFI(GSPC)#+RSI(Cl(GSPC))+mySAR(GSPC)+runMean(Cl(GSPC))#+runSD(Cl(GSPC)))Tdata.train<-as.data.frame(modelData(data.model,data.window=c('2007-01-09','2012-01-09')))Tdata.eval<-na.omit(as.data.frame(modelData(data.model,data.window=c('2012-01-10','2016-06-2'))))Tform<-as.formula('T.ind.GSPC~.')######################################################ThePredictionModels###################################################set.seed(1234)library(nnet)norm.data<-scale(Tdata.train)nn1<-nn<-nnet(Tform,norm.data[1:700,],size=10,decay=0.01,maxit=1000,linout=T,trace=F)norm.preds<-predict(nn,norm.data[701:1000,])preds<-unscale(norm.preds,norm.data)sigs.nn<-trading.signals(preds,0.1,-0.1)true.sigs<-trading.signals(Tdata.train[701:1000,'T.ind.GSPC'],0.1,-0.1)sigs.PR(sigs.nn,true.sigs)summary(nn1)set.seed(1234)library(nnet)signals<-trading.signals(Tdata.train[,'T.ind.GSPC'],0.1,-0.1)norm.data<-data.frame(signals=signals,scale(Tdata.train[,-1]))nn2<-nn<-nnet(signals~.,norm.data[1:500,],size=10,decay=0.01,maxit=1000,trace=F)preds<-predict(nn,norm.data[501:1000,],type='class')sigs.PR(preds,norm.data[501:1000,1])summary(nn2)library(e1071)sv<-svm(Tform,Tdata.train[1:700,],gamma=0.001,cost=100)s.preds<-predict(sv,Tdata.train[701:1000,])sigs.svm<-trading.signals(s.preds,0.1,-0.1)true.sigs<-trading.signals(Tdata.train[701:1000,'T.ind.GSPC'],0.1,-0.1)sigs.PR(sigs.svm,true.sigs)summary(sv)library(kernlab)data<-cbind(signals=signals,Tdata.train[,-1])ksv<-ksvm(signals~.,data[1:500,],C=10)ks.preds<-predict(ksv,data[501:1000,])sigs.PR(ks.preds,data[501:1000,1])summary(ksv)library(earth)e<-earth(Tform,Tdata.train[1:500,])e.preds<-predict(e,Tdata.train[501:1000,])sigs.e<-trading.signals(e.preds,0.1,-0.1)true.sigs<-trading.signals(Tdata.train[501:1000,'T.ind.GSPC'],0.1,-0.1)sigs.PR(sigs.e,true.sigs)summary(e)######################################################FromPredictionsintoActions###################################################policy.1<-function(signals,market,opened.pos,money,bet=0.2,hold.time=10,f=0.025,max.loss=0.05){d<-NROW(market)#thisistheIDoftodayorders<-NULLnOs<-NROW(opened.pos)#nothingtodo!if(!nOs&&signals[d]=='h')return(orders)#Firstletscheckifwecanopennewpositions#i)longpositionsif(signals[d]=='b'&&!nOs){quant<-round(bet*money/market[d,'Close'],0)if(quant>0)orders<-rbind(orders,data.frame(order=c(1,-1,-1),order.type=c(1,2,3),val=c(quant,market[d,'Close']*(1+f),market[d,'Close']*(1-max.loss)),action=c('open','close','close'),posID=c(NA,NA,NA)))#ii)shortpositions}elseif(signals[d]=='s'&&!nOs){#thisisthenrofstockswealreadyneedtobuy#becauseofcurrentlyopenedshortpositionsneed2buy<-sum(opened.pos[opened.pos[,'pos.type']==-1,"N.stocks"])*market[d,'Close']quant<-round(bet*(money-need2buy)/market[d,'Close'],0)if(quant>0)orders<-rbind(orders,data.frame(order=c(-1,1,1),order.type=c(1,2,3),val=c(quant,market[d,'Close']*(1-f),market[d,'Close']*(1+max.loss)),action=c('open','close','close'),posID=c(NA,NA,NA)))}#Nowletscheckifweneedtoclosepositions#becausetheirholdingtimeisoverif(nOs)for(iin1:nOs){if(d-opened.pos[i,'Odate']>=hold.time)orders<-rbind(orders,data.frame(order=-opened.pos[i,'pos.type'],order.type=1,val=NA,action='close',posID=rownames(opened.pos)[i]))orders}policy.2<-function(signals,market,opened.pos,money,bet=0.2,f=0.025,max.loss=0.05){d<-NROW(market)#thisistheIDoftodayorders<-NULLnOs<-NROW(opened.pos)#nothingtodo!if(!nOs&&signals[d]=='h')return(orders)#Firstletscheckifwecanopennewpositions#i)longpositionsif(signals[d]=='b'){quant<-round(bet*money/market[d,'Close'],0)if(quant>0)orders<-rbind(orders,data.frame(order=c(1,-1,-1),order.type=c(1,2,3),val=c(quant,market[d,'Close']*(1+f),market[d,'Close']*(1-max.loss)),action=c('open','close','close'),posID=c(NA,NA,NA)))#ii)shortpositions}elseif(signals[d]=='s'){#thisisthemoneyalreadycommittedtobuystocks#becauseofcurrentlyopenedshortpositionsneed2buy<-sum(opened.pos[opened.pos[,'pos.type']==-1,"N.stocks"])*market[d,'Close']quant<-round(bet*(money-need2buy)/market[d,'Close'],0)if(quant>0)orders<-rbind(orders,data.frame(order=c(-1,1,1),order.type=c(1,2,3),val=c(quant,market[d,'Close']*(1-f),market[d,'Close']*(1+max.loss)),action=c('open','close','close'),posID=c(NA,NA,NA)))}orders}#Trainandtestperiodsstart<-1len.tr<-500len.ts<-500tr<-start:(start+len.tr-1)ts<-(start+len.tr):(start+len.tr+len.ts-1)#gettingthequotesforthetestingperiod###data(GSPC)date<-rownames(Tdata.train[start+len.tr,])market<-GSPC[paste(date,'/',sep='')][1:len.ts]#learningthemodelandobtainingitssignalpredictionslibrary(e1071)s<-svm(Tform,Tdata.train[tr,],cost=10,gamma=0.01)p<-predict(s,Tdata.train[ts,])sig<-trading.signals(p,0.1,-0.1)#nowusingthesimulatedtradert1<-trading.simulator(market,sig,'policy.1',list(f=0.05,bet=0.2,hold.time=30))t1summary(t1)tradingEvaluation(t1)plot(t1,market,theme='white',name='CLI')t2<-trading.simulator(market,sig,'policy.2',list(f=0.05,bet=0.3))summary(t2)tradingEvaluation(t2)start<-100len.tr<-500len.ts<-500tr<-start:(start+len.tr-1)ts<-(start+len.tr):(start+len.tr+len.ts-1)s<-svm(Tform,Tdata.train[tr,],cost=10,gamma=0.01)p<-predict(s,Tdata.train[ts,])sig<-trading.signals(p,0.1,-0.1)t2<-trading.simulator(market,sig,'policy.2',list(f=0.05,bet=0.3))summary(t2)tradingEvaluation(t2)######################################################ModelEvaluationandSelection###################################################MC.svmR<-function(form,train,test,b.t=0.1,s.t=-0.1,...){require(e1071)t<-svm(form,train,...)p<-predict(t,test)trading.signals(p,b.t,s.t)}MC.svmC<-function(form,train,test,b.t=0.1,s.t=-0.1,...){require(e1071)tgtName<-all.vars(form)[1]train[,tgtName]<-trading.signals(train[,tgtName],b.t,s.t)t<-svm(form,train,...)p<-predict(t,test)factor(p,levels=c('s','h','b'))}MC.nnetR<-function(form,train,test,b.t=0.1,s.t=-0.1,...){require(nnet)t<-nnet(form,train,...)p<-predict(t,test)trading.signals(p,b.t,s.t)}MC.nnetC<-function(form,train,test,b.t=0.1,s.t=-0.1,...){require(nnet)tgtName<-all.vars(form)[1]train[,tgtName]<-trading.signals(train[,tgtName],b.t,s.t)t<-nnet(form,train,...)p<-predict(t,test,type='class')factor(p,levels=c('s','h','b'))}MC.earth<-function(form,train,test,b.t=0.1,s.t=-0.1,...){require(earth)t<-earth(form,train,...)p<-predict(t,test)trading.signals(p,b.t,s.t)}singleModel<-function(form,train,test,learner,policy.func,...){p<-do.call(paste('MC',learner,sep='.'),list(form,train,test,...))eval.stats(form,train,test,p,policy.func=policy.func)}slide<-function(form,train,test,learner,relearn.step,policy.func,...){real.learner<-learner(paste('MC',learner,sep='.'),pars=list(...))p<-slidingWindowTest(real.learner,form,train,test,relearn.step)p<-factor(p,levels=1:3,labels=c('s','h','b'))eval.stats(form,train,test,p,policy.func=policy.func)}grow<-function(form,train,test,learner,relearn.step,policy.func,...){real.learner<-learner(paste('MC',learner,sep='.'),pars=list(...))p<-growingWindowTest(real.learner,form,train,test,relearn.step)p<-factor(p,levels=1:3,labels=c('s','h','b'))eval.stats(form,train,test,p,policy.func=policy.func)}eval.stats<-function(form,train,test,preds,b.t=0.1,s.t=-0.1,...){#SignalsevaluationtgtName<-all.vars(form)[1]test[,tgtName]<-trading.signals(test[,tgtName],b.t,s.t)st<-sigs.PR(preds,test[,tgtName])dim(st)<-NULLnames(st)<-paste(rep(c('prec','rec'),each=3),c('s','b','sb'),sep='.')#Tradingevaluationdate<-rownames(test)[1]market<-GSPC[paste(date,"/",sep='')][1:length(preds),]trade.res<-trading.simulator(market,preds,...)c(st,tradingEvaluation(trade.res))}pol1<-function(signals,market,op,money)policy.1(signals,market,op,money,bet=0.2,f=0.025,max.loss=0.05,hold.time=10)pol2<-function(signals,market,op,money)policy.1(signals,market,op,money,bet=0.2,f=0.05,max.loss=0.05,hold.time=20)pol3<-function(signals,market,op,money)policy.2(signals,market,op,money,bet=0.5,f=0.05,max.loss=0.05)#ThelistoflearnerswewilluseTODO<-c('svmR','svmC','earth','nnetR','nnetC')#ThedatasetsusedinthecomparisonDSs<-list(dataset(Tform,Tdata.train,'SP500'))#MonteCarlo(MC)settingsusedMCsetts<-mcSettings(20,#20repetitionsoftheMCexps0.5,#~10yearsfortraining0.4,#~5yearsfortesting1234)#randomnumbergeneratorseed#VariantstotryforalllearnersVARS<-list()VARS$svmR<-list(cost=c(