R语言期末实验报告预测股票河北钢铁

SAS与数据分析方法实验课程

实验报告

班级姓名学号

数学一班杨燕2012102184

数学一班胡露霞2012102175

实验地点：81008

实验时间：2015年1月11日

成绩：______

实验成绩评定表

项目分值优秀(100>x≥

90)

良好(90>x≥

80)

中等(80>x

≥70)

及格(70>x≥

60)

不及格

（x＜60)参考标准参考标准参考标准参考标准参考标准

实验态度15

实验态度认

真，科学作风

严谨，严格保

证实验时间并

按指导书中规

定的进度开展

各项工作。

实验态度比

较认真，科学

作风良好，能

按期圆满完

成指导书规

定的任务。

实验态度尚

好，遵守组

织纪律，保

证实验时

间，能按期

完成各项工

作。

实验态度尚

可，能遵守组

织纪律，能按

期完成任务。

学习马虎，纪

律涣散，工作

作风不严谨，

不能保证实验

时间和进度。

实验准备20

认真预习实验

指导书，实验

步骤明确，实

验准备充分。

预习实验指

导书，实验步

骤较明确，实

验准备较充

分。

预习实验指

导书，实验

步骤基本明

确。

预习实验指

导书，实验主

要步骤基本

明确。

不预习实验指

导书，实验主

要步骤不明

确。

实验能力40

实验方法正

确，数据处理

及数据符合

专业要求，实

验数据真实、

完整、准确，

对相关问题有

很强的分析能

力和概括能

力。

实验方法正

确，实验数据

真实较完整，

个别数据不

合格，数据处

理及数据较

符合专业要

求，对相关问

题有较强的

分析能力和

概括能力。

实验方法正

确，实验数

据基本完

整，个别数

据不合格，

数据处理及

数据基本符

合专业要

求，对相关

问题有一定

的分析能力

和概括能

力。

实验方法基

本正确，实

验数据有较

多不合格，但

态度端正，

能基本按要

求操作，但

整体效果不

理想，综合

能力很差。

实验方法不正

确，不能完成

实验报告，缺

乏对相关问题

分析和概况的

基本能力。

报告撰写

质量及规范化25

实验报告内容

完整，数据处

理规范，作图

认真细致、表

格规范、清晰，

字迹工整，完

实验报告内

容较完整，数

据处理较规

范，作图较细

致、表格较规

范，字迹较工

实验报告内

容达到基本

要求，数据

处理基本规

范，作图较

细致、表格

有实验报告，

内容不完整。

有实验步骤，

不详细。实验

报告格式欠

规范，作图不

实验报告数

据记录不完

整，不作图，

不处理数据。

完全抄袭别

人报告等。

全符合规范化要求。整。较规范，字

迹基本工

整。

太规范等。

指导教师评定成绩：

指导教师签名：年月日

预测股票市场收益

问题描述与目标

股票交易市场绝对是数据挖掘技术的实战之地，事实上，数据挖掘技术在处理大量历史数据时要比人工检验更加优越，更加高效，另外，研究者声称，市场通过价格调整而发生快速变化，导致很难捕捉到持续盈利的投资方式，这就是所谓的有效市场假说，它表示在一个市场中如果价格能反映所有可获得信息，这种市场就成为有效市。但这一理论很快就被更加松弛的版本取代，因为有效市场并不一定永远是有效的，所以这也为交易机随时留有余地。

一般股票交易的目的就是通过买单卖单来维持一组证券资产。长期目标就是通过交易行为获得尽可能多的利润。本章我们的主要内容就限定在这个领域，基于既定的有价证券全要求和初始资本，我们会尽量扩大某一交易方式下未来的交易收益。我们的交易策略就是把数据挖掘的结果作为基本的决策信号，用于指导交易行为。这里的数据挖掘首先是基于历史报价构建模型，然后通过模型预测指数未来的变化趋势。我们的模型会被整合成一个交易系统，用来依据预测结果进行交易决策。另外，我们通过交易系统的表现作为评价准则，简单的说就是系统的交易行为所带来的获利和损失，以及其他一些投资者感兴趣的统计值。因此，通过数据挖掘过程所获得知识进行交易，这些交易获得的最终收益才是我们评价系统的标准，而不仅仅是模型预测的准确性

二、从网站上获取数据

获取数据的一种方法是使用Yahoo财经网站提供的免费服务。

