• Archive by category "R语言实践"

Blog Archives

当R语言遇上Docker

R的极客理想系列文章,涵盖了R的思想,使用,工具,创新等的一系列要点,以我个人的学习和体验去诠释R的强大。

R语言作为统计学一门语言,一直在小众领域闪耀着光芒。直到大数据的爆发,R语言变成了一门炙手可热的数据分析的利器。随着越来越多的工程背景的人的加入,R语言的社区在迅速扩大成长。现在已不仅仅是统计领域,教育,银行,电商,互联网….都在使用R语言。

要成为有理想的极客,我们不能停留在语法上,要掌握牢固的数学,概率,统计知识,同时还要有创新精神,把R语言发挥到各个领域。让我们一起动起来吧,开始R的极客理想。

关于作者:

  • 张丹(Conan), 程序员R,Nodejs,Java
  • weibo:@Conan_Z
  • blog: http://blog.fens.me
  • email: bsspirit@gmail.com

转载请注明出处:
http://blog.fens.me/r-docker

r-docker

前言

R语言作为数据分析的工具,已经广泛被大家所接受并使用。但要把R语言项目工程化,部署到生产环境,提供在线用户使用却是难度很大的。主要原因就是R本身是单线程的,不支持并行处理。

当R遇到上了Docker会发生什么呢?本文将做详细的解释。

目录

  1. 当R遇到上了Docker
  2. 用Docker来管理R的程序

1. 当R遇到上了Docker

前言中提到,R运行时环境是单线程的,不支持并行处理,所以我们很难把R直接应用到生产环境中。当R遇到上了Docker,就出现了一个可以解决上面问题的方案。

通过Docker的容器化技术,把R的应用Docker化。每当用户发出请求,程序可以自动地在线启动一个Docker化的容器,来装载R的任务,部署,运行,计算,并返回结果。

r-docker2

从极端的情况考虑,如果要面对100万次并发的请求,我们需要启动100万个Docker的容器,每次容器单独执行自己的任务。但这种情况是要避免的,因为R本身来说,是做数据任务的,并不善于处理web是请求。如果可以把用户的大批量请求,转换成少量的数据计算的任务,那么这个设计就完美地解决了R由于并发而不能被工程化的问题。

r-docker3

比如,针对大量用户的重复性计算,把R的计算结果保存在缓存池中。

2. 用Docker来管理R的程序

设计方案定好,接下来就是就是动手实践了。

操作过程分成4步:

  • 1. 要有Docker的环境
  • 2. 找到第三方成熟的R语言的Docker镜像
  • 3. 把我们的R程序装进去
  • 4. 打包,运行,上传

1. Docker的环境。

安装Docker环境,就不在本文中介绍了,Docker环境的安装,请参考文章在Ubuntu中安装Docker

2. 找到第三方成熟的R语言的Docker镜像。

在docker hub中,搜索关键字 r, 共有535条结果。我们直接选用,排在第一位的r-base做为Docker容器的基础就行了。

docker-r

从仓库中,下载r-base镜像。


# 下载r-base镜像,大概300mb要下一会儿
~ sudo docker pull r-base
Using default tag: latest
latest: Pulling from library/r-base
9cd73496e13f: Pull complete 
f10af350cd29: Pull complete 
eea7b33eea97: Pull complete 
c91475e50472: Pull complete 
1e5e5f6785b4: Pull complete 
8c4091261ff6: Pull complete 
Digest: sha256:5f06e5a89cc64cbc513d02a8c650ea8bcbf0499795add57d18793069795c6f8d
Status: Downloaded newer image for r-base:latest

# 查看本地镜像列表
~ sudo docker images
REPOSITORY          TAG                 IMAGE ID            CREATED             SIZE
bsspirit/fensme     latest              8496b10e857a        2 hours ago         182.8 MB
ubuntu              latest              f8d79ba03c00        2 weeks ago         126.4 MB
r-base              latest              e2abe45e47d7        3 weeks ago         959.9 MB

3. 把我们的R程序装进去。

把R程序放进去之前,我们要先通过命令交互的方法,看看r-base容器中,是什么样子的。

运行r-base容器,会直接打开一个R的命令行窗口。


~ sudo docker run -ti --rm r-base

R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> 

我们通过执行R语言程序,了解一下Docker环境的信息。


# R程序启动路径
> getwd()
[1] "/"

# 当前路径中的目录
> dir()
 [1] "bin"   "boot"  "dev"   "etc"   "home"  "lib"   "lib64" "media" "mnt"  
[10] "opt"   "proc"  "root"  "run"   "sbin"  "srv"   "sys"   "tmp"   "usr"  
[19] "var"  

# 用户身份
> system('whoami')
root

# 系统信息
> sessionInfo()
R version 3.3.1 (2016-06-21)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Debian GNU/Linux stretch/sid

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base    

# R语言版本
> R.Version()
$platform
[1] "x86_64-pc-linux-gnu"

$arch
[1] "x86_64"

$os
[1] "linux-gnu"

$system
[1] "x86_64, linux-gnu"

$status
[1] ""

$major
[1] "3"

$minor
[1] "3.1"

$year
[1] "2016"

$month
[1] "06"

$day
[1] "21"

$`svn rev`
[1] "70800"

$language
[1] "R"

$version.string
[1] "R version 3.3.1 (2016-06-21)"

$nickname
[1] "Bug in Your Hair"

通过上面的几条命令,r-base容器的系统环境已经掌握。接下来,我们可以写一个R的算法,让这个程序在r-base的容器中运行。退出容器。

新建项目目录


~ mkdir ret && cd ret
~ pwd
/home/conan/ret

我们用a.r写一个计算万科(WANK)000002.SZ股票收益率的程序。数据从yahoo财经进行采集,R语言用于收益率计算,计算结果通过在控制台打印。

wanke

新建R语言算法文件,a.r。


~ vi a.r

install.packages(c('quantmod','PerformanceAnalytics'))
library(quantmod)
library(PerformanceAnalytics)
VANKE<-getSymbols("000002.SZ",auto.assign = FALSE, from = '2010-10-10')
close<-VANKE$'000002.SZ.Close'
ret<-CalculateReturns(close, method = "discrete")
cumret<-cumprod((ret+1)[-1])-1
VANKE_ret<-merge(close,ret,cumret)
names(VANKE_ret)<-c('close','ret','cumret')
print(tail(VANKE_ret))

我们先在本机中运行这段代码。


> 安装类库
> install.packages(c('quantmod','PerformanceAnalytics'))
> # 装载类库
> library(quantmod)
> library(PerformanceAnalytics)
> 
> # 获得VANKE每K线数据
> VANKE<-getSymbols("000002.SZ",auto.assign = FALSE, from = '2010-10-10')
>
> # 收盘价
> close<-VANKE$'000002.SZ.Close'
> 
> # 每日收益率 = (T日收盘价 - (T-1日收盘价))/T-1日收盘价
> ret<-CalculateReturns(close, method = "discrete")
> 
> # 每日累计收盘率 = (T日收益率+1)*(T+1日收益率+1)*...*(T+N日收益率+1)-1
> cumret<-cumprod((ret+1)[-1])-1
> 
> # 合并数据集
> VANKE_ret<-merge(close,ret,cumret)
> names(VANKE_ret)<-c('close','ret','cumret')
> 
> # 查看VANKE最近几日收益率
> print(tail(VANKE_ret))
           close          ret   cumret
2016-08-18 25.58 -0.010444874 1.893665
2016-08-19 24.59 -0.038702111 1.781674
2016-08-22 24.70  0.004473363 1.794118
2016-08-23 24.70  0.000000000 1.794118
2016-08-24 23.99 -0.028744939 1.713801
2016-08-25 23.54 -0.018757816 1.662896

接下来,编写Dockerfile通过加载外部文件的方法。


~ vi Dockerfile

FROM r-base
COPY . /usr/local/src/myscripts
WORKDIR /usr/local/src/myscripts
CMD ["Rscript", "a.r"]

4. 打包,运行,上传。

打包,生成Docker的镜像文件a.r。


~ sudo docker build -t a.r .
[sudo] password for conan: 
Sending build context to Docker daemon 3.072 kB
Step 1 : FROM r-base
 ---> e2abe45e47d7
Step 2 : COPY . /usr/local/src/myscripts
 ---> e6ef215d3683
Removing intermediate container aaabfdfe92ab
Step 3 : WORKDIR /usr/local/src/myscripts
 ---> Running in e3f2c65b947a
 ---> c667baee06bf
Removing intermediate container e3f2c65b947a
Step 4 : CMD Rscript a.r
 ---> Running in dc040bbdd3b9
 ---> 9a48d6dc02fe
Removing intermediate container dc040bbdd3b9
Successfully built 9a48d6dc02fe

启动r-base容器,运行a.r的脚本。


~  sudo docker run a.r

看着大段的日志从眼前飞过,计算出了万科的收益率的结果。

docker-r2

最后一步,不忘上传到docker hub,仓库地址为:https://hub.docker.com/r/bsspirit/ret/

上传镜像的操作命令:


~ sudo docker tag 9a48d6dc02fe bsspirit/ret
~ sudo docker push bsspirit/ret

如果你有Docker的环境,你可以直接用下面的命令,进行容器下载和运行。


~ sudo docker run bsspirit/ret

R和Docker的相遇,给R提供了并行计算施展的空间。Docker和R的相遇,也让Docker能够切入数据处理领域,有了更广阔的应用场景。感谢R和Docker给程序员的世界,带来了新的机会!!