数据的结构：交易的日期、开盘价、最高价、最低价、收盘价、交易量、调整的收盘价

为简单起见，用的是股票指数的数据

处理时间序列的包zoo,xts,后者是前者的拓展

表示的处理时间的类;POSIXct/POSIXIt,date

从CSV文件读数据

从网络读取数据

Quantmod包里的getSymbols()

从数据库读数据：包RODBC,RMySQL

获取数据如下：

library(quantmod)

library(tseries)

YT <- as.xts(get.hist.quote("000709.sz",start="2013-01-01",

quote=c("Open", "High", "Low", "Close","Volume","AdjClose"))) head(HB)

Date Open High Low Close Volume AdjClose

2001-01-01 7.95 8.02 7.95 8.00 0 1.57

2001-01-02 8.02 8.10 7.95 8.05 1579200 1.58

2001-01-03 8.06 8.12 8.04 8.10 2260700 1.59

2001-01-04 8.10 8.15 8.01 8.04 2144700 1.57

2001-01-05 8.02 8.06 7.98 8.01 3095100 1.57

2001-01-08 8.01 8.20 7.90 8.19 4344700 1.60

三、定义预测任务

1、预测什么

将学习的交易策略假设我们可以预测未来几天的市场变化趋势，如果这个预测在未来被验证时正确的，那么基于该预测下达的交易指令将是获利的。

事实上，我们需要预测的是在未来K天中价格总的动态变化，并不是预测某个特定时间的一个特定价格。

假设每天的平均价格可以有一下公式来近似：

其中,Ci、Hi、Li分别为第i天的收盘价、最高价和最低价。

设Vi代表未来k天的平均价格相对今天收盘价的百分比变化（通常称为算数收益）：Vi=

我们把动态变化绝对值超过目标收益p%的变化进行累加作为一个指标变化T：

指标变化T用来找出在k天内，日平均价格明显高于目标变化的那些日期的变化之和。大的正T值意味着有几天的日平均报价高于今天收盘价的p%，这种情况表明有潜在的机会发出买入指令，因为有良好的预期价格会上涨。另一方面，大的负T值标明价格可能下降，可能进行卖出行动。如果T值接近零，则可能是由于价格平稳波动或价格涨跌互现，正的变化和负的变化互相抵消。

下面的函数实现这个简单的指标：

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))

for(x in 1:n.days) r[,x] <- Next(Delt(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)) else x

}

为了更好的理解指标T的性质，绘图的代码如下：

candleChart(last(HB,'3 months'),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()

函数绘制股票价格的K线图。K线图用一个彩色框和竖直的柱条来代表每日的报价情况。柱条代表当天的最高，最低价格，而框代表开盘价和收盘价。框的颜色用来表示框的顶部所代表的价格（开盘价和收盘价），即在一天中价格是下降的是橘色，上升的是绿色。

本案例的方法假设在时刻t的正确交易行为是和我们对未来k天价格的变化预期相关的。总之，我忙呢想要预测时刻t的正确交易信号。对于历史价格数据，我们将是计算每天的T 值，并用给定的界限值，用上面给出的方法得到每一天的正确信号。

2、预测的变量是什么

由于买入/卖出决策是在每天交易结束后，所以这里将集中于收盘价的分析。初始的输入变量由多个基于收盘价的过去收益构成。可以用下列公式来计算h天（算术）的收益，或百分比收益

这里Ci是第i天的收盘价。

TTR的技术指标：

平均真实范围（Average true range，ATR）

随机动量指数 ( Stochastic Momentum Index ，SMI)

威尔斯-威尔德（Welles Wilder）定向运动指数（ADX）

Aroon指标

布林带（Bollinger Bands）指标

蔡金波动（Chaikin Volatility）指标

收盘价位置指标（Close Location Value,CLV）

阿姆氏简易波动指标（Ease of Movement Value,EMV）

MACD指标

资金流向指数（Money Flow Index,MFI）

抛物线止损反转和波动性指标

我们对这些TTR函数的输出做了后续处理以获得一个单一值。以下函数实现了这个过程：

ATR <- function(x) ATR(HLC(x))[,'atr']

SMI <- function(x) SMI(HLC(x))[,'SMI']

ADX <- function(x) ADX(HLC(x))[,'ADX']

Aroon <- function(x) aroon(x[,c('High','Low')])$oscillator

BB <- function(x) BBands(HLC(x))[,'pctB']

ChaikinVol <- function(x) Delt(chaikinVolatility(x[,c("High","Low")]))[,1]

CLV <- function(x) EMA(CLV(HLC(x)))[,1]

EMV <- function(x) EMV(x[,c('High','Low')],x[,'Volume'])[,2]

MACD <- function(x) MACD(Cl(x))[,2]

MFI<- function(x) MFI(x[,c("High","Low","Close")], x[,"Volume"])

SAR <- function(x) SAR(x[,c('High','Close')]) [,1]

Volat <- function(x) volatility(OHLC(x),calc="garman")[,1]

上面描述的变量给出了预测指标T的未来取值的初始预测变量集合。

为了避免出现有偏差的结果，最终的测试数据不用于变量选择过程。用训练集数据构建随机森林模型，代码如下；

YT=na.omit(HB)

library(randomForest)

colnames(HB) <- c("Open", "High", "Low", "Close","Volume","Adjusted") data.model <- specifyModel(T.ind(HB) ~ Delt(Cl(HB),k=1:10)

+ATR(HB)+ADX(HB)+ Aroon(HB)+BB(HB)+ ChaikinVol(HB)

+CMO(Cl(HB)) + EMA(Delt(Cl(HB)))+Volat(HB)

+RSI(Cl(HB))+SAR(HB)+ runMean(Cl(HB)))

library(randomForest)

set.seed(1234)

rf <- buildModel(data.model,method='randomForest',

training.per=c(start(HB),index(HB["2007-12-31"])),

ntree=50,importance=T)

varImpPlot(rf@fitted.model,type=1)

如果要使交易系统适应更新的报价信息，这种方法将十分便利。

然后，我们用给定的界限值来筛选用于建模的变量子集，我们使用界限值4：imp <- importance(rf@fitted.model,type=1)

rownames(imp)[which(imp > 10)]

[1] "ADX.HB.ADX" "Volat.HB"

我们将在我们的模拟尝试中使用这些变量。利用变量重要性信息，我们可以得到最终用于建立模型的数据集：

colnames(HB) <- c("Open", "High", "Low", "Close","Volume","Adjusted")

data.model <- specifyModel(T.ind(HB) ~ Delt(Cl(HB),k=1)

+ATR(HB)+ADX(HB)+ Aroon(HB)+BB(HB)+ ChaikinVol(HB)

+CMO(Cl(HB)) + EMA(Delt(Cl(HB)))+Volat(HB)

+RSI(Cl(HB))+SAR(HB)+ runMean(Cl(HB)))

3、预测任务

用解释变量来预测T（回归问题），然后计算信号signal

Signal=

下面的代码构造本节下面部分两个预测任务的预测模型所应用的数据结构。

Tdata.train <- as.data.frame(modelData(data.model,

data.window=c('2001-01-01','2007-12-31'))) Tdata.eval <- na.omit(as.data.frame(modelData(data.model,

data.window=c('2007-01-01','2015-01-10')))) Tform <- as.formula('T.ind.HB~ .')

调用函数na.train（）是必要的，它可以避免在数据结束部分由于没有未来数据计算T值所导致的NA值。

4、模型评价准则

可以通过测量错误率来衡量信号预测，错误率的定义为

error.rate=,)

这里是第i个测试数据的模型预测值而真实的类标签为，是0/1损失函数，其定义为：,)=

有时候经常采用上面损失函数的对立测度，即精确度，定义为

1-error.rate。

通过表1，本案例的决策精确度和回溯精确度指标正式定义如下：

Prec=

表1 预测交易信号的分类矩阵

预测结果

买入持有卖出

卖出

真实结果持有买入买入

也可以单独对某个特定的信号（例如买信号或者卖信号）计算其独立的决策精确度和回溯精确度，例如：

经常把决策精确度和回溯精确度合并为一个统计量，称为F度量，合并后的统计量为：F=

这里0，它控制回溯精确度相对决策精确度的相对重要性。

四、预测模型

1、如何应用训练集数据来建模

我们考虑每一个模型应用会用到以下三个不同的方法：1）所有的测试时段都使用一个模型；2）每个w天更新数据增长窗口；3)用每隔w天的数据滑动窗口。

2、建模工具

①、人工神经网络

通过下面的公式，可以很容易地对数据集的每一列应用下列公式进行转换：，其中是原始变量X的均值，是变量X的标准偏差。