转载请注明出处:
http://blog.fens.me/r-docker

打赏作者

R语言解读自回归模型

R的极客理想系列文章,涵盖了R的思想,使用,工具,创新等的一系列要点,以我个人的学习和体验去诠释R的强大。

R语言作为统计学一门语言,一直在小众领域闪耀着光芒。直到大数据的爆发,R语言变成了一门炙手可热的数据分析的利器。随着越来越多的工程背景的人的加入,R语言的社区在迅速扩大成长。现在已不仅仅是统计领域,教育,银行,电商,互联网….都在使用R语言。

要成为有理想的极客,我们不能停留在语法上,要掌握牢固的数学,概率,统计知识,同时还要有创新精神,把R语言发挥到各个领域。让我们一起动起来吧,开始R的极客理想。

关于作者:

  • 张丹(Conan), 程序员R,Nodejs,Java
  • weibo:@Conan_Z
  • blog: http://blog.fens.me
  • email: bsspirit@gmail.com

转载请注明出处:
http://blog.fens.me/r-ar

r-ar

前言

时间序列是金融分析中常用到的一种数据格式,自回归模型是分析时间序列数据的一种基本的方法。通过建立自回归模型,找到数据自身周期性的规律,从而帮助我们理解金融市场的发展变化。

在时间序列分析中,有一个常用的模型包括AR,MA,ARMA,ARIMA,ARCH,GARCH,他们的主要区别是适用条件不同,且是层层递进的,后面的一个模型解决了前一个模型的某个固有问题。本文以AR模型做为开始,将对时间序列分析体系,进行完整的介绍,并用R语言进行模型实现。

由于本文为非统计的专业文章,所以当出现与教课书不符的描述,请以教课书为准。本文力求用简化的语言,来介绍自回归模型的知识,同时配合R语言的实现。

目录

  1. 自回归模型介绍
  2. 用R语言构建自回归模型
  3. 模型识别ACF/PACF
  4. 模型预测

1. 自回归模型(AR)

自回归模型(Autoregressive model),简称AR模型,是统计上一种处理时间序列的方法,用来描述当前值与历史值之间的关系,用变量自身的历史时间数据对自身进行预测,自回归模型必须满足平稳性的要求。比如,时间序列数据集X 的历史各期数据从X1至Xt-1,假设它们为线性关系,可以对当期Xt的表现进行预测。X的当期值等于一个或数个落后期的线性组合,加常数项,加随机误差。

p阶自回归过程的公式定义:

ar-forumla

字段解释:

  • Xt是当期的X的表现
  • c是常数项
  • p是阶数,i为从1到p的值
  • φi是自相关系数
  • t为时间周期
  • εt是均值为0,标准差为δ 的随机误差,同时δ是独立于t的

对于一阶自回模型,用AR(1)来表示,简化后的公式为:

ar1-forumla

自回归是从线性回归分析中发展而来,只是把自变量x对因变量y的分析,变成自变量x对自身的分析。如果你需要了解线性回归的知识,请参考文章R语言解读一元线性回归模型

自回归模型的限制

自回归模型是用自身的数据来进行预测,但是这种方法受到一定的限制:

  • 必须具有平稳性,平稳性要求随机过程的随机特征不随时间变化。
  • 必须具有自相关性,如果自相关系数(φi)小于0.5,则不宜采用,否则预测结果极不准确。
  • 自回归只适用于预测与自身前期相关的现象,即受自身历史因素影响较大的现象。对于受其他因素影响的现象,不宜采用自回归,可以改用向量自回归模型。

平稳性时间序列的特点

平稳性要求产生时间序列Y的随机过程的随机特征不随时间变化,则称过程是平稳的;假如该随机过程的随机特征随时间变化,则称过程是非平稳的。

平稳性是由样本时间序列所得到的拟合曲线,在未来的一段期间内能顺着现有的形态能一直地延续下去;如果数据非平稳,则说明样本拟合曲线的形态不具有延续的特点,也就是说拟合出来的曲线将不符合当前曲线的形态。

  • 随机变量Yt的均值和方差均与时间t无关
  • 随机变量Yt和Ys的协方差只与时间差(步长)t-s有关
  • 对于平稳时间序列在数学上有比较丰富的处理手段,非平稳的时间序列通过差分等手段转化为平稳时间序列处理

2. 用R语言构建自回归模型

了解了自回归模型的定义,我们就可以用R语言来模拟一下自回归模型的构建和计算过程。

生成一个随机游走的数据集,满足平稳性的要求。


# 随机游走的数据集
> set.seed(0)
> x<-w<-rnorm(1000)       # 生成符合正态分布N(0,1)的数据
> for(t in 2:1000) x[t]<-x[t-1]+w[t]
> tsx<-ts(x)              # 生成ts时间序列的数据集

# 查看数据集
> head(tsx,15)
 [1] 1.2629543 0.9367209 2.2665202 3.5389495 3.9535909 2.4136409
 [7] 1.4850739 1.1903534 1.1845862 3.5892396 4.3528331 3.5538238
[13] 2.4061668 2.1167053 1.8174901

> plot(tsx)            # 生成可视化图形 
> a<-ar(tsx);a         # 进行自回归建模

Call:
ar(x = tsx)

Coefficients:
     1  
0.9879  

Order selected 1  sigma^2 estimated as  1.168

数据的如图展示:
01

自相关系数为0.9879 ,这是一个非常强的自相关性,所以上述的数列符合自相关的特性。

R语言中ar()函数提供了多种自相关系数的估计,包括"yule-walker", "burg", "ols", "mle", "yw",默认是用yule-walker方法,常用的方法还有最小二乘法(ols),极大似然法(mle)。

我们用最小二乘法,来进行参数估计。


> b<-ar(tsx,method = "ols");b

Call:
ar(x = tsx, method = "ols")

Coefficients:
     1  
0.9911  

Intercept: -0.017 (0.03149) 

Order selected 1  sigma^2 estimated as  0.9906

用最小二乘法的计算结果,则自相关系统数为0.9911,截距为-0.017。只有使用最小二乘法进行参数估计的时候,才会有截距。

我们用极大似然法,来进行参数估计。


> d<-ar(tsx,method = "mle");d

Call:
ar(x = tsx, method = "mle")

Coefficients:
     1  
0.9904  

Order selected 1  sigma^2 estimated as  0.9902

用极大似然法计算结果,则自相关系统数为0.9904。对于上面3种估计方法,自相关系数的值都是很接近的。

3. 模型识别ACF/PACF

在上面的例子中,我们默认是用一阶的自回归模型AR(1),进行程序实现的。在实际应用中,自回归模型AR时间序列的阶数P是未知的,必须根据实际数据来决定,就要对AR模型定阶数。常的方法就是利用自相关函数(ACF)和偏自相关函数(PACF)来确定自回归模型的阶数。在ACF/PACF不能确定的情况下,还需要用AIC(Aikaike info Criterion)、BIC(Bayesian information criterion)的信息准则函数来确定阶数。

自回归模型的确立过程,是通过确定阶数,参数估计,再次确定阶数的方法进行判断。自相关函数ACF,用来确定采用自回归模型是否合适。如果自相关函数具有拖尾性,则AR模型为合适模型。偏自相关函数PACF用来确定模型的阶数,如果从某个阶数之后,偏自相关函数的值都很接近0,则取相应的阶数作为模型阶数,偏自相关函数通过截尾性确定阶数。

1. 自相关函数ACF(autocorrelation function)

将一个有序的随机变量序列与其自身相比较,这就是自相关函数在统计学中的定义。每个不存在相位差的序列,都与其自身相似,即在此情况下,自相关函数值最大。如果序列中的组成部分相互之间存在相关性(不再是随机的),则由以下相关值方程所计算的值不再为零,这样的组成部分为自相关。

自相关函数反映了同一序列在不同时序的取值之间的相关程序。

ACF的公式为:

acf-forumla

字段解释

  • Pk,为ACF的标准误差
  • t,为数据集的长度
  • k,为滞后,取值从1到t-1,表示相距 k个时间间隔的序列值之间的相关性
  • Yt,为样本在t时期的值
  • Yt-k,为样本在t-k时期的值
  • μ,为样本的均值

所得的自相关值Pk的取值范围为[-1,1],1为最大正相关值,-1则为最大负相关值,0为不相关。

根据上面公式,我们可以手动计算出tsx数据集的ACF值


> u<-mean(tsx)  #均值
> v<-var(tsx)   #方差

> # 1阶滞后
> p1<-sum((x[1:length(tsx)-1]-u)*(x[2:length(tsx)]-u))/((length(tsx)-1)*v);p1
[1] 0.9878619
> # 2阶滞后
> p2<-sum((x[1:(length(tsx)-2)]-u)*(x[3:length(tsx)]-u))/((length(tsx)-1)*v);p2
[1] 0.9760271
> # 3阶滞后
> p3<-sum((x[1:(length(tsx)-3)]-u)*(x[4:length(tsx)]-u))/((length(tsx)-1)*v);p3
[1] 0.9635961

同时,我们可以用R语言中的acf()函数来计算,会打印前30个滞后的ACF值。


> acf(tsx)$acf
, , 1

           [,1]
 [1,] 1.0000000
 [2,] 0.9878619
 [3,] 0.9760271
 [4,] 0.9635961
 [5,] 0.9503371
 [6,] 0.9384022
 [7,] 0.9263075
 [8,] 0.9142540
 [9,] 0.9024862