举例说明：

set.seed(1234)

library(LPCM)

library(DMwR)

library(nnet)

norm.data <- scale(Tdata.train)

nn <-nnet(Tform,norm.data[1:1000,],size=10,decay=0.1,maxit=1000,linout=T,trace=F) norm.preds <- predict(nn,norm.data[1001:1500,])

preds <- unscale(norm.preds,norm.data)

我们可以将数值型预测值转换成信号，然后使用统计方法。

sigs.nn <- trading.signals(preds,0.1,-0.1)

true.sigs <- trading.signals(Tdata.train[1001:1500,'T.ind.HB'],0.1,-0.1)

sigs.PR(sigs.nn,true.sigs)

结果：

precision recall

s 0.3521127 0.5681818

b 0.2883436 0.2685714

s+b 0.3244681 0.3973941

使用nnet()函数来完成这个任务与使用该函数解决回归问题很相似。使用训练数据并应用下面的代码进行演示：

set.seed(1234)

library(nnet)

signals <- trading.signals(Tdata.train[,'T.ind.HB'],0.1,-0.1)

norm.data <- data.frame(signals=signals,scale(Tdata.train[,-1]))

nn <- nnet(signals ~ .,norm.data[1:1000,],size=10,decay=0.1,maxit=1000,trace=F)

preds <- predict(nn,norm.data[1001:1500,],type='class')

②、支持向量机

在支持向量回归中，这个过程很相似，主要区别在于误差和相关损失的计算。这通过相助于所谓的不敏感损失函数，该函数形式如下：

下面，我们将提供使用R中这类模型的简单例子。我们从使用添加包e1071中的函数进行回归任务开始，代码如下：

library(e1071)

sv <- svm(Tform,Tdata.train[1:1000,],gamma=0.001,cost=10)

s.preds <- predict(sv,Tdata.train[1001:1500,])

sigs.svm <- trading.signals(s.preds,0.1,-0.1)

true.sigs <- trading.signals(Tdata.train[1001:1500,'T.ind.HB'],0.1,-0.1)

sigs.PR(sigs.svm,true.sigs)

结果：

precision recall

s 0.50000 0.007575758

b 0.60000 0.102857143

s+b 0.59375 0.061889251

我们可以观察到，支持向量机模型的决策精确度值比人工神经网络好的多，尽管回溯精确度值较低，下一步，我们考虑分类任务：

library(kernlab)

data <- cbind(signals=signals,Tdata.train[,-1])

ksv <- ksvm(signals ~ .,data[1:1000,],C=10)

ks.preds <- predict(ksv,data[10-1:1500,])

sigs.PR(ks.preds,data[1001:1500,1])

结果：

precision recall

s 0.5238095 0.25000000

b 0.1415929 0.09142857

s+b 0.2784091 0.15960912

③、多元自适应回归样条

多元自适应回归样条是自适应回归模型的一个例子。一个多元自适应回归样条模型具有以下一般形式：

mars(x)=

其中是常数，是基函数。

下面应用函数earth()进行回归的代码：

library(DMwR)

library(earth)

e <- earth(Tform,Tdata.train[1:1000,])

e.preds <- predict(e,Tdata.train[1001:1500,])

sigs.e <- trading.signals(e.preds,0.1,-0.1)

true.sigs <- trading.signals(Tdata.train[1001:1500,'T.ind.HB'],0.1,-0.1)

sigs.PR(sigs.e,true.sigs)

结果：

precision recall

s 0.2771084 0.1742424

b 0.4368421 0.4742857

s+b 0.3882784 0.3452769

五、从预测到实践

1、如何应用预测模型

从客观的角度看，我们的交易系统将可以采取两种交易头寸：多头头寸（也称为多头仓位）和空头头寸（也称为空头仓位）。简略的说，当认为价格会下降时开空头仓位，认为价格会上涨时开多头仓位。下面将描述我们的模型实验中将应用并进行比较的一些具有可信性的策略。我们只考虑着两个主要的交易策略，这两个策略所应用的参数可以有细微的变化（例如，持有期、预期收益率，或在每个仓位上投资的资金量）。

2、与交易相关的评价准则

我们将使用以下3个指标来衡量交易的整体结果：1）初始资本与测试期期末资本之间的净差额；2）净差额所代表的百分比收益率；3）买入并持有策略的超额回报。这个策略包括在测试期开始多头仓位和等待到最后平仓。用购买并持有的收益来衡量我们的交易策略和这个简单策略之间的差异。

对于与风险相关的测量，我们将使用夏普比率系数测量每单位风险的回报，风险由收益的标准偏差来衡量。我们也将算出跌幅最大值，以衡量出模型的最大连续累积损失。

最后，根据测试期所持有仓位的数、每个仓位的平均收益、盈利仓位的百分比以及其他相关性不大的绩效指标来评估它们的效果。

3、模型集成：仿真交易

在给出这种类型的仿真实例之前，我们需要提供用于模拟器的交易策略函数的更多细节。应该使用某种协议来编写这些函数，也就是说，它们应该意识到模拟器将如何调用它们，同时它们如何返回模拟器所期望的信息。

下面是用户定义的交易策略的实例：

policy.1 <- function(signals,market,opened.pos,money,

bet=0.2,hold.time=10,

exp.prof=0.025, max.loss= 0.05

)

{

d <- NROW(market) # this is th

e ID o

f today

orders <- NULL

nOs <- NROW(opened.pos)

# nothing to do!

if (!nOs && signals[d] == 'h') return(orders)

# First lets check if we can open new positions

# i) long positions

if (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+exp.prof),

market[d,'Close']*(1-max.loss) ),

action = c('open','close','close'),

posID = c(NA,NA,NA) ) )

# ii) short positions

} else if (signals[d] == 's' && !nOs) {

# this is the nr of stocks we already need to buy

# because of currently opened short positions

need2buy <- 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-exp.prof),

market[d,'Close']*(1+max.loss) ),

action = c('open','close','close'),

posID = c(NA,NA,NA) ) ) }

# Now lets check if we need to close positions

# because their holding time is over

if (nOs)

for(i in 1: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,exp.prof=0.025, max.loss= 0.05）

{

d <- NROW(market) # this is th

e ID o

f today

orders <- NULL

nOs <- NROW(opened.pos)

# nothing to do!

if (!nOs && signals[d] == 'h') return(orders)

# First lets check if we can open new positions

# i) long positions

if (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+exp.prof),

market[d,'Close']*(1-max.loss) ),

action = c('open','close','close'),

posID = c(NA,NA,NA)）)

# ii) short positions

} else if (signals[d] == 's') {

# this is the money already committed to buy stocks

# because of currently opened short positions

need2buy <- 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-exp.prof),

market[d,'Close']*(1+max.loss) ),

action = c('open','close','close'),

posID = c(NA,NA,NA)))}

}

定义了交易策略函数之后，我们就准备好了来尝试交易模拟器。为了演示方便，我们将选择数据集中的一个较小数据子集，建立支持向量机模型，然后用该模型来获取之后一个时间段的预测值。在某个交易策略下，用预测值来调用交易模拟器，得到利用支持向量机的交易信号所获得的交易结果。代码如下：

# Train and test periods

start <- 1

len.tr <- 1000

len.ts <- 500

tr <- start:(start+len.tr-1)

ts <- (start+len.tr):(start+len.tr+len.ts-1)

# getting the quotes for the testing period

data(HB)

date <- rownames(Tdata.train[start+len.tr,])

market <- HB[paste(date,'/',sep='')][1:len.ts]

# learning the model and obtaining its signal predictions

library(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)

# now using the simulated trader

t1 <- trading.simulator(market,sig,'policy.1',list(exp.prof=0.05,bet=0.2,hold.time=30))

在调用交易模拟器时，我们采用第一种交易策略，提供了一些不同的值给它的某些参数。我们使用默认的交易成本、默认的初始成本。指令的结果是一个tradeRecord类对象。我们可以检查该返回对象的内容，如下所示：

t1

Object of class tradeRecord with slots:

trading:

positions:

init.cap : 1e+06

trans.cost : 5

policy.func : policy.1

policy.pars :

summary(t1)

结果：

== Summary of a Trading Simulation with 500 days ==

Trading policy function : policy.1

Policy function parameters:

exp.prof = 0.05

bet = 0.2

hold.time = 30

Transaction costs : 5

Initial Equity : 1e+06

Final Equity : 973769.3 Return : -2.62 %

Number of trading positions: 49

Use function "tradingEvaluation()" for further stats on this simulation.

函数tradingEvaluation()用于获得在模拟交易期间的一系列表示交易效果的经济指标：

tradingEvaluation(t1)