[10,] 0.8914740
[11,] 0.8809663
[12,] 0.8711005
[13,] 0.8628609
[14,] 0.8544984
[15,] 0.8462270
[16,] 0.8384758
[17,] 0.8301834
[18,] 0.8229206
[19,] 0.8161523
[20,] 0.8081941
[21,] 0.8009467
[22,] 0.7942255
[23,] 0.7886249
[24,] 0.7838154
[25,] 0.7789733
[26,] 0.7749697
[27,] 0.7709313
[28,] 0.7662547
[29,] 0.7623381
[30,] 0.7604101
[31,] 0.7577333

比较前3个值的计算结果,与我们自己的计算结果是一样的,同时可以用R语言进行可视化输出。


> acf(tsx)

02

从上图中看出,数据的ACF为拖尾,存在很严重的自相关性。接下来,这时候我们用偏自相关函数确定一下AR的阶数。

2. 偏自相关函数(PACF)(partial autocorrelation function)

偏自相关函数是有自相关函数推到而来。对于一个平稳AR(p)模型,求出滞后k自相关系数p(k)时,实际上得到并不是x(t)与x(t-k)之间单纯的相关关系。因为x(t)同时还会受到中间k-1个随机变量x(t-1)、x(t-2)、……、x(t-k+1)的影响,而这k-1个随机变量又都和x(t-k)具有相关关系,所以自相关系数p(k)里实际掺杂了其他变量对x(t)与x(t-k)的影响。

为了能单纯测度x(t-k)对x(t)的影响,引进偏自相关系数的概念。对于平稳时间序列{x(t)},所谓滞后k偏自相关系数指在给定中间k-1个随机变量x(t-1)、x(t-2)、……、x(t-k+1)的条件下,或者说,在剔除了中间k-1个随机变量x(t-1)、x(t-2)、……、x(t-k+1)的干扰之后,x(t-k)对x(t)影响的相关程度。

简单来说,就是自相关系数ACF还包含了其他变量的影响,而偏自相关系数PACF是严格这两个变量之间的相关性。在ACF中存在着线性关系和非线性关系,偏自相关函数就是把线性关系从自动关系性中消除。当PACF近似于0,表明两个时间点之间的关系性是完全由线性关系所造成的。

通过R语言的pacf()函数来进行偏自相关函数计算。


> pacf(tsx)$acf
, , 1

              [,1]
 [1,]  0.987861891
 [2,]  0.006463542
 [3,] -0.030541593
 [4,] -0.041290415
 [5,]  0.047921168
 [6,] -0.009774246
 [7,] -0.006267004
 [8,]  0.002146693
 [9,]  0.028782423
[10,]  0.014785187
[11,]  0.019307564
[12,]  0.060879259
[13,] -0.007254278
[14,] -0.004139848
[15,]  0.015707900
[16,] -0.018615370
[17,]  0.037067452
[18,]  0.019322565
[19,] -0.048471479
[20,]  0.023388065
[21,]  0.027640953
[22,]  0.051177900
[23,]  0.028063875
[24,] -0.003957142
[25,]  0.034030631
[26,]  0.004270416
[27,] -0.029613088
[28,]  0.033715973
[29,]  0.092337583
[30,] -0.031264028

# 可视化输出 
> pacf(tsx)

03

从上面的这个结果分析,当滞后为1时AR模型显著,滞后为其他值是PACF的值接近于0不显著。所以,对于数据集tsx来说,数据满足AR(1)的自回归模型。对于上文中参数估计出的1阶自相关系数值是可以用的。

4. 模型预测

通过模型识别,我们已经确定了数据集tsx是符合AR(1)的建模条件的,同时我们也创建了AR(1)模型。接下来,就可以利用这个自回测的模型的进行预测,通过规律发现价值。在R语言中,我们可以用predict()函数,实现预测的计算。

使用AR(1)模型进行预测,并保留前5个预测点。


> predict(a,10,n.ahead=5L)
$pred
Time Series:
Start = 2 
End = 6 
Frequency = 1 
[1] 9.839680 9.681307 9.524855 9.370303 9.217627

$se
Time Series:
Start = 2 
End = 6 
Frequency = 1 
[1] 1.080826 1.519271 1.849506 2.122810 2.359189

上面结果中,变量$pred表示预测值,变量$se为误差。

我可以生成可视化的图,更直观的看到预测的结果。


# 生成50个预测值 
> tsp<-predict(a,n.ahead=50L)

# 把原数据画图 
> plot(tsx)

# 把预测值和误差画出来
> lines(tsp$pred,col='red')                
> lines(tsp$pred+tsp$se,col='blue')
> lines(tsp$pred-tsp$se,col='blue')

04

图中,黑色线为原始数据的,红色线为预测值,蓝色线为预测值的范围。这样我们就利用AR(1)模型,实现了对规律的预测计算。

上面关于预测和可视化的过程,我们是通过原生的predict()函数和plot()函数完成的。在R语言中,可以用forecast包来简化上面的操作过程,让代码更少,操作更便捷。


# 加载forecast包
> library('forecast')

# 生成模型AR(1) 
> a2 <- arima(tsx, order=c(1,0,0))
> tsp2<-forecast(a2, h=50)
> plot(tsp2)

05

查看forecast()计算后的预测结果。


> tsp2
     Point Forecast     Lo 80      Hi 80     Lo 95       Hi 95
1001      -15.71590 -16.99118 -14.440628 -17.66627 -13.7655369
1002      -15.60332 -17.39825 -13.808389 -18.34843 -12.8582092
1003      -15.49181 -17.67972 -13.303904 -18.83792 -12.1456966
1004      -15.38136 -17.89579 -12.866932 -19.22685 -11.5358726
1005      -15.27197 -18.06994 -12.474000 -19.55110 -10.9928432
1006      -15.16362 -18.21425 -12.112996 -19.82915 -10.4980922
1007      -15.05631 -18.33593 -11.776682 -20.07206 -10.0405541
1008      -14.95001 -18.43972 -11.460312 -20.28705  -9.6129750
1009      -14.84474 -18.52891 -11.160567 -20.47919  -9.2102846
1010      -14.74046 -18.60591 -10.875013 -20.65216  -8.8287673
1011      -14.63718 -18.67257 -10.601802 -20.80877  -8.4655994
1012      -14.53489 -18.73030 -10.339486 -20.95121  -8.1185723
1013      -14.43357 -18.78024 -10.086905 -21.08123  -7.7859174
1014      -14.33322 -18.82333  -9.843112 -21.20026  -7.4661903
1015      -14.23383 -18.86034  -9.607319 -21.30947  -7.1581923
1016      -14.13538 -18.89190  -9.378864 -21.40985  -6.8609139

通过forecast()函数,直接生成了Forecast值,80%概率的预测值范围,和95%概率的预测值范围。

在明白了整个自回归模型的设计思路、建模过程、检验条件、预测计算、可视化展示的完整操作后,我们就可以真正地把自回归模型用到实际的业务中。发现规律,发现价值!!

自回归模型只是开始,下一篇继续介绍移动平均模型(MA)的建模和使用过程。

转载请注明出处:
http://blog.fens.me/r-ar

打赏作者

R语言量化投资常用包总结

用IT技术玩金融系列文章,将介绍如何使用IT技术,处理金融大数据。在互联网混迹多年,已经熟练掌握一些IT技术。单纯地在互联网做开发,总觉得使劲的方式不对。要想靠技术养活自己,就要把技术变现。通过“跨界”可以寻找新的机会,创造技术的壁垒。

金融是离钱最近的市场,也是变现的好渠道!今天就开始踏上“用IT技术玩金融”之旅!

关于作者:

  • 张丹(Conan), 程序员Java,R,Nodejs
  • weibo:@Conan_Z
  • blog: http://blog.fens.me
  • email: bsspirit@gmail.com

转载请注明出处:
http://blog.fens.me/r-quant-packages/

quant-packages

前言

总是被很多的人问,为什么用R语言做量化投资,R、Python、Matlab比起来哪个更好?其实,答案很简单,你哪个用的熟就用哪个,工具是用来提升效率的,结果才是你会得到的。认准一门语言,坚持把它做好你就会成长。

每个领域,每种编程语言都用推动它前进的人,跟上牛人的脚步,你慢慢地也会变牛。

目录

  1. 为什么用R语言做量化投资?
  2. 常用量化投资工具包

1. 为什么用R语言做量化投资?

R做量化投资到底有哪些优势呢?最主要的一点,就是R语言有很多第三方包的支持。通常编程语言的设计,都是为了解决软件开发和程序实现的问题。但R语言在开始时,就被设计为主要解决数据的问题。量化投资就是对数据进行各种数据处理、数据分析,从而找到数据的规律。所以,有很多从事量化投资的人,把R语言用来构建量化交易的模型,进行回测,风险管理等,最后把研究成果开源并贡献给R语言的社区,为后面的人提供了非常大的帮助。

相比Python来说也有很多的第三方包的支持,这些第三方大部分提供是Web开发,数据爬虫,系统管理,数据库调用,数学计算等,这些都是属于通用的软件需求,而非某个行业的数据需求。当某个Python大神,开始关注量化投资领域,并用Python实现了一套量化的程序库,后面的人就会进入这个领域,只是沿着大神的路线走,等待下一个大神的出现。所以本质上,Python是面向程序设计的语言,而R是面向数据的语言。

R语言在量化投资领域,已经有很多年的积累,很多的算法已经成型。从投资研究到交易分析,再到风险管理,有着完整的体系结构。我们同样可以沿着前人走出来的路,快速学习,快速搭建出量化投资的系统来。对于有IT但背景缺乏金融知识的人来说,有很多的部分知识上手比较困难,同时看不太懂各种统计指标,对学习造成了很大的阻力。这其实是你深入到具体地某个行业后,都会面临的问题。行业知识和数学知识才是最难的,只有突破了,你才能打开认知新领域的方法。

R语言让我们更接近数据,同时提供了各种数学统计的工具,又有大量由第三方贡献的行业知识库,所以我会选择R语言,我会把R语言作为最好的工具,进行量化投资的分析。

2. 常用量化投资工具包

R语言在金融领域提供了很多的金融计算框架和工具,当你具备金融理论知识和市场经验,你可以利用这些第三方提供的技术框架来构建自己的金融模型。我们可以从CRAN上找到各种的金融项目,访问R的官方网站 (https://cran.r-project.org/),找到Task Views 菜单里的 Finance标签。

task

金融领域涉及范围是非常广的,包括银行业、保险业、信托业、证券业、租赁业等。金融行业都具有指标性、垄断性、高风险性、效益依赖性和高负债经营性的特点。量化投资是证券投资的一个很细分的专业领域,涉及到的金融工具包其实并不是太多。我们其实能把这些工具包研究好了,就可以方便地做量化的模型和交易了。

如果我们想用R构建自己的量化交易系统,你需要用到5方面的R语言工具包:数据管理、指标计算、回测交易、投资组合、风险管理。

quant-lib

  • 数据管理:包括数据集抓取、存储、读取、时间序列、数据处理等,涉及R包有 zoo(时间序列对象), xts(时间序列处理), timeSeries(Rmetrics系时间序列对象) timeDate(Rmetrics系时间序列处理), data.table(数据处理), quantmod(数据下载和图形可视化), RQuantLib(QuantLib数据接口), WindR(Wind数据接口), RJDBC(数据库访问接口), rhadoop(Hadoop访问接口), rhive(Hive访问接口), rredis(Redis访问接口), rmongodb(MongoDB访问接口), SparkR(Spark访问接口),fImport(Rmetrics系数据访问接口)等。
  • 指标计算:包括金融市场的技术指标的各种计算方法,涉及R包有 TTR(技术指标), TSA(时间序列计算), urca(单位根检验), fArma(Rmetrics系ARMA计算), fAsianOptions(Rmetrics系亚洲期权定价), fBasics(Rmetrics系计算工具), fCopulae(Rmetrics系财务分析), fExoticOptions(Rmetrics系期权计算), fGarch(Rmetrics系Garch模型), fNonlinear(Rmetrics系非线模型), fOptions(Rmetrics系期权定价), fRegression(Rmetrics系回归分析), fUnitRoots(Rmetrics系单位根检验) 等。
  • 回测交易:包括金融数据建模,并验证用历史数据验证模型的可靠性,涉及R包有 FinancialInstrument(金融产品), quantstrat(策略模型和回测), blotter(账户管理), fTrading(Rmetrics系交易分析)等。
  • 投资组合:对多策略或多模型进行管理和优化,涉及R包有 PortfolioAnalytics(组合分析和优化), stockPortfolio(股票组合管理), fAssets(Rmetrics系组合管理)等
  • 风险管理:对持仓进行风险指标的计算和风险提示,涉及R包有 PerformanceAnalytics(风险分析),fPortfolio(Rmetrics系组合优化), fExtremes(Rmetrics系数据处理)等。

基于上文中列出的R包,我们可以选择使用独立地第三方R包来构建我们的量化交易的系统,也可以选用完整的Rmetrics体系来构建量化交易的系统。这两类R包也可以混合使用,如果在混用时,由于他们基于的时间序列的底层对象是不一样的,那么类型转换的时候,可以你需要花点功夫处理一下。

上文中列出的R语言,并不是所有的R语言量化投资的R包,仅仅我关注的一些包。还有很多其他的,比如用于配对交易的包PairTrading;在Github上发布的,我并没有发现的R包等。

对于我自己来说,倾向于用独立地第三方R包来做量化交易系统,会用到其中的几个独立的R包。这样选择的主要原因有2个,一是中国市场比较特别,很多规则并不完全符合世界的标准。比如,股票T+1交易就是全球唯一的。另外一点是第三方的开源包,有一些可能有错误,所以你不应该把程序完全依赖于第三方包,要有独立的思考和判断,第三方包只是给我们提供了便利性。

那么常用的第三方R包的组合为:zoo, xts, TTR, quantmod, FinancialInstrument, quantstrat, blotter, PortfolioAnalytics, PerformanceAnalytics。这其中的任何一个包,都可以被替换或自己实现,从而保证自己量化交易系统的独特性。引用国外量化的教材上的一张图,国外用R来研究量化交易已经体系。

quantitative-analysis

图片摘自Introduction to Trading Systems,作者Guy Yollin。

本系列文章,稍后将对整个量化体系的金融R包进行全面的介绍,并加上我自己的理解。量化相关R包介绍的相关文章列表,持续更新中。。。

数据管理

策略模型

量化交易一条程序员可以利用技术优势,突破自己过上幸福生活的一条路,很艰难也很兴奋。我会一直坚持,希望路上的朋友一起加油!

转载请注明出处:
http://blog.fens.me/r-quant-packages/

打赏作者

R语言跨界调用C++

R的极客理想系列文章,涵盖了R的思想,使用,工具,创新等的一系列要点,以我个人的学习和体验去诠释R的强大。

R语言作为统计学一门语言,一直在小众领域闪耀着光芒。直到大数据的爆发,R语言变成了一门炙手可热的数据分析的利器。随着越来越多的工程背景的人的加入,R语言的社区在迅速扩大成长。现在已不仅仅是统计领域,教育,银行,电商,互联网….都在使用R语言。

要成为有理想的极客,我们不能停留在语法上,要掌握牢固的数学,概率,统计知识,同时还要有创新精神,把R语言发挥到各个领域。让我们一起动起来吧,开始R的极客理想。

关于作者:

  • 张丹(Conan), 程序员R,Nodejs,Java
  • weibo:@Conan_Z
  • blog: http://blog.fens.me
  • email: bsspirit@gmail.com

转载请注明出处:
http://blog.fens.me/r-cpp-rcpp

rcpp

前言

使用R语言已经很多年了,对很多的R包都已经了解,唯独没有碰和C++相关的部分,这可能很大的原因和我长期使用Java的背景有关。但随着多语言的发展,跨语言应用的流行,打通各语言界限的方法也已经是成熟。让R和C++实现通信,已经变得很简单。

跟上跨语言的步伐,打开R和C++的通道,让C++来解决R性能的诟病吧。

目录

  1. Rcpp的简单介绍
  2. 5分钟上手
  3. 数据类型转换

1. Rcpp的简单介绍

Rcpp包是一个打通R语言和C++语言的通信组件包,提供了R语言和C++函数的相互调用。R语言和C++语言的数据类型通过Rcpp包进行完整的映射。

Rcpp的官方网站:https://cran.r-project.org/web/packages/Rcpp/index.html

本文做为入门教程,只是简单介绍,如何能打通R语言和C++的通信通道,并不做深入地探讨。R语言和其他语言也有类似的通信实现,R语言和JAVA的调用,请参考文章解惑rJava R与Java的高速通道;R语言和Nodejs的调用,请参考文章Nodejs与R跨平台通信

2. 5分钟上手

做为5分钟上手的教程,我们只讲例子不讲API。

本文的系统环境

  • Win10 64bit
  • R version 3.2.3 (2015-12-10)

由于Windows系统的环境下需要Rtools支持,所以要手动下载对应版本的Rtoosl包,下载地址。我的R语言版本是3.2.3,所以我需要安装Rtools33.exe。安装EXE程序就不多说了,双击完成即可。

rtools

下载Rcpp的程序包,进行安装,一行代码搞定。


> install.packages("Rcpp")
trying URL 'https://mirrors.tuna.tsinghua.edu.cn/CRAN/bin/windows/contrib/3.2/Rcpp_0.12.6.zip'
Content type 'application/zip' length 3221864 bytes (3.1 MB)
downloaded 3.1 MB

package ‘Rcpp’ successfully unpacked and MD5 sums checked
Warning in install.packages :
  cannot remove prior installation of package ‘Rcpp’

The downloaded binary packages are in
	C:\Users\tinkpad\AppData\Local\Temp\RtmpKkg8zo\downloaded_packages

2.1 从hello world开始

从一个简单程序hello world开始吧,让R语言程序调用C++中的hello()函数。我用需要新建2个文件,放在同一个目录中。

  • demo.cpp,C++程序的源文件
  • demo.r,R程序源文件

首先,编辑demo.cpp,定义hello()函数。


~ notepad demo.cpp

#include <Rcpp.h>
#include <string>  

using namespace std;
using namespace Rcpp;

//[[Rcpp::export]]
string hello(string name) {
  cout << "hello " << name << endl;  
  return name;
}

/*** R
hello('world')
hello('Conan')
*/

上面Rcpp的代码,我们可以从3部分来看。

  • #include和using部分: 为包引用和命名空间的声明。<Rcpp.h>和namespace Rcpp是必要要加载的,另外由于使用了string的类型作为参数和返回值,所以需要<string>和namespace std。
  • 功能函数部分: 们定义了一个 hello(string name) 函数,有一个参数是string类型,返回值也为string类型。需要强调的是,对R开放的函数必须增加 //[[Rcpp::export]] 的注释声明
  • 代码执行: 用/*** R 和 */ 包含的部分,为R语言的代码,会默认被执行。

编辑demo.r,用来调用demo.cpp的hello()函数。


~ notepad demo.r

library(Rcpp)

sourceCpp(file='demo.cpp')
hello('R')

执行R语言的代码


# 加载Rcpp包
> library(Rcpp)

# 编译和加载demo.cpp文件
> sourceCpp(file='demo.cpp')

# 执行封装在demo.cpp中的R代码
> hello('world')
hello world
[1] "world"

> hello('Conan')
hello Conan
[1] "Conan"

# 执行hello函数
> hello('R')
hello [1]R
 "R"

一个非常简单的helloworld程序,就这样子完成了。

2.2 R和Rcpp的混写代码

上面2行代码,就完成了R对C++程序的调用,sourceCpp()函数真是强大。其实,sourceCpp()函数还提供了一种代码混写的方法,就是在R的代码中,直接嵌入C++代码。


sourceCpp(code='
  #include >Rcpp.h<
  #include >string<
  
  using namespace std;
  using namespace Rcpp;
  
  //[[Rcpp::export]]
  string hello(string name) {
    cout << "hello " << name << endl;  
    return name;
  }
')
hello('R2')

运行代码


> sourceCpp(code='
+   #include >Rcpp.h<
+   #include >string<
+   
+   using namespace std;
+   using namespace Rcpp;
+   
+   //[[Rcpp::export]]
+   string hello(string name) {
+     cout << "hello " << name << endl;  
+     return name;
+   }
+ ')

> hello('R2')
hello R2
[1] "R2"

这种多语言混写的语法虽然不太推荐,但对于这只有几行代码来说,还是很方便的。

2.2 用RStudioIDE生成cpp文件

如果你使用的RStudio IDE,开发起来将会非常方便,可以直接新建C++程序,生成一段标准的代码模板。

rstudio-cpp

生成的代码模板如下


#include <Rcpp.h>
using namespace Rcpp;

// This is a simple example of exporting a C++ function to R. You can
// source this function into an R session using the Rcpp::sourceCpp 
// function (or via the Source button on the editor toolbar). Learn
// more about Rcpp at:
//
//   http://www.rcpp.org/
//   http://adv-r.had.co.nz/Rcpp.html
//   http://gallery.rcpp.org/
//

// [[Rcpp::export]]
NumericVector timesTwo(NumericVector x) {
  return x * 2;
}

// You can include R code blocks in C++ files processed with sourceCpp
// (useful for testing and development). The R code will be automatically 
// run after the compilation.
//

/*** R
timesTwo(42)
*/

通过RStudio可以快速生成一段标准的代码模板,改改马上就能用了。

3. 数据类型转换

上面的例子中,我们测试了字符串类型的调用。R语言有多种的数据类型,我接下来都测试一下!

3.1 基本类型

基本类型,C++对应R语言的默认映射关系。C++的代码部分,如下所示:


// [[Rcpp::export]]
char char_type(char x){
  return x;
}

// [[Rcpp::export]]
int int_type(int x){
  return x;
}

// [[Rcpp::export]]
double double_type(double x){
  return x;
}

// [[Rcpp::export]]
bool bool_type(bool x){
  return x;
}

// [[Rcpp::export]]
void void_return_type(){
  Rprintf( "return void" );
}

执行R语言调用


# char类型
> a1<-char_type('a')
> a1;class(a1)         # 默认对应R的character类型
[1] "a"
[1] "character"
> char_type('bbii')    # 只处理字符串的第一个字节
[1] "b"

# int类型 
> a2<-int_type(111)
> a2;class(a2)         # 默认对应R的integer类型
[1] 111
[1] "integer" 
> int_type(111.1)      # 直接去掉小数位
[1] 111

# double类型 
> a3<-double_type(111.1)
> a3;class(a3)         # 默认对应R的numeric类型
[1] 111.1
[1] "numeric"
> double_type(111)
[1] 111

# boolean类型 
> a4<-bool_type(TRUE)
> a4;class(a4)        # 默认对应R的logical类型
[1] TRUE
[1] "logical"
> bool_type(0)        # 0为FALSE
[1] FALSE
> bool_type(1)        # 非0为TRUE
[1] TRUE

# 无参数无返回值 的函数
> a5<-void_return_type()
return void
> a5;class(a5)         # 返回值为NULL
NULL
[1] "NULL"

3.2 向量类型

向量类型,C++对应R语言的默认映射关系。C++的代码部分,如下所示:


// [[Rcpp::export]]
CharacterVector CharacterVector_type(CharacterVector x){
  return x;
}

// [[Rcpp::export]]
StringVector StringVector_type(StringVector x){
  return x;
}

// [[Rcpp::export]]
NumericVector NumericVector_type(NumericVector x) {
  return x;
}

// [[Rcpp::export]]
IntegerVector IntegerVector_type(IntegerVector x) {
  return x;
}

// [[Rcpp::export]]
DoubleVector DoubleVector_type(DoubleVector x){
  return x;
}

// [[Rcpp::export]]
LogicalVector LogicalVector_type(LogicalVector x){
  return x;
}

// [[Rcpp::export]]
DateVector DateVector_type(DateVector x){
  return x;
}

// [[Rcpp::export]]
DatetimeVector DatetimeVector_type(DatetimeVector x){
  return x;
}

执行R语言调用


# Character向量
> a6<-CharacterVector_type(c('abc','12345'))     
> a6;class(a6)                                    # 默认对应R的character类型
[1] "abc"   "12345"
[1] "character" 
> CharacterVector_type(c('abc',123.5, NA, TRUE))  # NA不处理
[1] "abc"   "123.5" NA      "TRUE" 

# String向量,完全同Character向量
> a7<-StringVector_type(c('abc','12345'))
> a7;class(a7)                                    # 默认对应R的character类型
[1] "abc"   "12345"
[1] "character"
> StringVector_type(c('abc',123.5, NA, TRUE))
[1] "abc"   "123.5" NA      "TRUE" 

# Numeric向量
> a8<-NumericVector_type(rnorm(5))
> a8;class(a8)                                    # 默认对应R的numeric类型
[1] -0.2813472 -0.2235722 -0.6958443 -1.5322172  0.5004307
[1] "numeric"
> NumericVector_type(c(rnorm(5),NA,TRUE))         # NA不处理,TRUE为1
[1]  0.1700925  0.5169612 -0.3622637  1.0763204 -0.5729958
[6]         NA  1.0000000

# Integer向量
> a9<-IntegerVector_type(c(11,9.9,1.2))           # 直接去掉小数位
> a9;class(a9)                                    # 默认对应R的integer类型
[1] 11  9  1
[1] "integer"
> IntegerVector_type(c(11,9.9,1.2,NA,TRUE))       # NA不处理,TRUE为1
[1] 11  9  1 NA  1

# Double向量,同Numeric向量
> a10<-DoubleVector_type(rnorm(5))
> a10;class(a10)                                  # 默认对应R的numeric类型
[1]  0.9400947 -0.8976913  0.2744319 -1.5278219  1.2010569
[1] "numeric"
> DoubleVector_type(c(rnorm(5),NA,TRUE))          # NA不处理,TRUE为1
[1]  2.0657148  0.2810003  2.1080900 -1.2783693  0.2198551
[6]         NA  1.0000000

# Logical向量 
> a11<-LogicalVector_type(c(TRUE,FALSE))
> a11;class(a11)                                  # 默认对应R的logical类型
[1]  TRUE FALSE
[1] "logical"
> LogicalVector_type(c(TRUE,FALSE,TRUE,0,-1,NA))  # NA不处理,0为FALSE, 非0为TRUE
[1]  TRUE FALSE  TRUE FALSE  TRUE    NA

 # Date向量 
> a12<-DateVector_type(c(Sys.Date(),as.Date('2016-10-10')))
> a12;class(a12)                                  # 默认对应R的Date类型
[1] "2016-08-01" "2016-10-10"
[1] "Date"
> DateVector_type(c(Sys.Date(),as.Date('2016-10-10'),NA,TRUE,FALSE))   # NA不处理,TRUE为1970-01-02, FALSE为1970-01-01
[1] "2016-08-01" "2016-10-10" NA           "1970-01-02"
[5] "1970-01-01"

 # Datetime向量 
> a13<-DatetimeVector_type(c(Sys.time(),as.POSIXct('2016-10-10')))
> a13;class(a13)                                  # 默认对应R的POSIXct类型
[1] "2016-08-01 20:05:25 CST" "2016-10-10 00:00:00 CST"
[1] "POSIXct" "POSIXt" 
> DatetimeVector_type(c(Sys.time(),as.POSIXct('2016-10-10'),NA,TRUE,FALSE))  # NA不处理
[1] "2016-08-01 20:05:25 CST" "2016-10-10 00:00:00 CST"
[3] NA                        "1970-01-01 08:00:01 CST"
[5] "1970-01-01 08:00:00 CST"

3.3 矩阵类型

矩阵类型,C++对应R语言的默认映射关系。C++的代码部分,如下所示:


// [[Rcpp::export]]
CharacterMatrix CharacterMatrix_type(CharacterMatrix x){
  return x;
}

// [[Rcpp::export]]
StringMatrix StringMatrix_type(StringMatrix x){
  return x;
}

// [[Rcpp::export]]
NumericMatrix NumericMatrix_type(NumericMatrix x){
  return x;
}

// [[Rcpp::export]]
IntegerMatrix IntegerMatrix_type(IntegerMatrix x){
  return x;
}

// [[Rcpp::export]]
LogicalMatrix LogicalMatrix_type(LogicalMatrix x){
  return x;
}

// [[Rcpp::export]]
ListMatrix ListMatrix_type(ListMatrix x){
  return x;
}

执行R语言调用


# Character矩阵
> a14<-CharacterMatrix_type(matrix(LETTERS[1:20],ncol=4))
> a14;class(a14)
     [,1] [,2] [,3] [,4]
[1,] "A"  "F"  "K"  "P" 
[2,] "B"  "G"  "L"  "Q" 
[3,] "C"  "H"  "M"  "R" 
[4,] "D"  "I"  "N"  "S" 
[5,] "E"  "J"  "O"  "T" 
[1] "matrix"                        

# String矩阵,同Character矩阵
> a15<-StringMatrix_type(matrix(LETTERS[1:20],ncol=4))
> a15;class(a15)
     [,1] [,2] [,3] [,4]
[1,] "A"  "F"  "K"  "P" 
[2,] "B"  "G"  "L"  "Q" 
[3,] "C"  "H"  "M"  "R" 
[4,] "D"  "I"  "N"  "S" 
[5,] "E"  "J"  "O"  "T" 
[1] "matrix"

# Numeric矩阵
> a16<-NumericMatrix_type(matrix(rnorm(20),ncol=4))
> a16;class(a16)
           [,1]       [,2]       [,3]       [,4]
[1,]  1.2315498  2.3234269  0.5974143  0.9072356
[2,]  0.3484811  0.3814024 -0.2018324  0.8717205
[3,] -0.2025285  2.1076947 -0.3433948  1.1523710
[4,] -1.4948252 -0.7724951 -0.7681800 -0.5406494
[5,]  0.4815904  1.4930873 -1.1444258  0.2537099
[1] "matrix"

# Integer矩阵 
> a17<-IntegerMatrix_type(matrix(seq(1,10,length.out = 20),ncol=4))
> a17;class(a17)
     [,1] [,2] [,3] [,4]
[1,]    1    3    5    8
[2,]    1    3    6    8
[3,]    1    4    6    9
[4,]    2    4    7    9
[5,]    2    5    7   10
[1] "matrix"

# Logical矩阵 
> a18<-LogicalMatrix_type(matrix(c(rep(TRUE,5),rep(FALSE,5),rnorm(10)),ncol=4))
> a18;class(a18)
     [,1]  [,2] [,3] [,4]
[1,] TRUE FALSE TRUE TRUE
[2,] TRUE FALSE TRUE TRUE
[3,] TRUE FALSE TRUE TRUE
[4,] TRUE FALSE TRUE TRUE
[5,] TRUE FALSE TRUE TRUE
[1] "matrix"

# List矩阵,支持多类型的矩阵
> a19<-ListMatrix_type(matrix(rep(list(a=1,b='2',c=NA,d=TRUE),10),ncol=5))
> a19;class(a19)
     [,1] [,2] [,3] [,4] [,5]
[1,] 1    1    1    1    1   
[2,] "2"  "2"  "2"  "2"  "2" 
[3,] NA   NA   NA   NA   NA  
[4,] TRUE TRUE TRUE TRUE TRUE
[5,] 1    1    1    1    1   
[6,] "2"  "2"  "2"  "2"  "2" 
[7,] NA   NA   NA   NA   NA  
[8,] TRUE TRUE TRUE TRUE TRUE
[1] "matrix"

3.4 其他数据类型

其他数据类型包括了,R语言特有的数据类型数据框(data.frame),环境空间(Environment)S3,S4,RC等的对象类型。


// [[Rcpp::export]]
Date Date_type(Date x){
  return x;
}

// [[Rcpp::export]]
Datetime Datetime_type(Datetime x){
  return x;
}

// [[Rcpp::export]]
S4 S4_type(S4 x){
  return x;
}

// [[Rcpp::export]]
RObject RObject_type(RObject x){
  return x;
}

// [[Rcpp::export]]
SEXP SEXP_type(SEXP x){
  return x;
}

// [[Rcpp::export]]
Environment Environment_type(Environment x){
  return x;
}

执行R语言调用


# data.frame类型
> a19<-DataFrame_type(data.frame(a=rnorm(3),b=1:3))
> a19;class(a19)
           a b
1 -1.8844994 1
2  0.6053935 2
3 -0.7693985 3
[1] "data.frame"

# list类型 
> a20<-List_type(list(a=1,b='2',c=NA,d=TRUE))
> a20;class(a20)
$a
[1] 1
$b
[1] "2"
$c
[1] NA
$d
[1] TRUE
[1] "list"

# Date类型
> a21<-Date_type(Sys.Date())
> a21;class(a21)
[1] "2016-08-01"
[1] "Date"
> Date_type(Sys.time())                # 不能正确处理POSIXct类型的数据
[1] "4026842-05-26"

# POSIXct类型 
> a22<-Datetime_type(Sys.time())
> a22;class(a22)
[1] "2016-08-01 20:27:37 CST"
[1] "POSIXct" "POSIXt" 
> Datetime_type(Sys.Date())            # 不能正确处理Date类型的数据
[1] "1970-01-01 12:43:34 CST"

# S3面向对象类型,对应S4的类型定义
> setClass("Person",slots=list(name="character",age="numeric"))
> s4<-new("Person",name="F",age=44)
> a23<-S4_type(s4)
> a23;class(a23)
An object of class "Person"
Slot "name":
[1] "F"
Slot "age":
[1] 44
[1] "Person"
attr(,"package")
[1] ".GlobalEnv"

# S3面向对象类型 ,没有对应的类型,通过RObject来传值
> s3<-structure(2, class = "foo")
> a24<-RObject_type(s3)
> a24;class(a24)
[1] 2
attr(,"class")
[1] "foo"
[1] "foo"

# RObject也可以处理S4对象
> a25<-RObject_type(s4)
> a25;class(a25)
An object of class "Person"
Slot "name":
[1] "F"
Slot "age":
[1] 44
[1] "Person"
attr(,"package")
[1] ".GlobalEnv"

# RObject也可以处理RC对象 
> User<-setRefClass("User",fields=list(name="character"))
> rc<-User$new(name="u1")
> a26<-RObject_type(rc)
> a26;class(a26)
Reference class object of class "User"
Field "name":
[1] "u1"
[1] "User"
attr(,"package")
[1] ".GlobalEnv"

# RObject也可以处理function类型
> a27<-RObject_type(function(x) x+2)
> a27;class(a27)
function(x) x+2
[1] "function"

# environment类型
> a28<-Environment_type(new.env())
> a28;class(a28)
<environment: 0x0000000015350a80>
[1] "environment"

# SEXP为任意类型,通过具体调用时再进行类型判断
> SEXP_type('fdafdaa')
[1] "fdafdaa"

> SEXP_type(rc)
Reference class object of class "User"
Field "name":
[1] "u1"

> SEXP_type(data.frame(a=rnorm(3),b=1:3))
           a b
1 -0.5396140 1
2  0.1694799 2
3 -1.8818596 3

> SEXP_type(function(x) x+2)
function(x) x+2

最后总结一下,R和Rcpp中类型对应的关系。

C++类型 R类型
char character
int integer
double numeric
bool logical
Rcpp::Date Date
Rcpp::Datetime POSIXct
Rcpp::CharacterVector character
Rcpp::StringVector character
Rcpp::NumericVector numeric
Rcpp::IntegerVector integer
Rcpp::DoubleVector numeric
Rcpp::LogicalVector logical
Rcpp::DateVector Date
Rcpp::DatetimeVector POSIXct
Rcpp::CharacterMatrix matrix
Rcpp::StringMatrix matrix
Rcpp::NumericMatrix matrix
Rcpp::IntegerMatrix matrix
Rcpp::LogicalMatrix matrix
Rcpp::ListMatrix matrix
Rcpp::DataFrame data.frame
Rcpp::List list
Rcpp::S4 S4
Rcpp::Environment environment
Rcpp::RObject 任意类型
Rcpp::SEXP 任意类型

本文简单地介绍了通过R语言Rcpp包调用C++程序的一种方法,调用的关键点就在于数据类型的匹配,而从保证R语言和C++之间的数据传输。从上面测试来看,R语言中的所有数据类型,都可以通过Rcpp包进行映射到C++的程序中。接下来,我们就可以根据自己的需求,把一些更关注的性能的程序放到C++中来实现,从而提高计算效率。

转载请注明出处:
http://blog.fens.me/r-cpp-rcpp

打赏作者

R语言解读多元线性回归模型

R的极客理想系列文章,涵盖了R的思想,使用,工具,创新等的一系列要点,以我个人的学习和体验去诠释R的强大。

R语言作为统计学一门语言,一直在小众领域闪耀着光芒。直到大数据的爆发,R语言变成了一门炙手可热的数据分析的利器。随着越来越多的工程背景的人的加入,R语言的社区在迅速扩大成长。现在已不仅仅是统计领域,教育,银行,电商,互联网….都在使用R语言。

要成为有理想的极客,我们不能停留在语法上,要掌握牢固的数学,概率,统计知识,同时还要有创新精神,把R语言发挥到各个领域。让我们一起动起来吧,开始R的极客理想。

关于作者:

  • 张丹(Conan), 程序员R,Nodejs,Java
  • weibo:@Conan_Z
  • blog: http://blog.fens.me
  • email: bsspirit@gmail.com

转载请注明出处:
http://blog.fens.me/r-multi-linear-regression

reg-multi-liner

前言

本文接上一篇R语言解读一元线性回归模型。在许多生活和工作的实际问题中,影响因变量的因素可能不止一个,比如对于知识水平越高的人,收入水平也越高,这样的一个结论。这其中可能包括了因为更好的家庭条件,所以有了更好的教育;因为在一线城市发展,所以有了更好的工作机会;所处的行业赶上了大的经济上行周期等。要想解读这些规律,是复杂的、多维度的,多元回归分析方法更适合解读生活的规律。

由于本文为非统计的专业文章,所以当出现与教课书不符的描述,请以教课书为准。本文力求用简化的语言,来介绍多元线性回归的知识,同时配合R语言的实现。

目录

  1. 多元线性回归介绍
  2. 元线性回归建模
  3. 模型优化
  4. 案例:黑色系期货日K线数据验证

1. 多元线性回归介绍

对比一元线性回归,多元线性回归是用来确定2个或2个以上变量间关系的统计分析方法。多元线性回归的基本的分析方法与一元线性回归方法是类似的,我们首先需要对选取多元数据集并定义数学模型,然后进行参数估计,对估计出来的参数进行显著性检验,残差分析,异常点检测,最后确定回归方程进行模型预测。

由于多元回归方程有多个自变量,区别于一元回归方程,有一项很重要的操作就是自变量的优化,挑选出相关性最显著的自变量,同时去除不显著的自变量。在R语言中,有很方便地用于优化函数,可以很好的帮助我们来改进回归模型。

下面就开始多元线性回归的建模过程。

2. 多元线性回归建模

做过商品期货研究的人,都知道黑色系品种是具有产业链上下游的关系。铁矿石是炼钢的原材料,焦煤和焦炭是炼钢的能源资源,热卷即热轧卷板是以板坯为原料经加热后制成的钢板,螺纹钢是表面带肋的钢筋。

由于有产业链的关系,假设我们想要预测螺纹钢的价格,那么影响螺纹钢价格的因素可以会涉及到原材料,能源资源和同类材料等。比如,铁矿石价格如果上涨,螺纹钢就应该要涨价了。

2.1 数据集和数学模型

先从数据开始介绍,这次的数据集,我选择的期货黑色系的品种的商品期货,包括了大连期货交易所的 焦煤(JM),焦炭(J),铁矿石(I),上海期货交易所的 螺纹钢(RU) 和 热卷(HC)。

数据集为2016年3月15日,当日白天开盘的交易数据,为黑色系的5个期货合约的分钟线的价格数据。


# 数据集已存在df变量中
> head(df,20)
                       x1    x2    x3   x4    y
2016-03-15 09:01:00 754.5 616.5 426.5 2215 2055
2016-03-15 09:02:00 752.5 614.5 423.5 2206 2048
2016-03-15 09:03:00 753.0 614.0 423.0 2199 2044
2016-03-15 09:04:00 752.5 613.0 422.5 2197 2040
2016-03-15 09:05:00 753.0 615.5 424.0 2198 2043
2016-03-15 09:06:00 752.5 614.5 422.0 2195 2040
2016-03-15 09:07:00 752.0 614.0 421.5 2193 2036
2016-03-15 09:08:00 753.0 615.0 422.5 2197 2043
2016-03-15 09:09:00 754.0 615.5 422.5 2197 2041
2016-03-15 09:10:00 754.5 615.5 423.0 2200 2044
2016-03-15 09:11:00 757.0 616.5 423.0 2201 2045
2016-03-15 09:12:00 756.0 615.5 423.0 2200 2044
2016-03-15 09:13:00 755.5 615.0 423.0 2197 2042
2016-03-15 09:14:00 755.5 615.0 423.0 2196 2042
2016-03-15 09:15:00 756.0 616.0 423.5 2200 2045
2016-03-15 09:16:00 757.5 616.0 424.0 2205 2052
2016-03-15 09:17:00 758.5 618.0 424.0 2204 2051
2016-03-15 09:18:00 759.5 618.5 424.0 2205 2053
2016-03-15 09:19:00 759.5 617.5 424.5 2206 2053
2016-03-15 09:20:00 758.5 617.5 423.5 2201 2050

数据集包括有6列:

  • 索引, 为时间
  • x1, 为焦炭(j1605)合约的1分钟线的报价数据
  • x2, 为焦煤(jm1605)合约的1分钟线的报价数据
  • x3, 为铁矿石(i1605)合约的1分钟线的报价数据
  • x4, 为热卷(hc1605)合约的1分钟线的报价数据
  • y, 为螺纹钢(rb1605)合约的1分钟线的报价数据

假设螺纹钢的价格与其他4个商品的价格有线性关系,那么我们建立以螺纹钢为因变量,以焦煤、焦炭、铁矿石和热卷的为自变量的多元线性回归模型。用公式表示为:

y = a + b * x1 + c * x2 + d * x3 + e * x4 + ε
  • y,为因变量,螺纹钢
  • x1,为自变量,焦煤
  • x2,为自变量,焦炭
  • x3,为自变量,铁矿石
  • x4,为自变量,热卷
  • a,为截距
  • b,c,d,e,为自变量系数
  • ε, 为残差,是其他一切不确定因素影响的总和,其值不可观测。假定ε服从正态分布N(0,σ^2)。

通过对多元线性回归模型的数学定义,接下来让我们利用数据集做多元回归模型的参数估计。

2.2. 回归参数估计

上面公式中,回归参数 a, b, c, d,e都是我们不知道的,参数估计就是通过数据来估计出这些参数,从而确定自变量和因变量之前的关系。我们的目标是要计算出一条直线,使直线上每个点的Y值和实际数据的Y值之差的平方和最小,即(Y1实际-Y1预测)^2+(Y2实际-Y2预测)^2+ …… +(Yn实际-Yn预测)^2 的值最小。参数估计时,我们只考虑Y随X自变量的线性变化的部分,而残差ε是不可观测的,参数估计法并不需要考虑残差。

类似于一元线性回归,我们用R语言来实现对数据的回归模型的参数估计,用lm()函数来实现多元线性回归的建模过程。


# 建立多元线性回归模型
> lm1<-lm(y~x1+x2+x3+x4,data=df)

# 打印参数估计的结果
> lm1

Call:
lm(formula = y ~ x1 + x2 + x3 + x4, data = df)

Coefficients:
(Intercept)           x1           x2           x3           x4  
   212.8780       0.8542       0.6672      -0.6674       0.4821  

这样我们就得到了y和x关系的方程。

y = 212.8780 + 0.8542 * x1 + 0.6672 * x2 - 0.6674 * x3 + 0.4821 * x4

2.3. 回归方程的显著性检验

参考一元线性回归的显著性检验,多元线性回归的显著性检验,同样是需要经过 T检验,F检验,和R^2(R平方)相关系统检验。在R语言中这三种检验的方法都已被实现,我们只需要把结果解读,我们可以summary()函数来提取模型的计算结果。


> summary(lm1)

Call:
lm(formula = y ~ x1 + x2 + x3 + x4, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.9648 -1.3241 -0.0319  1.2403  5.4194 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 212.87796   58.26788   3.653 0.000323 ***
x1            0.85423    0.10958   7.795 2.50e-13 ***
x2            0.66724    0.12938   5.157 5.57e-07 ***
x3           -0.66741    0.15421  -4.328 2.28e-05 ***
x4            0.48214    0.01959  24.609  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.028 on 221 degrees of freedom
Multiple R-squared:  0.9725,	Adjusted R-squared:  0.972 
F-statistic:  1956 on 4 and 221 DF,  p-value: < 2.2e-16
  • T检验:所自变量都是非常显著***
  • F检验:同样是非常显著,p-value < 2.2e-16
  • 调整后的R^2:相关性非常强为0.972

最后,我们通过的回归参数的检验与回归方程的检验,得到最后多元线性回归方程为:


y = 212.87796 + 0.85423 * x1 + 0.66724 * x2 - 0.66741 * x3 + 0.48214 * x4

即

螺纹钢 = 212.87796 + 0.85423 * 焦炭 + 0.66724 * 焦煤 - 0.66741 * 铁矿石 + 0.48214 * 热卷

2.4 残差分析和异常点检测

在得到的回归模型进行显著性检验后,还要在做残差分析(预测值和实际值之间的差),检验模型的正确性,残差必须服从正态分布N(0,σ^2)。直接用plot()函数生成4种用于模型诊断的图形,进行直观地分析。


> par(mfrow=c(2,2))
> plot(lm1)

m01

  • 残差和拟合值(左上),残差和拟合值之间数据点均匀分布在y=0两侧,呈现出随机的分布,红色线呈现出一条平稳的曲线并没有明显的形状特征。
  • 残差QQ图(右上),数据点按对角直线排列,趋于一条直线,并被对角直接穿过,直观上符合正态分布。
  • 标准化残差平方根和拟合值(左下),数据点均匀分布在y=0两侧,呈现出随机的分布,红色线呈现出一条平稳的曲线并没有明显的形状特征。
  • 标准化残差和杠杆值(右下),没有出现红色的等高线,则说明数据中没有特别影响回归结果的异常点。

结论,没有明显的异常点,残差符合假设条件。

2.5. 模型预测

我们得到了多元线性回归方程的公式,就可以对数据进行预测了。我们可以用R语言的predict()函数来计算预测值y0和相应的预测区间,并把实际值和预测值一起可视化化展示。


> par(mfrow=c(1,1))  #设置画面布局

# 预测计算
> dfp<-predict(lm1,interval="prediction")

# 打印预测时
> head(dfp,10)
                fit      lwr      upr
2014-03-21 3160.526 3046.425 3274.626
2014-03-24 3193.253 3078.868 3307.637
2014-03-25 3240.389 3126.171 3354.607
2014-03-26 3228.565 3114.420 3342.710
2014-03-27 3222.528 3108.342 3336.713
2014-03-28 3262.399 3148.132 3376.666
2014-03-31 3291.996 3177.648 3406.344
2014-04-01 3305.870 3191.447 3420.294
2014-04-02 3275.370 3161.018 3389.723
2014-04-03 3297.358 3182.960 3411.755

# 合并数据
> mdf<-merge(df$y,dfp)	 

# 画图
> draw(mdf)

m02

图例说明

  • y, 实际价格,红色线
  • fit, 预测价格,绿色线
  • lwr,预测最低价,蓝色线
  • upr,预测最高价,紫色线

从图中看出,实际价格y和预测价格fit,在大多数的时候都是很贴近的。我们的一个模型就训练好了!

3. 模型优化

上文中,我们已经很顺利的发现了一个非常不错的模型。如果要进行模型优化,可以用R语言中update()函数进行模型的调整。我们首先检查一下每个自变量x1,x2,x3,x4和因变量y之间的关系。

pairs(as.data.frame(df))

m03

从图中,我们可以发现x2与Y的关系,可能是最偏离线性的。那么,我们尝试对多元线性回归模型进行调整,从原模型中去掉x2变量。



# 模型调整
> lm2<-update(lm1, .~. -x2)

> summary(lm2)

Call:
lm(formula = y ~ x1 + x3 + x4, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.0039 -1.3842  0.0177  1.3513  4.8028 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 462.47104   34.26636   13.50  < 2e-16 ***
x1            1.08728    0.10543   10.31  < 2e-16 ***
x3           -0.40788    0.15394   -2.65  0.00864 ** 
x4            0.42582    0.01718   24.79  < 2e-16 ***
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.142 on 222 degrees of freedom
Multiple R-squared:  0.9692,	Adjusted R-squared:  0.9688 
F-statistic:  2330 on 3 and 222 DF,  p-value: < 2.2e-16

当把自变量x2去掉后,自变量x3的T检验反而变大了,同时Adjusted R-squared变小了,所以我们这次调整是有问题的。

如果通过生产和原材料的内在逻辑分析,焦煤与焦炭属于上下游关系。焦煤是生产焦炭的一种原材料,焦炭是焦煤与其他炼焦煤经过配煤焦化形成的产品,一般生产 1 吨焦炭需要1.33 吨炼焦煤,其中焦煤至少占 30% 。

我们把焦煤 和 焦炭的关系改变一下,增加x1*x2的关系匹配到模型,看看效果。


# 模型调整
> lm3<-update(lm1, .~. + x1*x2)
> summary(lm3)

Call:
lm(formula = y ~ x1 + x2 + x3 + x4 + x1:x2, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.8110 -1.3501 -0.0595  1.2019  5.3884 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) 7160.32231 7814.50048   0.916    0.361    
x1            -8.45530   10.47167  -0.807    0.420    
x2           -10.58406   12.65579  -0.836    0.404    
x3            -0.64344    0.15662  -4.108 5.63e-05 ***
x4             0.48363    0.01967  24.584  < 2e-16 ***
x1:x2          0.01505    0.01693   0.889    0.375    
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.029 on 220 degrees of freedom
Multiple R-squared:  0.9726,	Adjusted R-squared:  0.972 
F-statistic:  1563 on 5 and 220 DF,  p-value: < 2.2e-16

从结果中发现,增加了x1*x2列后,原来的x1,x2和Intercept的T检验都不显著。继续调整模型,从模型中去掉x1,x2两个自变量。


# 模型调整
> lm4<-update(lm3, .~. -x1-x2)
> summary(lm4)

Call:
lm(formula = y ~ x3 + x4 + x1:x2, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.9027 -1.2516 -0.0167  1.2748  5.8683 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  6.950e+02  1.609e+01  43.183  < 2e-16 ***
x3          -6.284e-01  1.530e-01  -4.108 5.61e-05 ***
x4           4.959e-01  1.785e-02  27.783  < 2e-16 ***
x1:x2        1.133e-03  9.524e-05  11.897  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.035 on 222 degrees of freedom
Multiple R-squared:  0.9722,	Adjusted R-squared:  0.9718 
F-statistic:  2588 on 3 and 222 DF,  p-value: < 2.2e-16

从调整后的结果来看,效果还不错。不过,也并没有比最初的模型有所提高。

对于模型调整的过程,如果我们手动调整测试时,一般都会基于业务知识来操作。如果是按照数据指标来计算,我们可以用R语言中提供的逐步回归的优化方法,通过AIC指标来判断是否需要参数优化。


#对lm1模型做逐步回归
> step(lm1)
Start:  AIC=324.51
y ~ x1 + x2 + x3 + x4

       Df Sum of Sq    RSS    AIC
               908.8 324.51
- x3    1     77.03  985.9 340.90
- x2    1    109.37 1018.2 348.19
- x1    1    249.90 1158.8 377.41
- x4    1   2490.56 3399.4 620.65

Call:
lm(formula = y ~ x1 + x2 + x3 + x4, data = df)

Coefficients:
(Intercept)           x1           x2           x3           x4  
   212.8780       0.8542       0.6672      -0.6674       0.4821  

通过计算AIC指标,lm1的模型AIC最小时为324.51,每次去掉一个自变量都会让AIC的值变大,所以我们还是不调整比较好。

对刚才的lm3模型做逐步回归的模型调整。


#对lm3模型做逐步回归
> step(lm3)
Start:  AIC=325.7               #当前AIC
y ~ x1 + x2 + x3 + x4 + x1:x2

        Df Sum of Sq    RSS    AIC
- x1:x2  1      3.25  908.8 324.51
                905.6 325.70
- x3     1     69.47  975.1 340.41
- x4     1   2487.86 3393.5 622.25

Step:  AIC=324.51               #去掉x1*x2项的AIC
y ~ x1 + x2 + x3 + x4

       Df Sum of Sq    RSS    AIC
               908.8 324.51
- x3    1     77.03  985.9 340.90
- x2    1    109.37 1018.2 348.19
- x1    1    249.90 1158.8 377.41
- x4    1   2490.56 3399.4 620.65

Call:
lm(formula = y ~ x1 + x2 + x3 + x4, data = df)

Coefficients:
(Intercept)           x1           x2           x3           x4  
   212.8780       0.8542       0.6672      -0.6674       0.4821  

通过AIC的判断,去掉X1*X2项后AIC最小,最后的检验结果告诉我们,还是原初的模型是最好的。

4. 案例:黑色系期货日K线数据验证

最后,我们用上面5个期货合约的日K线数据测试一下,找到多元回归关系。


> lm9<-lm(y~x1+x2+x3+x4,data=df)  # 日K线数据
> summary(lm9)

Call:
lm(formula = y ~ x1 + x2 + x3 + x4, data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-173.338  -37.470    3.465   32.158  178.982 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 386.33482   31.07729  12.431  < 2e-16 ***
x1            0.75871    0.07554  10.045  < 2e-16 ***
x2           -0.62907    0.14715  -4.275 2.24e-05 ***
x3            1.16070    0.05224  22.219  < 2e-16 ***
x4            0.46461    0.02168  21.427  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 57.78 on 565 degrees of freedom
Multiple R-squared:  0.9844,	Adjusted R-squared:  0.9843 
F-statistic:  8906 on 4 and 565 DF,  p-value: < 2.2e-16

数据集的基本统计信息。


> summary(df)
     Index                           x1               x2       
 Min.   :2014-03-21 00:00:00   Min.   : 606.5   Min.   :494.0  
 1st Qu.:2014-10-21 06:00:00   1st Qu.: 803.5   1st Qu.:613.1  
 Median :2015-05-20 12:00:00   Median : 939.0   Median :705.8  
 Mean   :2015-05-21 08:02:31   Mean   : 936.1   Mean   :695.3  
 3rd Qu.:2015-12-16 18:00:00   3rd Qu.:1075.0   3rd Qu.:773.0  
 Max.   :2016-07-25 00:00:00   Max.   :1280.0   Max.   :898.0  

       x3              x4             y       
 Min.   :284.0   Min.   :1691   Min.   :1626  
 1st Qu.:374.1   1st Qu.:2084   1st Qu.:2012  
 Median :434.0   Median :2503   Median :2378  
 Mean   :476.5   Mean   :2545   Mean   :2395  
 3rd Qu.:545.8   3rd Qu.:2916   3rd Qu.:2592  
 Max.   :825.0   Max.   :3480   Max.   :3414  

m04

对于日K线数据,黑色系的5个品种,同样具有非常强的相关关系,那么我们就可以把这个结论应用到实际的交易中了。

本文通过多元回归的统计分析方法,介绍多元回归在金融市场的基本应用。我们通过建立因变量和多个自变量的模型,从而发现生活中更复杂的规律,并建立有效的验证指标。让我们我们的技术优势,去金融市场抢钱吧。

转载请注明出处:
http://blog.fens.me/r-multi-linear-regression

打赏作者