• Archive by category "程序算法"

Blog Archives

线性和非线性的最小二乘回归

算法为王系列文章,涵盖了计算机算法,数据挖掘(机器学习)算法,统计算法,金融算法等的多种跨学科算法组合。在大数据时代的背景下,算法已经成为了金字塔顶的明星。一个好的算法可以创造一个伟大帝国,就像Google。

算法为王的时代正式到来….

关于作者:

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

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

前言

很多经典的机器学习模型都是在最小二乘法的基础上发展起来的,理解最小二乘法的原理,对于我们掌握机器学习的高级算法是非常重要的,尤其是回归模型中。最小二乘法是主要是用来做函数拟合,或者求函数极值等问题。

本文以回归模型为例,介绍在最小二乘法在线性和非线性模型中的应用。

目录

  1. 最小二乘法介绍
  2. 最小二乘求解
  3. 线性最小二乘回归
  4. 非线性最小二乘回归

1. 最小二乘法介绍

最小二乘法(Least Squares Method,简记为LSE),源于天文学和测地学上的应用需要,是勒让德( A. M. Legendre)于1805年在其著作《计算慧星轨道的新方法》中提出的,它的主要思想是通过计算理论值与观测值之差的平方和最小,求解未知的参数。

最小二乘法的公式为:

公式解读:

  • H, 目标值,观测值与理论值之差的平方和
  • H’, 最小化的H值
  • y, 观测值
  • yi, 理论值
  • argmin, 最小化函数

我们先了把术语进行统一定义,观测值通常说的是我们的样本数据,理论值就是假设拟合函数产生的预测数据,目标函数是判断一个模型的好不好的标准,进行最小二乘法建立目标函数,让目标函数最小化。

以上图为线性拟合为例,观测值为y1,y2,…,y7,理论值为y1′,y2′,….,y7’,是线上的点并与观测值的x相同,红色线就是我们要做的拟合直线,我们可以画出无数条这样的直线,判断标准是让这条线上的点的理论值与观测值的误差平方和最小,即最小二乘法。

SSE=sum((y1-y1')^2+(y2-y2')^2+(y3-y3')^2+(y4-y4')^2+(y5-y5')^2+(y6-y6')^2+(y7-y7')^2)

利用最小二乘法,可以对数据进行线性和非线性拟合,使误差平方和(SSE)或残差平方和最小。如果观测到的误差近似正态分布时,这种方法是非常有效的。我们在线性回归模型中,要进行残差的正态分布检验,就是基于这个目的的,可以参考文章R语言解读一元线性回归模型

2. 最小二乘求解

在线性方程的求解或是数据曲线拟合中,利用最小二乘法求得的解则被称为最小二乘解。最小二乘法(又称最小平方法)是一种数学优化技术,它通过最小化误差的平方和寻找数据的最佳函数匹配。其他一些优化问题也可通过最小化能量或最大化熵用最小二乘法来表达。

定义数据集x和y,两个变量。

> x<-c(6.19,2.51,7.29,7.01,5.7,2.66,3.98,2.5,9.1,4.2)
> y<-c(5.25,2.83,6.41,6.71,5.1,4.23,5.05,1.98,10.5,6.3)

可视化,画出x,y变量的散点图。

> plot(x,y)

对数据进行最小二乘求解,找到最小二乘解。


> fit<-lsfit(x,y);fit
$coefficients
Intercept         X
0.8310557 0.9004584

$residuals
[1] -1.1548933 -0.2612063 -0.9853975 -0.4332692 -0.8636686  1.0037250  0.6351198 -1.1022017
[9]  1.4747728  1.6870190

$intercept
[1] TRUE

$qr
$qt
[1] -17.1901414   6.2044421  -0.7047339  -0.1530724  -0.5856560   1.2766692   0.9102648
[8]  -0.8295242   1.7584540   1.9625308

$qr
Intercept           X
[1,] -3.1622777 -16.1718880
[2,]  0.3162278   6.8903149
[3,]  0.3162278  -0.2782874
[4,]  0.3162278  -0.2376506
[5,]  0.3162278  -0.0475287
[6,]  0.3162278   0.3936703
[7,]  0.3162278   0.2020970
[8,]  0.3162278   0.4168913
[9,]  0.3162278  -0.5409749
[10,]  0.3162278   0.1701682

$qraux
[1] 1.316228 1.415440

$rank
[1] 2

$pivot
[1] 1 2

$tol
[1] 1e-07

attr(,"class")
[1] "qr"

结果解读:

  • coefficients, 系数的最小二乘估计,包括截距Intercept和自变量X
  • residuals, 残差
  • intercept, 有截距
  • qt, 矩阵QR分解。QR分解法是将矩阵分解成一个正规正交矩阵与上三角形矩阵。关于矩阵的详细操作,请参考文章R语言中的矩阵计算

我们看到lsfit()函数使用矩阵的QR分解的方法,计算得出了最小二乘估计的截距Intercept和自变量系数。

目标:y=a*x+b,用最小二乘法估计出参数a和b。

根据最小二乘法的矩阵计算公式,我直接跳过了推到过程,得出结论。

公式解释:

  • X,为自变量矩阵
  • Y,为因变量矩阵
  • B,为参数,需要求解
  • XT,X的转置矩阵
  • *, 矩阵乘法
  • -1, 计算逆矩阵

手动计算过程:


# 因变量矩阵
> y1<-matrix(y,ncol=1)  

# 自变量矩阵,增加截距列
> x1<-matrix(c(x,rep(1,length(x))),ncol=2)

# 求解最小二乘解
> solve(t(x1)%*%x1) %*% t(x1) %*% y1
          [,1]
[1,] 0.9004584
[2,] 0.8310557

对应到一元线性回归方程中,y = 0.9004584 * x + 0.8310557,与lsfit()的结果一致。

3. 线性最小二乘回归

接下来,我们先进行线性回归分析,以上面的数据集作为样本,只包含x和y。以普通最小二乘法进行线性回归,计算出预测值并和观测值进行比较,计算预测值和观测值的误差平方和,使得误差函数最小。线性回归模型的详细介绍,可以参考文章R语言解读一元线性回归模型, R语言解读多元线性回归模型

建立线性回归模型


# 建立回归模型
> line<-lm(y~x)

# 查看回归模型
> summary(line)
Call:
lm(formula = y ~ x)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.1549 -0.9550 -0.3472  0.9116  1.6870 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.8311     0.9440   0.880 0.404338    
x             0.9005     0.1698   5.302 0.000726 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.17 on 8 degrees of freedom
Multiple R-squared:  0.7785,	Adjusted R-squared:  0.7508 
F-statistic: 28.12 on 1 and 8 DF,  p-value: 0.0007263

# 残差平方和
> sum(line$residuals^2)
[1] 10.95334

通过查看模型的结果数据,我们可以发现通过T检验自变量x是非常显著,但截距Intercept不显著。通过F检验判断出整个模型的自变量方差是非常显著。通过R^2的相关系数检验可以判断自变量和因变量是不是高度相关的,残差平方和(residual sum-of-squares)为10.95334,效果不好。

接下来,使用模型进行预测,生成预测值。


# 设置x的取值区间
> new <- data.frame(x = seq(min(x),max(x),len = 100))

# 预测结果
> pred<-predict(line,newdata=new)

进行可视化,把观测值与预测值进行展示,以散点表现观测值,以线表示预测值。


> plot(x,y)
> lines(new$x,pred,col="red",lwd=2)

从图中我们可以很容易看出,拟合的效果并不好,与R^2检验不通过是一致的。那么,如果线性模型表现不好,我们可以试试非线性的模型,还是用最小二乘法进行非线性的回归。

4. 非线性最小二乘回归

由于散点的分布并不是线性的,我们再尝试用非线性的函数,进行最小二乘回归的拟合,包括一元二次函数,一元三次函数,指数函数,倒数函数。在R语言中,我们可以使用nls()函数,进行非线性最小二乘回归。

4.1 一元二次函数非线性拟合
定义公式 y ~ a*x+b*x^2+c ,其中x, x^2是一次函数和二次函数,a,b,c分别是参数,我们需要给定参数的初始值。


# 建立回归模型
> nline1 <- nls(y ~ a*x+b*x^2+c, start = list(a = 1,b = 1,c=2))

# 查看模型结果
> nline1
Nonlinear regression model
  model: y ~ a * x + b * x^2 + c
   data: parent.frame()
       a        b        c 
0.008987 0.082188 2.850389 
 residual sum-of-squares: 9.758

Number of iterations to convergence: 1 
Achieved convergence tolerance: 3.024e-06

# 分析模型
> summary(nline1)

Formula: y ~ a * x + b * x^2 + c

Parameters:
  Estimate Std. Error t value Pr(>|t|)
a 0.008987   0.977890   0.009    0.993
b 0.082188   0.088760   0.926    0.385
c 2.850389   2.379763   1.198    0.270

Residual standard error: 1.181 on 7 degrees of freedom

Number of iterations to convergence: 1 
Achieved convergence tolerance: 3.024e-06

通过查看模型的结果数据,我们可以发现通过T检验相关系数a,b,c都不显著,残差平方和(residual sum-of-squares)为9.758,数值过大,效果不好。

把观测值与预测值进行可视化展示,以散点表现观测值,以线表示预测值。


# 模型预测
> npred1<-predict(nline1,newdata=new)

# 可视化
> plot(x,y)
> lines(new$x,npred1,col="blue",lwd=2)

从图中我们可以看出,拟合的是一条曲线,效果不太好。接下来,我们再试试一元三次函数。

4.2 一元三次函数非线性拟合

定义公式 y ~ a*x+b*x^2+c*x^3+d ,其中x, x^2, x^3是一次函数,二次函数和三次函数,a,b,c,d分别是参数,我们需要给定参数的初始值。

使用一元三次函数,进行非线性拟合。


> nline2 <- nls(y ~ a*x+b*x^2+c*x^3+d, start = list(a = 1,b = 1,c=2,d=1))

# 查看模型结果
> nline2
Nonlinear regression model
  model: y ~ a * x + b * x^2 + c * x^3 + d
   data: parent.frame()
      a       b       c       d 
 10.011  -1.822   0.110 -12.516 
 residual sum-of-squares: 3.453

Number of iterations to convergence: 2 
Achieved convergence tolerance: 3.983e-08

# 分析模型
> summary(nline2)

Formula: y ~ a * x + b * x^2 + c * x^3 + d

Parameters:
   Estimate Std. Error t value Pr(>|t|)  
a  10.01051    3.08630   3.244   0.0176 *
b  -1.82152    0.57797  -3.152   0.0198 *
c   0.10998    0.03323   3.310   0.0162 *
d -12.51586    4.88778  -2.561   0.0429 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7586 on 6 degrees of freedom

Number of iterations to convergence: 2 
Achieved convergence tolerance: 3.983e-08

通过查看模型的结果数据,我们可以发现通过T检验相关系数a,b,c,d在0.01的置信区间显著,残差平方和(residual sum-of-squares)为3.453,数值较小,效果应该算是还可以的,能够表达拟合的特征。

把观测值与预测值进行可视化展示,以散点表现观测值,以线表示预测值。


> npred2<-predict(nline2,newdata=new)
> plot(x,y)
> lines(new$x,npred2,col="green",lwd=2)

从图中我们可以看出,拟合的曲线基本穿过观测值的区域,直观上说明效果还不错,与数据检验给出的结论是一致的。

4.3 指数函数非线性拟合

定义公式 y ~ a*exp(b*x)+c ,其中exp(x)是指数函数,a,b,c分别是参数,我们需要给定参数的初始值。


> nline3 <- nls(y ~ a*exp(b*x)+c, start = list(a = 1,b = 1,c=2))

> nline3
Nonlinear regression model
  model: y ~ a * exp(b * x) + c
   data: parent.frame()
     a      b      c 
0.6720 0.2718 2.2087 
 residual sum-of-squares: 8.966

Number of iterations to convergence: 17 
Achieved convergence tolerance: 7.136e-06

> summary(nline3)

Formula: y ~ a * exp(b * x) + c

Parameters:
  Estimate Std. Error t value Pr(>|t|)
a   0.6720     1.2127   0.554    0.597
b   0.2718     0.1764   1.541    0.167
c   2.2087     2.2447   0.984    0.358

Residual standard error: 1.132 on 7 degrees of freedom

Number of iterations to convergence: 17 
Achieved convergence tolerance: 7.136e-06

通过查看模型的结果数据,我们可以发现通过T检验相关系数a,b,c都不显著,残差平方和(residual sum-of-squares)为8.966,数值过大,效果不好。


> npred3<-predict(nline3,newdata=new)
> plot(x,y)
> lines(new$x,npred3,col="yellow",lwd=2)

从图中我们可以看出,指数函数的拟合曲线,和二次函数表现差不多,效果不太好。

4.4 倒数函数非线性拟合

定义公式 y ~ y ~ b/x+c ,其中1/x是倒数函数,b,c分别是参数,我们需要给定参数的初始值。

使用倒数函数,进行拟合。

> nline4 <- nls(y ~ b/x+c, start = list(b = 1,c=2))

# 查看模型结果
> nline4
Nonlinear regression model
  model: y ~ b/x + c
   data: parent.frame()
    b     c 
-17.0   9.5 
 residual sum-of-squares: 15.74

Number of iterations to convergence: 1 
Achieved convergence tolerance: 1.259e-07

> summary(nline4)

Formula: y ~ b/x + c

Parameters:
  Estimate Std. Error t value Pr(>|t|)    
b  -17.003      4.108  -4.139  0.00326 ** 
c    9.500      1.077   8.817 2.15e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.403 on 8 degrees of freedom

Number of iterations to convergence: 1 
Achieved convergence tolerance: 1.259e-07

通过查看模型的结果数据,我们可以发现通过T检验相关系数b,c非常显著,残差平方和(residual sum-of-squares)为15.74,数值较大,效果不好。

把观测值与预测值进行可视化展示,以散点表现观测值,以线表示预测值。


> npred4<-predict(nline4,newdata=new)
> plot(x,y)
> lines(new$x,npred4,col="gray",lwd=2)


从图中我们可以看出,倒数函数的拟合曲线,确实也效果不太好,而且最后一个点有个想离群点,如果没有最后一个点这个效果应该是更好的。

4.5 模型优化
根据上面的几种模型验证的结果,我手动进行了优化,发现高次函数进行拟合效果会更好。定义公式 y ~ a*x^3+b*x^5+c*x^7+e*d^9+e*x^11 ,其中高次函数定义为11次,并去掉了截距。


> nline5 <- nls(y ~ a*x^3+b*x^5+c*x^7+d*x^9+e*x^11,
+               start = list(a=1,b=2,c=1,d=1,e=1))

# 查看模型结果
> nline5
Nonlinear regression model
  model: y ~ a * x^3 + b * x^5 + c * x^7 + d * x^9 + e * x^11
   data: parent.frame()
         a          b          c          d          e 
 2.959e-01 -2.064e-02  5.837e-04 -7.324e-06  3.367e-08 
 residual sum-of-squares: 2.583

Number of iterations to convergence: 3 
Achieved convergence tolerance: 2.343e-07

# 模型分析
> summary(nline5)

Formula: y ~ a * x^3 + b * x^5 + c * x^7 + d * x^9 + e * x^11

Parameters:
    Estimate Std. Error t value Pr(>|t|)   
a  2.959e-01  5.132e-02   5.767   0.0022 **
b -2.064e-02  6.176e-03  -3.341   0.0205 * 
c  5.837e-04  2.468e-04   2.365   0.0644 . 
d -7.324e-06  3.940e-06  -1.859   0.1222   
e  3.367e-08  2.139e-08   1.574   0.1763   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7187 on 5 degrees of freedom

Number of iterations to convergence: 3 
Achieved convergence tolerance: 2.343e-07

通过查看模型的结果数据,我们可以发现通过T检验相关系数a,b,c是显著,残差的误差平方和(residual sum-of-squares)为2.583,数值是最小的,而且比三次函数增长了2上系数,效果应该算不错的,能够表达拟合的特征。

把观测值与预测值进行可视化展示,以散点表现观测值,以线表示预测值。


> npred5<-predict(nline5,newdata=new)
> plot(x,y)
> lines(new$x,npred5,col="purple3",lwd=2)

从图中我们可以看出,拟合的曲线穿过观测值的区域,并且非常贴近部分数据点,直观上说明效果很好,与数据检验给出的结论是一致的。

4.6 总结

我们把上面的1种线性函数和5种非线性的拟合函数的模型评价指标,整理到一个表格中。其中,三次函数和高次函数的拟合效果是不错的,而其他的几种拟合效果,都是不太好的。判断好坏的标准,就是残差的误差平方和,同时参数的T检验显著,才是比较好的拟合模型。

拟合函数 残差平方和 参数T检验显著性 解释性 颜色
一次函数 10.95334 不显著 红色
二次函数 9.758 不显著 蓝色
三次函数 3.453 显著 不好 绿色
指数函数 8.966 不显著 不好 黄色
倒数函数 15.74 显著 灰色
高次函数 2.583 显著 无法解释 紫色

最后,合并所有拟合曲线,画到同一张图上,进行对比。绿色的三次函数拟合曲线,和紫色的高次函数拟合曲线效果最好。

但是,高次函数也带来了的负面问题,一是过拟合,二是业务的可解释性非常差。如果x和y分别带入到业务场景,比如x是货物的重量和y是货物的单价,那么x的11次方,完全是无法理解的,等同于一个黑盒的模型。所以,真正落地项目在使用的时候,要根据业务场景的要求,在准确度,可解释性,复杂度等标准之下,进行权衡。

本文利用最小二乘法,用R语言进行了线性回归和非线性回归的实验。验证了最小二乘法在处理回归问题上是非常有效的,除了回归问题,最小二乘法还能解决分类的问题。最小二乘法作为机器学习的底层算法,是非常值得我们学习和掌握的,明白原理就能掌握一类的高级模型。

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

打赏作者

R语言实现聚类kmeans

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

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

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

关于作者:

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

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

前言

聚类属于无监督学习中的一种方法,k-means作为数据挖掘的十大算法之一,是一种最广泛使用的聚类算法。我们使用聚类算法将数据集的点,分到特定的组中,同一组的数据点具有相似的特征,而不同类中的数据点特征差异很大。PAM是对k-means的一种改进算法,能降低异常值对于聚类效果的影响。

聚类可以帮助我们认识未知的数据,发现新的规律。

目录

  1. k-means实现
  2. PAM实现
  3. 可视化和段剖面图

1. k-means实现

k-means算法,是一种最广泛使用的聚类算法。k-means以k作为参数,把数据分为k个组,通过迭代计算过程,将各个分组内的所有数据样本的均值作为该类的中心点,使得组内数据具有较高的相似度,而组间的相似度最低。

k-means工作原理:

  1. 初始化数据,选择k个对象作为中心点。
  2. 遍历整个数据集,计算每个点与每个中心点的距离,将它分配给距离中心最近的组。
  3. 重新计算每个组的平均值,作为新的聚类中心。
  4. 上面2-3步,过程不断重复,直到函数收敛,不再新的分组情况出现。

k-means聚类,适用于连续型数据集。在计算数据样本之间的距离时,通常使用欧式距离作为相似性度量。k-means支持多种距离计算,还包括maximum, manhattan, pearson, correlation, spearman, kendall等。各种的距离算法的介绍,请参考文章R语言实现46种距离算法

1.1 kmeans()函数实现

在R语言中,我们可以直接调用系统中自带的kmeans()函数,就可以实现k-means的聚类。同时,有很多第三方算法包也提供了k-means的计算函数。当我们需要使用kmeans算法,可以使用第三方扩展的包,比如flexclust, amap等包。

本文的系统环境为:

  • Win10 64bit
  • R: 3.4.4 x86_64-w64-mingw32

接下来,让我们做一个k-means聚类的例子。首先,创建数据集。

# 创建数据集
> set.seed(0)
> df <- rbind(matrix(rnorm(100, 0.5, 4.5), ncol = 2),
+             matrix(rnorm(100, 0.5, 0.1), ncol = 2))
> colnames(df) <- c("x", "y")
> head(df)
              x          y
[1,]  6.1832943  1.6976181
[2,] -0.9680501 -1.1951622
[3,]  6.4840967 11.4861408
[4,]  6.2259319 -3.0790260
[5,]  2.3658865  0.2530514
[6,] -6.4297752  1.6256360

使用stats::kmeans()函数,进行聚类。


> cl <- kmeans(df,2); cl
K-means clustering with 2 clusters of sizes 14, 86

Cluster means:                   # 中心点坐标
          x         y
1  5.821526 2.7343127
2 -0.315946 0.1992429

Clustering vector:               # 分组的索引
  [1] 1 2 1 1 2 2 2 2 2 1 1 2 2 2 2 2 2 2 2 2 2 2 2 1 2 1 1 2 2 2 2 2 2 2 2 1 1 2 1 2 1 2 2 2 2 2 2 1 1 2
 [51] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2

Within cluster sum of squares by cluster:   
[1] 316.0216 716.4009                       # withinss,分组内平方和  
 (between_SS / total_SS =  34.0 %)          # 组间的平方和/总平方和,用于衡量点聚集程度

Available components:            # 对象属性
[1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss" "betweenss"   
[7] "size"         "iter"         "ifault"      

# 查看数据分组情况,第1组86个,第2组14个
> cl$size
[1] 86 14

对象属性解读:

  • cluster,每个点的分组
  • centers,聚类的中心点坐标
  • totss,总平方和
  • withinss,每个分组内的平方和
  • tot.withinss,分组总和,sum(withinss)
  • betweenss,组间的平方和,totss – tot.withinss
  • size,每个组中的数据点数量
  • iter,迭代次数。
  • ifault,可能有问题的指标

1.2 kcca()函数实现
我们再使用flexclust::kcca()函数,进行聚类。


# 安装flexclust包
> # install.packages("flexclust")
> library(flexclust)

# 进行聚类
> clk<-kcca(df,k=2);clk
kcca object of family ‘kmeans’ 

call:
kcca(x = df, k = 2)

cluster sizes:  # 聚类的分组大小
 1  2 
84 16 

# 聚类的中心
> clk@centers
              x         y
[1,] -0.3976465 0.2015319
[2,]  5.4832702 2.4054118

# 查看聚类的概览信息
> summary(clk)
kcca object of family ‘kmeans’ 

call:
kcca(x = df, k = 2)

cluster info:         # 每个组的基本信息,包括分组数量,平均距离、最大距离、分割值
  size  av_dist max_dist separation
1   84 2.102458 9.748136   3.368939
2   16 3.972920 9.576635   3.189891

convergence after 5 iterations                   # 5次迭代
sum of within cluster distances: 240.1732        # 聚类距离之和

我们比较2个不同包的k-means算法,所得到的分组数据都是一样的,中心点位置略有一点偏差。接下来,我们可以把聚类画图。

> plot(df, col = cl$cluster, main="Kmeans Cluster")
> points(cl$centers, col = 1:3, pch = 10, cex = 4) # 画出kmeans()函数效果

从上图中看到k-means的总分2组,每个组的中心点分别用红色十字圆圈和黑色十字圆圈表示,为组内的所有数据样本的均值。再叠加上kcca()函数聚类后的中心点画图。

> points(clk@centers, col = 3:4, pch = 10, cex = 4)  # 画出kcca()函数效果


新的中心点,分别用别用绿色十字圆圈和蓝色十字圆圈表示。虽然我们使用了相同的算法,分组个数也相同,但中心点还有一些不同的。

这里其实就要对聚类的稳定性进行判断了,有可能是聚类迭代次数过少,就会出现不同的聚类结果,就需要增加迭代次数,达到每次计算结果是一致的。也有可能是因为不同的包,实现的代码有所区别导致的。

k-means算法,也有一些缺点就是对于孤立点是敏感的,会被一些极端值影响聚类的效果。一种改进的算法是PAM,用于解决这个问题。PAM不使用分组平均值作为计算的参照点,而是直接使用每个组内最中心的对象作为中心点。

2. PAM实现

PAM(Partitioning Around Medoids),又叫k-medoids,它可以将数据分组为k个组,k为数量是要事前定义的。PAM与k-means一样,找到距离中心点最小点组成同一类。PAM对噪声和异常值更具鲁棒性,该算法的目标是最小化对象与其最接近的所选对象的平均差异。PAM可以支持混合的数据类型,不仅限于连续变量。

PAM算法分为两个阶段:

  1. 第1阶段BUILD,为初始集合S选择k个对象的集合。
  2. 第2阶段SWAP,尝试用未选择的对象,交换选定的中心点,来提高聚类的质量。

PAM的工作原理:

  1. 初始化数据集,选择k个对象作为中心。
  2. 遍历数据点,把每个数据点关联到最近中心点m。
  3. 随机选择一个非中心对象,与中心对象交换,计算交换后的距离成本
  4. 如果总成本增加,则撤销交换的动作。
  5. 上面2-4步,过程不断重复,直到函数收敛,中心不再改变为止。

优点与缺点:

  • 消除了k-means算法对于孤立点的敏感性。
  • 比k-means的计算的复杂度要高。
  • 与k-means一样,必须设置k的值。
  • 对小的数据集非常有效,对大数据集效率不高。

在R语言中,我们可以通过cluster包来使用pam算法函数。cluster包的安装很简单,一条命令就安装完了。


> install.packages("cluster")
> library(cluster)

pam()函数定义:


pam(x, k, diss = inherits(x, "dist"), metric = "euclidean",
    medoids = NULL, stand = FALSE, cluster.only = FALSE,
    do.swap = TRUE,
    keep.diss = !diss && !cluster.only && n < 100,
    keep.data = !diss && !cluster.only,
    pamonce = FALSE, trace.lev = 0)

参数列表:

  • x,数据框或矩阵,允许有空值(NA)
  • k,设置分组数量
  • diss,为TRUE时,x为距离矩阵;为FALSE时,x是变量矩阵。默认为FALSE
  • metric,设置距离算法,默认为euclidean,距离矩阵忽略此项
  • medoids,指定初始的中心,默认为不指定。
  • stand,为TRUE时进行标准化,距离矩阵忽略此项。
  • cluster.only,为TRUE时,仅计算聚类结果,默认为FALSE
  • do.swap,是否进行中心点交换,默认为TRUE;对于超大的数据集,可以不进行交换。
  • keep.diss,是否保存距离矩阵数据
  • keep.data,是否保存原始数据
  • pamonce,一种加速算法,接受值为TRUE,FALSE,0,1,2
  • trace.lev,日志打印,默认为0,不打印

我们使用上面已创建好的数据集df,进行pam聚类,设置k=2。

> kclus <- pam(df,2)

# 查看kclus对象
> kclus
Medoids:                                     # 中心点
     ID         x         y
[1,] 27 5.3859621 1.1469717
[2,] 89 0.4130217 0.4798659

Clustering vector:                           # 分组
  [1] 1 2 1 1 2 2 2 2 2 1 1 2 2 2 2 2 2 2 2 2 2 2 2 1 2 1 1 2 2 2 2 2 2 2 1 1 1 2 1 2 1 2 2 2 2 2 2 1 1 2
 [51] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2

Objective function:                          # 目标函数的局部最小值
   build     swap                           
2.126918 2.124185 

Available components:                        # 聚类对象的属性
 [1] "medoids"    "id.med"     "clustering" "objective"  "isolation"  "clusinfo"   "silinfo"   
 [8] "diss"       "call"       "data"      

> kclus$clusinfo        # 聚类的分组数量,每个组的平均距离、最大距离、分割值
     size  max_diss  av_diss diameter separation
[1,]   15 10.397323 4.033095 17.35984   1.556862
[2,]   85  9.987604 1.787318 15.83646   1.556862

属性解读:

  • medoids,中心点的数据值
  • id.med,中心点的索引
  • clustering,每个点的分组
  • objective,目标函数的局部最小值
  • isolation,孤立的聚类(用L或L*表示)
  • clusinfo,每个组的基本信息
  • silinfo,存储各观测所属的类、其邻居类以及轮宽(silhouette)值
  • diss,不相似度
  • call,执行函数和参数
  • data,原始数据集

把聚类画图输出。

# 画图
> plot(df, col = kclus$clustering, main="Kmedoids Cluster")
> points(kclus$medoids, col = 1:3, pch = 10, cex = 4)

图中,PAM聚类后分为2组,红色一组,黑色一组,用十字圆圈表示2个中心点,可以清晰地看到中心点就是数据点。

我们可以在开始计算时,设置聚类的中心点,为索引1,2坐标点,打印聚类的日志,查看计算过程。


# 设置聚类的中心为1,2
> kclus2<-pam(df,2,medoids=c(1,2),trace.lev=20)
C pam(): computing 4951 dissimilarities from  100 x 2  matrix: [Ok]
pam()'s bswap(*, s=21.837, pamonce=0): medoids given
  after build: medoids are   1   2
  and min.dist dysma[1:n] are
      0      0   9.79   4.78   3.63   6.15   5.23  0.929   8.44   8.59
   2.29   2.69   4.48   1.19   1.98   2.81   5.39    4.2   3.72   4.56
   1.84   3.99    2.4    2.7   4.84   5.08  0.969   2.01   4.94   5.06
   1.94    7.4   5.19   1.62   3.94   3.12   3.51   0.65   4.46   4.61
   5.16   4.57   1.82   3.21   5.79   4.01   5.59   5.38   1.95    6.2
   2.41   2.09    2.2   2.43   2.24   2.26   2.09   2.39   2.21   2.33
   2.24   2.14   2.45   2.37    2.2   2.37   2.13   2.33   2.25   2.18
   2.38   2.19   2.15   2.14    2.1   2.39   2.24   2.24   2.12   2.14
   2.34   2.18   2.25   2.26   2.33   2.17   2.18   2.12   2.17   2.27
   2.29   2.26   2.38   2.12   2.25   2.33   2.09   2.21   2.24   2.13
   swp new  89 <->   2 old; decreasing diss. 306.742 by -93.214
   swp new  27 <->   1 old; decreasing diss. 213.528 by -1.10916
end{bswap()}, end{cstat()}

# 查看中心
> kclus2$id.med
[1] 27 89

通过日志查看,我们可以清楚地看到,2个中心的选择过程,分别用89替换1,距离成本减少93.214,用27替换2,距离成本减少1.1。

PAM作为k-means的一种改进算法,到底结果是否更合理,还要看最终哪种结果能够准确地表达业务的含义,被业务人员所认可,就需要不断地和业务人员来沟通。

3. 可视化和段剖面图

我们实现了聚类计算后,通常需要把复杂的数据逻辑,用简单的语言和图形来解释给业务人员,聚类的可视化就很重要的。如果数据量不太大,参与聚类的指标维度不太多的时候,我们可以用2维散点图,把指标两两画出来。

我们对iris数据集,进行k-means聚类分成3组,画出聚类后的2维散点图结果。

> res <- kmeans(iris[,1:4], centers=3)
> pairs(iris, col = res$cluster + 1)


每2个维度就会生成一张图, 我们可以全面直观的看到聚类的效果。

高级画图工具,使用GGally包中的ggpairs()函数。

> library(GGally)
> ggpairs(iris,columns = 1:5,mapping=aes(colour=as.character(res$cluster)))


图更漂亮了而且包含更多的信息,除了2维散点图,还包括了相关性检查,分布图,分箱图,频率图等。用这样的可视化效果图与业务人员沟通,一定会非常愉快的。

但是如果数据维度,不止3个而是30个,数据量也不是几百个点,而是几百万个点,再用2维散点图画出来就会很难看了,而且也表达不清,还会失去重点,计算的复杂度也是非常的高。

当数据量和数据维度多起来,我们就需要用段剖面图来做展示了,放弃个体特征,反应的群体特征和规律。

使用flexclust包中的barchart()函数,画出段剖面图,我们还是用iris数据集进行举例。


> library(flexclust)
> clk2 <- cclust(iris[,-5], k=3);clk2
kcca object of family ‘kmeans’ 

call:
cclust(x = iris[, -5], k = 3)

cluster sizes:
 1  2  3 
39 61 50 

# 画出段剖面图
> barchart(clk2,legend=TRUE)

如上图所示,每一区块是一个类别,每行是不同的指标。红点表示均值,柱状是这个类别每个指标的情况,透明色表示不重要指标。

查看段剖面图,可以清楚的看到,每个分组中特征是非常明显的。

  • Cluster1中,有39个数据点占26%,Sepal.Width指标在均值附近,其他指标都大于均值。
  • Cluster2中,有61个数据点占41%,Sepal.Width指标略小于均值,其他指标在均值附近。
  • Cluster3中,有50个数据点占33%,Sepal.Width略大于均值,其他指标都小于均值。

从段剖面图,我们可以一眼就能直观地发现数据聚类后的每个分组的总体特征,而不是每个分组中数据的个体特征,对于数据的解读是非常有帮助的。

对于段剖面图,原来我并不知道是什么效果。在和业务人员沟通中,发现他们使用SAS软件做出了很漂亮的段剖面图,而且他们都能理解,后来我发现R语言也有这个工具函数,图确实能极大地帮助进行数据解读,所以写了这篇文章记录一下。

本文介绍了k-means的聚类计算方法和具体的使用方法,也是对最近做了一个聚类模型的总结。作为数据分析师,我们不仅自己能发现数据的规律,还要让业务人员看明白你的思路,看懂数据的价值,这也是算法本身的价值。

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

打赏作者

R语言轻巧的时间包hms

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

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

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

关于作者:

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

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

前言
时间是数据的基本维度,是在做数据处理的时候,必须要掌握的技术。根据时间周期的不同,通常把时间分为,年、月、日、时、分、秒、毫秒等。对于年月日的数据是最常见的,也有很多的处理工具,时分秒的数据通常也会用处理日期的工具,这样有时候就不太方便。

hms包,很小很轻,专注于时、分、秒的时间数据处理。

目录

  1. hms包介绍
  2. hms包的使用

1. hms包介绍

hms包,用于存储和格式化时间,基于difftime类型,使用S3的面向对象数据结构。

本文的系统环境为:

  • Win10 64bit
  • R: 3.4.2 x86_64-w64-mingw32

安装hms包,非常简单,一条命令就可以了。


~ R
> install.packages("hms")
> library(hms)

函数列表:

  • hms: 创建一个hms类型对象
  • is.hms: 判断是否是hms类型
  • parse_hm: 解析hm值
  • parse_hms: 解析hms值
  • round_hms:四舍五入对齐
  • trunc_hms:裁剪对齐
  • as.hms: hms泛型函数,S3类型,用于隐式调用它的继承函数
  • as.hms.character: character转型hms,用于as.hms的继承调用
  • as.hms.default: hms转型hms,用于as.hms的继承调用
  • as.hms.difftime:difftime转型hms,用于as.hms的继承调用
  • as.hms.numeric: numeric转型hms,用于as.hms的继承调用
  • as.hms.POSIXlt: POSIXlt转型hms,用于as.hms的继承调用
  • as.hms.POSIXt: POSIXt转型hms,用于as.hms的继承调用
  • as.character.hms: hms转型character,用于as.character的继承调用
  • as.data.frame.hms: hms转型data.frame,用于as.data.frame的继承调用
  • as.POSIXct.hms: hms转型POSIXct,用于as.POSIXct的继承调用
  • as.POSIXlt.hms: hms转型POSIXlt,用于as.POSIXlt的继承调用
  • format.hms: 格式化hms,用于format的继承调用
  • print.hms: 打印hms对像,用于print的继承调用

从函数列表可以看到,hms包的功能很单一,就在做数据类型和数据变型,是底层的数据结构包,设计思路与zoo包的设计思路一致。zoo包的详细介绍,请参考文章R语言时间序列基础库zoo

hms包中,有大量的as.xxx()函数、format.hms()函数和print.hms()函数,是被用于S3类型函数继承调用的,是不需要我们在写程序的时候显示调用的。S3数据结构详细介绍,请参考文章R语言基于S3的面向对象编程

2. hms包的使用

接下来,我们找几个重点函数进行介绍。

2.1 hms()函数
hms()函数,用于创建一个hms类型的对象。
函数定义:

hms(seconds = NULL, minutes = NULL, hours = NULL, days = NULL)

hms()函数,接收4个参数,分别对应秒,分,时,日。

创建hms对象


# 创建12:34:56的时间对象
> a1<-hms(56, 34, 12);a1
12:34:56

# 创建 10日12:34:56的时间对象
> a2<-hms(56, 34, 12,10);a2
252:34:56

打印结果的第一位252=10*24+12。

2.2 is.hms: 判断是否是hms类型


# 判断是否是hms类型
> is.hms(a1)
[1] TRUE

# 查看hms类型,父类是difftime类型
> class(a1)
[1] "hms"      "difftime"

# 查看hms的属性
> attributes(a1)
$units
[1] "secs"

$class
[1] "hms"      "difftime"

# 查看hms对象的静态数据结构
> str(a1)
Classes 'hms', 'difftime'  atomic [1:1] 45296
  ..- attr(*, "units")= chr "secs"

# 查看面向对象类型
> library(pryr)
> otype(a1)
[1] "S3"

2.3 as.xxx.hms:把hms转型到其他类型


# 默认转型
> as.hms(a1)
12:34:56

# hms转型character,实际会隐式调用as.character.hms()函数
> as.character(a1)
[1] "12:34:56"

# hms转型POSIXct
> as.POSIXct(a1)
[1] "1970-01-01 12:34:56 UTC"

# hms转型POSIXlt
> as.POSIXlt(a1)
[1] "1970-01-01 12:34:56 UTC"

由于我们没有定义as.Date.hms()函数,所以as.Date()函数,不能认识hms类型的转换。


# hms转型Date
> as.Date(a1)
Error in as.Date.default(a1) : 不知如何将'a1'转换成“Date”类别
Error during wrapup: cannot open the connection

自己定义一个as.Date.hms()函数,仅用于转型实验,没有实际业务意义。


# 函数定义
> as.Date.hms<-function(hms){
+   s<-paste(Sys.Date(),' ',hms,sep="")
+   as.Date(s)
+ }

# 显示调用函数
> as.Date.hms(a1)
[1] "2018-12-14"

# 隐式调用函数
> as.Date(a1)
[1] "2018-12-14"

2.4 as.hms.xxx:把其他类型转型到hms


# 把字符串转hms
> as.hms('19:13:14')
19:13:14
# 非法时间字符串转型
> as.hms('19:78:14')
NA

# 数字转型
> as.hms(111312)
30:55:12

# 时间转型
> as.hms(Sys.time())
14:22:59.462795

# 日期转型,同样发生了错误
> as.hms(Sys.Date())
Error: Can't convert object of class Date to hms.
Error during wrapup: cannot open the connection

2.5 parse_hms()/parse_hm()字符串解析

parse_hms对字符串进行转型,对比parse_hms()与as.hms()结果一样的。


# 执行parse_hms
> parse_hms("12:34:56.789")
12:34:56.789
> as.hms("12:34:56.789")
12:34:56.789

# 执行parse_hm
> parse_hm("12:34")
12:34:00
> as.hms("12:34")
NA

打印parse_hms 函数名,查看源代码实现。


> parse_hms 
function (x) {
as.hms(as.difftime(as.character(x), format = "%H:%M:%OS",
units = "secs"))
}
>environment: namespace:hms<

parse_hms()函数,实际就是调用了as.hms()函数。

2.6 round_hms/trunc_hms
round_hms()函数,是把时间进行四舍五入对齐。


# 按秒,以5的倍数进行对齐,四舍五入
> round_hms(as.hms("12:34:51"), 5)
12:34:50
> round_hms(as.hms("12:34:54"), 5)
12:34:55
> round_hms(as.hms("12:34:56"), 5)
12:34:55
> round_hms(as.hms("12:34:59"), 5)
12:35:00

# 按秒,以60的倍数对齐
> round_hms(as.hms("12:34:56"), 60)
12:35:00

trunc_hms()函数,是把时间进行裁剪对齐。


# 按秒去掉末位,以5的倍数进行对齐
> trunc_hms(as.hms("12:34:01"), 5)
12:34:00
> trunc_hms(as.hms("12:34:44"), 5)
12:34:40
> trunc_hms(as.hms("12:34:56"), 60)
12:34:00

2.7 在data.frame中插入hms列


# 创建data.frame
> df<-data.frame(hours = 1:3, hms = hms(hours = 1:3))
> df
  hours      hms
1     1 01:00:00
2     2 02:00:00
3     3 03:00:00

# 查看df的静态结构
> str(df)
'data.frame':	3 obs. of  2 variables:
 $ hours: int  1 2 3
 $ hms  :Classes 'hms', 'difftime'  atomic [1:3] 3600 7200 10800
  .. ..- attr(*, "units")= chr "secs"

hms包很轻巧很简单,但却可以快速提高帮助我们处理时分秒数据,这些基础函数库是需要我们完全掌握和熟练运用的。

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

打赏作者

算法,如何改变命运

架构师的信仰系列文章,主要介绍我对系统架构的理解,从我的视角描述各种软件应用系统的架构设计思想和实现思路。

从程序员开始,到架构师一路走来,经历过太多的系统和应用。做过手机游戏,写过编程工具;做过大型Web应用系统,写过公司内部CRM;做过SOA的系统集成,写过基于Hadoop的大数据工具;做过外包,做过电商,做过团购,做过支付,做过SNS,也做过移动SNS。以前只用Java,然后学了PHP,现在用R和Javascript。最后跳出IT圈,进入金融圈,研发量化交易软件。

架构设计就是定义一套完整的程序规范,坚持架构师的信仰,做自己想做的东西。

关于作者:

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

转载请注明出处:
http://blog.fens.me/architect-algorithm/

前言

近年来,随着大数据的飞跃式的发展,已经越来越深地开始影响到我们的生活,社交有腾讯大数据,购物有阿里大数据,搜索有百度大数据,出行有滴滴大数据等等。当数据越来越多地被积累,就需要算法来挖掘出数据的价值。特别是进入到大数据时代,算法显得越来越重要。

让死的数据变得有价值,就是算法的力量。进入到全民大数据的时代后,数据已经不再是门槛儿,最重要的是算法,算法才是真正能够创造生产力的地方。算法工程师的价值也会越来越大,但是你们真的发掘出来你们的价值了吗?

目录

  1. 算法在各个行业的应用
  2. 投身于哪个行业好?
  3. 金融最靠谱

1. 算法在各个行业的应用

大数据的兴起冲击着各行各业,带来机遇也带来挑战,没有数据你就没有核心价值。当有了数据作为基础,你要继续需要思考如何让数据变的有价值。过去的2016年的投资市场很惨淡,唯有人工智能大火了一把。从深度挖掘(Deep Learning)技术在图像识别领域的精确识别,迭代决策树(GBDT)在数据挖掘算法比赛中频繁获奖,到AlphaGo在围棋领域打败在人类选手,百度小度机器人在最强大脑的舞台上挑战人类脑王等等,这些事件都是算法领域的突破。

算法,真的已经应用到了各行各业,在慢慢地改变着人们的生活和习惯,比如说图像识别,自动驾驶,用户行为,金融征信,量化投资等领域,都在发生着变化。

图像识别领域,深度学习算法异军突起,不仅可以进行准确的人脸识别、指纹识别,还可以进行复杂的图像对比。我深刻记得,2016年参加的光谷人工智能大会上,听西安电子科技大学公茂果教授分享的“深度神经网络稀疏特征学习与空时影像变化检测”主题,利用图像识别技术,对比汶川地震前后的卫星照片和光感照片,准确地找到了受到地震影响最严重的区域,即震前和震后地貌发生变化最大的区域,快速地为救援队定位到最需要帮助的地点,解救伤者,投放救援物资。

自动驾驶领域,可以通过识别路面的状况来实现自动驾驶、自动停车。Uber无人驾驶汽车已经在匹兹堡上路测试,自动驾驶汽车配备了各式传感器,包括雷达、激光扫描仪以及高分辨率摄像头,以便绘制周边环境的细节。自动驾驶汽车有望改善人类的生活质量,也可挽救数百万人的性命,为人们提供更多的出行方便。5年前,我在听Andrew Ng的斯坦福大学机器学习公开课的时候,就被当时的自动驾驶视频介绍所震撼,科幻电影中的世界就快变成现实了。

用户行为分析,人类有各种各样的行为和需求。衣食住行,吃喝玩乐,都是人的最基本的行为。大多数人的行为是共性的,商家可以收集这些行为数据,通过数据挖掘算法来找到人们行为共性的规律。根据用户的购物行为,商家可以为用户推荐喜欢的商品,这样就有了推荐系统; 根据用户对信息的查询行为,可以发现用户对信息的需求,这样就有了搜索引擎;根据用户位置的变化,可以发现用户的出行需求,这样就有了地图应用;针对用户个性化的行为,可以给用户打上标签,用来标注用户的特征或身份,这样就有了用户画像。用户行为分析,让商家了解用户习惯,同时也让用户了解自己,有巨大的商业价值。

在金融领域也有很多,算法应用的场景。

金融征信领域,传统信贷业务都是银行核心业务,但由于中国人数众多且小客户居多,银行无法负担为小客户服务的高成本,导致民间信贷的兴起。2014年底互联网金融P2P的开始爆发,贷款需求被满足的同时,却暴露出了违约风险。征信体系缺失,导致很多P2P公司坏账率很高,到2016年底P2P跑路的多达数千家。征信需求,变得非常迫切。比如,某个人想买车但现金不够,这时就需要进行贷款。商家给用户进行贷款时,通过信用风险的评级就能判断出这个用户的还款能力,从而来决定给他贷多少钱,以什么周期还款,减少违约风险。支付宝的芝麻信用分,是目前被市场一致认可的信用评分模型。

量化投资领域,我认为这个领域最复杂的,最有挑战性的,同时也是最有意思的。可以通过量化算法模型实现赚钱,是最容易变现的一种方法。在金融投资领域中,有各式各样的数据,反应的各种金融市场的规则,有宏观数据,经济数据,股票数据,债券数据,期货数据,还有新闻数据,情绪数据等等,金融宽客(Quant)通过分析各种各样的数据,判断出国家的经济形势和个股的走势,进行投资组合算法,实现投资的盈利。

看到这里,我想问问大家,你们脑子里那些聪明的想法,有没有被金融行业的魅力撩出些许的荷尔蒙?

2. 投身于哪个行业好?

从上面各个行业的算法应用来说,都有很广阔的应用前景。作为一个算法的研究者,那我们究竟投身到哪个行业更好呢?

这个其实要从多个方面进行考虑,我们的目标是个人价值最大化。那么,你要选择一个自己能够接触到的、完全竞争的、短流程的渠道,利用你的算法技术和对业务的理解实现变现的过程。

其实,满足个人可变现的渠道其实非常有限,你很难通过一个图像识别的算法,直接面向市场进行收钱,你需要有一个承载的产品,而产品研发的过程是非常漫长的。同样地,自动驾驶算法需要汽车生产场商的实验。用户行为分析算法,需要电子商务平台的以用户购买行为进行验证。

量化投资,可以用个人账号在中国二级投资交易市场,完成交易过程。这种方式没有很多的中间环节,你获得交易所的数据,自己编写算法模型,然后用自己的钱去交易,完全自己把握。只要算法有稳定的收益率,你就可以赚到钱。这种变现方法,其实就是量化投资,从金融的角度才是最靠谱的一种变现方法。

3. 金融最靠谱

作为IT人,我们懂编程,懂算法,只要再了解金融市场的规则,就能去金融市场抢钱了。中国的金融二级投资交易市场,是一个不成熟的市场,同时又是情绪化的市场。市场中,每天都存在着大量的交易机会,每天都会有“乌龙指”。量化投资的技术,可以帮助我们发现这些由于信息不对称出现的机会,赚取超额的收益。

那么到底怎么做量化投资呢?。

下面举个例子,一个私募基金,募集了1亿资金准备杀入金融市场。基金经理决定按照投资组合的建模思路,对各类金融资产进行组合配置。下图就反应了各类资产,以均值-方差的标准来创建投资组合,符合资本资产定价模型(CAPM)的原理。关于资本资产定价模型详细介绍,请参考文章R语言解读资本资产定价模型CAPM

图中,x轴为收益率的标准差,y轴为收益率的均值,图中的点构建了可投资区域,每个点代表一个可投资产品,每条虚线连接的点的集合,就是一个有效的投资组合。

对于,图中近百个点来说,假设每次要配置5种资产做投资组合,那么就是75287520种组合方法;如果配置10种资产,可选方案就是一个非常大的数字了。

我可以用R语言来计算一下,投资组合的数量。


# 100个选5个,做组合
> choose(100,5) 
[1] 75287520

# 100个选10个,做组合
> choose(100,10) 
[1] 1.731031e+13

对于金融市场来说,有非常多的金融资产可供我们来选择。中国A股股票有3000多只,基金2000多支, 债券3000多支,期货100多支,还有大综商品,货币市场产品,汇率产品,海外投资市场等。如果把这个多种的资产进行组合,将有无限多的投资组合可以进行选择,是一个无限大的计算量。我们需要利用算法进行组合优化,从而找到市场上最优的投资组合。算法本身,才是最能体现价值的部分。

那么传统的基金是如何进行投资组合的?大多都是靠投资经理的主观投资经验来完成的。在金融市场里,每支基金都配置了不同的资产做组合,我们随便找支基金看看,它的投资组合是如何配置的。比如,华夏成长(000001.OF)基金,它是股债混合型的。数据来源于万得, 2017年2月8日的数据。

从业绩表现来看,这支基金最辉煌的时代在2006-2007年,连续6个月回报101.49%,那么最低1年表现就比较差,为仅落后于沪深300指数,整体排名也都在后面。今年以来收益率0.58%,同类排名144/507;1年收益率-1.45%,同类排名400/487;3年收益11.67%,同类排名378/426;5年收益39.96%,同类排名290/352。

我们再来看一下,这支基金的组合成分,主要是股票和债券。

债券占比 :

证券名称 占净值比 近3月涨跌
12石化01 2.34%↑ -0.49%
116国泰君安CP008 2.12%↑ -0.03%
116农发01 1.91%↑ -0.08%
110营口港 1.70%↑ -1.59%
109常高新 1.62%↑ -0.65%

股票占比:

证券名称 占净值比 近3月涨跌
中工国际 4.09%↑ -0.95%
中国医药 3.85%↑ 0.34%
神雾环保 3.81%↑ 2.56%
东方网络 2.89%↑ -13.00%
立讯精密 1.52%↑ -1.82%
高能环境 1.42%↑ -14.96%
上汽集团 1.38%↑ 7.88%
田中精机 1.31%↑ -12.28%
上海医药 1.25%↑ 5.39%
中牧股份 1.21%↑ -4.25%

从市场上几千支的股票和债券中进行选择,并配置不同的权重,之前都是基金经理干的活,那么我们用算法一样也可以干,说不定用算法模型构建的组合业绩会更好。如果我们用算法模型,取代了年薪几百万的基金经理,那么你就能够获得这个收益。最终实现个人价值,从而用算法改变命运。所以,通过金融变现才是最靠谱的。

转载请注明出处:
http://blog.fens.me/architect-algorithm/

打赏作者

R语言解读资本资产定价模型CAPM

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

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

关于作者:

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

转载请注明出处:
http://blog.fens.me/finance-capm

前言

伴随2016年中国金融交易市场的跌宕起伏,风险越来越不确定,利率持续走低,理财等无风险资产收益持续下降的情况,唯有投资组合才能让我们的资产保值、增值。根据资本资产定价模型(CAPM),通过对金融数据的分析,构建投资组合,帮助我们在有效的市场中控制风险、稳定收益。

本文将深入浅出地介绍资本资产定价模型,从理论到建模,再到程序现实。资本资产定价模型反应的是资产的风险与期望收益之间的关系,风险越高,收益越高。当风险一样时,投资者会选择预期收益最高的资产;而预期收益一样时,投资者会选择风险最低的资产。

由于本文为非金融教材类文章,所以当出现与教课书不符的描述,请以教课书为准。本文力求用简化的语言,来介绍自资本资产定价模型的知识,同时配合R语言的实现。

目录

  1. 故事背景
  2. 资本市场线
  3. 资本资产定价模型
  4. 用R构建投资组合模型
  5. Beta VS Alpha

1. 故事背景

1952年,马科维茨(Markowitz)提出了投资组合选择理论,他认为最佳投资组合应当是,风险厌恶特征的投资者的无差异曲线和资产的有效边界线的交点。投资者在选择资产时会在收益和风险之间做出平衡:当风险一样时,会选择预期收益最高的资产;而预期收益一样时,会选择风险最低的资产。

01

图1 投资组合选择示意图

到1964年,威廉-夏普(William Sharp),约翰-林特纳(John Lintner)与简-莫森(Jan Mossin)则在马科维茨基础上提出的单指数模型,将市场组合引入均值-方差模型,极大地简化了计算,他们认为获得了市场任意资组合的收益与某个共同因素之间是有线性关系,最终将其发展为资本资产定价模型(Capital Asset Pricing Model, CAPM)。从马科维茨的投资组合选择理论,发展到资本资产定价模型经历了一个漫长的过程。

简单一句话概括,资本资产定价模型的核心思想,资产价格取决于其获得的风险价格补偿。

假设条件

资本资产定价模型,是基于一系列假设条件而成立的。但这些条件,可能并不符合现实的标准,资本资产定价模型也一度遭到质疑。

  • 资产可以无限分割。
  • 不存在交易成本和个人所得税。
  • 可以无限卖空。
  • 存在一种无风险利率,投资者在此利率水平下,可以无限制地贷出和借入任意数额的资金。
  • 投资者是价格接受者,市场是完全竞争的。
  • 投资者是理智的,通过比较资产的期望收益和方差来作出投资决策,在相同预期收益下会选择风险最小的资产。
  • 投资者在相同的投资期限出作出决策,而市场信息是公开免费的,并可以及时获得。
  • 投资者对市场中的经济变量有相同的预期,他们对任意资产的预期收益率、市场风险的看法是一致的。

资本资产定价模型的核心假设是认为市场满足完全、无摩擦和信息完会对称的条件,市场中的投资人都是Markowitz理论中的理性经济人。

2. 资本市场线

由于涉及到金融专业领域,有几个概念是我们应该提前知道的。

  • 风险资产:风险资产是指具有未来收益能力的资产,但收益率不确定且可能招致损失,比如股票、债券等。
  • 无风险资产:没有任何风险或者风险非常小的资产,有确定的收益率,并且不存在违约的风险。
  • 收益率:指从投资开始到投资结束时,所获得的投资回报率。
  • 无风险收益率:无风险资产,所产生的投资回报率。
  • 投资组合:由投资人或金融机构所持有的股票、债券、基金、衍生金融产品等组成的集合,目的在于分散风险。
  • 杠杆交易:就是利用小资金来进行数倍于原始金额的投资,以期望获取相对投资标的物波动的数倍收益率的盈利或亏损。

2.1 风险资产

对于风险资产来说,我们可以用预期收益和风险,通过二维的坐标来进行描述。

对上图的解释:

  • X轴,为风险
  • y轴,为收益率
  • 灰色区域,为金融资产可投资区域
  • 黑色线,为有效投资边界
  • A和B点,为2个风险资产

A和B有相同的x值,表示具有相同的风险。B点在A点上面,表示B的收益率高于A。对于理性的投资者来说,如果只在A点和B点之间做投资选择,那么大家都会投资到B点,而不投资于A点。

2.2 无风险资产

在下图中,我们加入无风险资产,来比较无风险资产和风险资产的关系。

对上图的解释:

  • B点,为1个风险资产,在有效投资边界上
  • C点,为无风险资产,在y轴上
  • X轴,为风险
  • y轴,为收益率
  • 灰色区域,金融资产为可投资区域
  • 黑色线,为有效投资边界

C点为无风险资产,他的位置在图示的y轴上,这时x为0,即风险为0。我们可以把投资,分配到C点或B点上。如果都投到C点,那么我们将获得的是R0部分的无风险收益;如果都投到B点,那么我们需要承担σB的风险,同时获得RB的风险收益。如果我们把资金,一部分投资到B点对应的风险资产上,另一部分投资到C点对应的无风险资产上,那么将构成一个由B和C资产组成的投资组合,而且风险和收益部分,将体现在B和C的连线上。

2.3 最优组合

那么,有没有最优的投资组合呢?收益最大、风险最小。下面就让我们来,发现这个最优的组合M。

对上图的解释:

  • M点,为最优组合的风险资产
  • B点,为1个风险资产,在有效投资边界上
  • C点,为无风险资产,在y轴上
  • X轴,为风险
  • y轴,为收益率
  • 灰色区域,金融资产为可投资区域
  • 黑色线,为有效投资边界

假设有最优的组合,在上图中M点处,当我们把C和M进行连线,使得CM的连线与灰色区域相切。从图上看,CM的连线会比任意的C与可投资区域点的连线斜率都要大,比如C和B的连线。我们取CB的连线的延长线,在CB的延长线上找到,与M具有相同x的点B’,这时M与B’风险相同,M点在B’点的上面,所以M点的收益率大。也就是说,当风险相同的时候,我们都会选择收益率最大的资产。

不论从可投资区域中怎么选取,M点都是斜率最大的点,那么我们可以认为,M点为市场上各资产的最优的投资组合.

对于最优的投资组合,其实不管投资者的收益风险的偏好是什么样子的,只要找到了最优的风险资产组合,再加上无风险的资产,就可以为投资者获得最佳的投资方案了。那么对于理性的投资者,如果发现了最优的组合,他们只会投资于这个组合,这时与收益和风险偏好无关。

M点构建的投资组合,一般是由所有可投资证券产品组成的,每种证券资产构成的比例,为证券的相对市值。无风险资产C,并没有包括在M中,人们都会选择CM的连接线进行投资,来构建最优的投资组合。

在实际的市场交易中,金融资产的价格会发生偏离,因为价格受市场的供需关系所影响,当价格发生偏离后,市场会自动修复会回均衡价格水平。

2.4 资本市场线

对于CM的连线,就是马科维茨提出了投资组合选择理论,风险厌恶特征的投资者的无差异曲线和资产的有效边界线的交点。这条线就叫,资本市场线(Capital Market Line)。

资本市场线是指表明有效组合的期望收益率和标准差之间的一种简单的线性关系。

资本市场线决定了证券的价格。因为资本市场线是证券有效组合条件下的风险与收益的均衡,如果脱离了这一均衡,则就会在资本市场线之外,形成另一种风险与收益的对应关系。

2.5 投资组合构建

资本市场线,就是我们最优的投资组合,当我们发现这个投资组合,所有资金都会投到这个组合上。通过对无风险资产C和风险资产M分配不同的投资权重,我们可以自己配置出自己想要的风险和收益来,同时可以利用金融工具来加杠杆放大风险和收益的范围。

如果我们把投资者分成,风险厌恶型和风险激进型。

对于风险厌恶型,他们对于资金安全有非常高的要求,不追求高收益但求本金安全,这些资金通常都是用来生活的。那么在为这些资金做资产配置方案的时候,可以把一部分资金配置无风险资产上,同时少量资金配置到M点的最优组合上,保证低风险并获得少量收益。

如图中CM1点,如果配置50%的风险资产M和50%的无风险资产C,来实现投资组合。公式如下:

CM1 = 0.5C + 0.5M

对于风险激进型,他们对于资金有非常高的收益要求,本金可以部分或全部损失,这些资金通常都是“闲钱”,就是用来进行投资活动的。那么在为这些资金做资配置方案时,可以全部都投到M上,再激进点,可以通过借钱、融资的方式,增加杠杆,把资金放大进行投资。这种操作风险会随着杠杆的放大剧增,当然同时你也会有更大的收益。

如图中CM2点,落在了CM的延长线上。我们可以配置150%的风险资产M,同时用50%的钱去抵押以无风险资产C的收益率去借钱。公式如下:

CM2 = -0.5C + 1.5M

2.6 风险和收益的关系

上面我们描述风险和收益的关系,主要是从思路上定性介绍,没有进行定量描述,那么究竟风险和收益从数学上怎么进行定义呢。

对上图的解释:

  • M点,为最优组合的风险资产
  • C点,为无风险资产,在y轴上
  • r0,为无风险资产的收益率
  • rM,为M点的收益率
  • x轴,σp为风险资产的收益率的方差
  • y轴,rp为收益率

根据威廉-夏普所引入的均值-方差模型,极大地简化了计算,就是解决了公式计算的问题。用方差来刻画风险,建立收益和风险的一元线性关系。可以用下面公式来表示:

公式

E(rm) – r0 = A * σM^2

公式解释:

  • E(rm):市场投资组合的预期收益率
  • r0:无风险收益率
  • E(rm)–r0, 市场投资组合的风险溢价
  • σM^2: 市场投资组合方差Var(rM)
  • A:风险厌恶水平

有了公式,我们就明确的知道了,风险和收益的定量关系,并且可以利用数据来进行计算。

3. 资本资产定价模型

对于市场的投资组合,风险溢价和市场投资组合的方差成线性关系。但对于单个资产来说,收益和风险是市场投资组合组成的一分部,受市场共同变化的影响。

3.1 单个资产风险溢价

对于单个资产的风险来说,在资本资产定价模型中,用β来进行表示。β是衡量单个金融资产与市场收益的共同变化程度,通过协方差来计算。单个资产的风险为,当前资产与投资组合收益率的协议差,除以投资组合收益率的方差。

单个资产的风险的计算公式:


βi = Cov(ri, rm) / Var(rm) 
   = Cov(ri, rm) /  σm^2

单个资产的风险溢价的计算公式:


E(ri) – rf = (Cov(ri, rm) / σm^2)*[E(rm) – rf] 
           =  βi  *  [E(rm) – rf]

对公式的解释:

  • E(ri),为风险资产i的预期收益
  • E(rm),为市场投资组合的预期收益
  • rf,为无风险资产收益
  • Cov(ri, rm),为风险资产收益率和市场投资组合收益率的协议差
  • Var(rm),为市场投资组合的收益率的方差

从公式可以看出,单个资产的风险溢价与市场投资组合M的风险溢价成正比,受β影响。

3.2 资本资产定价模型

资本资产定价模型,是现化金融学中的基石理论。在上述假设条件下,可以推到出资本资产定价模型的具体公式。整个和推到过程,就是上面文章介绍的过程,从后人学习的角度看,这个理论比较简单的,仅用到了简单地统计学知识,但是前人却花了很长的时间研究和探索。

判断单个资产的风险时,当β=1时,则说明当前资产与整个市场的趋势是完全保持一致的;当β为2时,代表高风险,其回报的变化将大于市场大盘的变化幅度;当β为0.5时,代表是低风险的资产配置。

3.3 2种风险

在资本资产定价模型,定义了2种风险,即系统性风险和非系统性风险。

系统性风险,就是由外部因素引起的风险,比如:通货膨胀,GDP,重大政治事件等等。这一类事件对于资产收益率的影响不能通过组合本身来消除的,所以这一类风险对于投资者来说是无法回避的。

非系统性风险,就是组合内部结构引起的风险,比如:A股与B股高度相关,A股的收益率出现大幅波动的时候,B股也会出现相似幅度的波动,波峰叠加或波谷叠加,就会增加整个组合的风险;反之,如果A与B为负相关,则A与B的波动就会相互抵消。这样,风险是由组合里的资产类型决定的,所以通过多样化分散的投资策略,无论在理论还是实际上,这种风险都是可以最小化甚至消除的。而这个消除的过程中,整个投资组合的收益率是不会下降的。

3.4 2种收益

与风险相对应是收益,我们承受了2种风险的同时,也获得了风险所带来的收益。一部分是与市场完全相关收益部分,即beta(β)收益;另一部分与市场不相关的收益部分,即alpha(α)收益。

  • beta收益,相对容易获得,例如,你看好一个市场,可以持有成本低廉的对应市场的指数基金,等待市场上涨。
  • alpha收益,比较难获得,alpha是体现投资水平的策略收益。

alpha是,投资组合的实际期望收益与预期收益之间的差。计算alpha的公式为:


E(ri) – rf = αi + βi  *  [E(rm) – rf]
αi         = [E(ri) – rf] -  βi * [E(rm) – rf]

alpha是衡量投资人投资水平的,我们举个例来说明。比如:市场收益率为14%,A证券的β=1.2,短期国债利率6%,投资者对这只股票的进行了交易,获得的实际收益为17%,那么我们怎么判断投资人的水平呢?

首先,先求出A证券的预期收益率 = 6% + 1.2*(14-6)% = 15.6%,再用投资者实际收益减去A证券预期收益 17% – 15.6% = 1.4%。最后获得的1.4%就是alpha,表示投资者能力,可以额外获得1.4%的收益。

3.5 资本资产定价模型的应用场景

进行组合投资分散风险:投资者可以按市场组合的构成比例分散持有多种风险资产,使持有的风险资产组合最大限度地接近市场组合,以达到消除非系统风险的目的。

调整收益风险比例:将无风险资产与风险资产市场组合进行再组合,以获得所希望的个性化的风险收益组合。

指数化投资:将资产配置在与某一指数相同的权重的投资方法,通过微调权重或成分,获得比指数更好的alpha。

资产定价:资本资产定价模型可以用来判断有价证券或其他金融资产的市场价格是否处于均衡水平,是否被高估或低估,以便通过套利活动获取超额收益。

基金购买:举一个贴近市场的例子,当我们要购买基金时,也可以用到资本资产定价模型帮我们分析。比如,基金A的期望收益率12%,风险β=1,基金B期望收益率13%,β=1.5。市场期望收益率11%,无风险资产收益率r0 = 5%。 那么哪只基金更值得买?

当你每天打开支付宝,看到里面的各种基金推荐。你就会发现这是一个实际的问题。如果你懂学了本文,按照资本资产定价模型的思路,其实就是求alpha,哪个基金的alpha高,就买哪个。

求alpha,我们就直接套用公式。


αA = 12 – 5 – 1 * [11 - 5] = 1%
αB = 13 – 5 – 1.5* [11 -5 ] = -1%

基金A的alpha为1%,而基金B的alpha为-1%。结论就很明显,基金A的管理人能力很好,超额收益1%;而基金B的管理人,就差一些,盈利低于市场1%。所以,我们会投资基金A,而不会投资基金B。

4. 用R构建投资组合模型

花了大量的篇幅介绍了资本资产定价模型的原理,对于程序实现其实是相当简单地。因为R语言中,已经把资本资产定价模型相关的计算函数都封包好了,我们仅仅是调用就能完成整个的计算过程。

R语言程序实现,我们主要会用到2个包,quantmod和PerformanceAnalytics。对于为什么要用R语言,可以参考文章R语言为量化而生

  • quantmod,用于下载数据。
  • PerformanceAnalytics,用于进行各种评价指标计算。

我们设计一个应用场景,假如我有10万美金想投资于美国的股市,我想获得比标普好(SP500)的投资收益,那么我应该如何购买股票。

首先,我们先想清楚,我的最终的目标是“比标普好的投资收益”。其次,我们基于资本资产定价模型理论基础,从投资组合角度思考投资策略,而不是技术指标的角度。比标普好,那么我们就需要以标普指数做为理想投资组合。然后,我们去市场上选择几个股票,分别计算出收益率,beta,alpha等指标,判断是否符合的预期,反复测试,直到找到合适的股票或股票组合。

本文只是案例介绍,用于说明投资思路和方法,不购成任何的股票推荐。

本文的系统环境

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

从yahoo下载IBM,GE(通用电器),YHOO(Yahoo)的3只股票,从2010年01月01日的日行情数据,同时下载标普指数(SP500)的日行情数据。

下面代码并不完整,但思路已经给出,请大家不要太随意地张嘴要数据和代码,毕竟写一篇文章非常辛苦。如果你想直接用我的代码,请扫文章下面二维码,请作者喝杯咖啡吧。 :_D

执行R语言程序。


# 加载程序包
> library(quantmod) 
> library(PerformanceAnalytics)

# 从yahoo下载3只股票的数据,和SP500的数据
> getSymbols(c('IBM','GE','YHOO','^GSPC'), from = '2010-01-01')

# 打印前6行和后6行数据
> head(GSPC)
              open    high     low   close     volume adjusted
2010-01-04 1116.56 1133.87 1116.56 1132.99 3991400000  1132.99
2010-01-05 1132.66 1136.63 1129.66 1136.52 2491020000  1136.52
2010-01-06 1135.71 1139.19 1133.95 1137.14 4972660000  1137.14
2010-01-07 1136.27 1142.46 1131.32 1141.69 5270680000  1141.69
2010-01-08 1140.52 1145.39 1136.22 1144.98 4389590000  1144.98
2010-01-11 1145.96 1149.74 1142.02 1146.98 4255780000  1146.98

> tail(GSPC)
              open    high     low   close     volume adjusted
2016-12-20 2266.50 2272.56 2266.14 2270.76 3298780000  2270.76
2016-12-21 2270.54 2271.23 2265.15 2265.18 2852230000  2265.18
2016-12-22 2262.93 2263.18 2256.08 2260.96 2876320000  2260.96
2016-12-23 2260.25 2263.79 2258.84 2263.79 2020550000  2263.79
2016-12-27 2266.23 2273.82 2266.15 2268.88 1987080000  2268.88
2016-12-28 2270.23 2271.31 2249.11 2249.92 2392360000  2249.92

# 画出SP500的K线图
> barChart(GSPC)

把4个品种的调整后的价格进行合并。


> # 改列名
> names(IBM)<-c("open","high","low","close","volume","adjusted")
> names(GE)<-c("open","high","low","close","volume","adjusted")
> names(YHOO)<-c("open","high","low","close","volume","adjusted")
> names(GSPC)<-c("open","high","low","close","volume","adjusted")

# 数据合并
> dat=merge(IBM$adjusted,GE$adjusted,YHOO$adjusted,GSPC$adjusted)
> names(dat)<-c('IBM','GE','YHOO','SP500')

# 打印前6行
> head(dat)
                IBM       GE  YHOO   SP500
2010-01-04 112.2859 12.27367 17.10 1132.99
2010-01-05 110.9295 12.33722 17.23 1136.52
2010-01-06 110.2089 12.27367 17.17 1137.14
2010-01-07 109.8274 12.90920 16.70 1141.69
2010-01-08 110.9295 13.18724 16.70 1144.98
2010-01-11 109.7680 13.31435 16.74 1146.98

计算每日收益率,合并收益率到dat_ret


> dat_ret=merge(IBM_ret,GE_ret,YHOO_ret,SP500_ret)
> names(dat_ret)<-c('IBM','GE','YHOO','SP500')
> head(dat_ret)
                    IBM           GE         YHOO        SP500
2010-01-04  0.009681385  0.015111695  0.009445041 0.0147147759
2010-01-05 -0.012079963  0.005177994  0.007602339 0.0031156762
2010-01-06 -0.006496033 -0.005151320 -0.003482298 0.0005455205
2010-01-07 -0.003461515  0.051779935 -0.027373267 0.0040012012
2010-01-08  0.010034759  0.021538462  0.000000000 0.0028817272
2010-01-11 -0.010470080  0.009638554  0.002395150 0.0017467554

定义无风险收益率为4%,计算4个资产的平均年化收益率。


# 无风险收益率
> Rf<-.04/12

# 计算平均年化收益率,平均年化标准差,平均年化Sharpe 
> results<-table.AnnualizedReturns(dat_ret,Rf=Rf)
> results
                               IBM      GE    YHOO   SP500
Annualized Return           0.0345  0.1108  0.1257  0.1055
Annualized Std Dev          0.1918  0.2180  0.3043  0.1555
Annualized Sharpe (Rf=84%) -2.8892 -2.3899 -1.6911 -3.3659

统计指标分析,每个资产有1760个样本点,没有NA值。日最小收益率,YHOO最小为-0.0871。日最大收益率,在GE为0.1080。算数平均,几何平均,方差,标准差都是YHOO最大。


# 计算统计指标
> stats
                      IBM        GE      YHOO     SP500
Observations    1760.0000 1760.0000 1760.0000 1760.0000
NAs                0.0000    0.0000    0.0000    0.0000
Minimum           -0.0828   -0.0654   -0.0871   -0.0666
Quartile 1        -0.0060   -0.0065   -0.0098   -0.0039
Median             0.0002    0.0004    0.0005    0.0005
Arithmetic Mean    0.0002    0.0005    0.0007    0.0004
Geometric Mean     0.0001    0.0004    0.0005    0.0004
Quartile 3         0.0067    0.0077    0.0112    0.0053
Maximum            0.0567    0.1080    0.1034    0.0474
SE Mean            0.0003    0.0003    0.0005    0.0002
LCL Mean (0.95)   -0.0004   -0.0001   -0.0002    0.0000
UCL Mean (0.95)    0.0008    0.0012    0.0015    0.0009
Variance           0.0001    0.0002    0.0004    0.0001
Stdev              0.0121    0.0137    0.0192    0.0098
Skewness          -0.5876    0.3084    0.0959   -0.3514
Kurtosis           4.6634    4.7294    2.9990    4.0151

画出IBM股票,日收益和月收益的图,4个资的累积收益率图,并对4个资产做相关性分析。

IBM股票,每日收益图

IBM股票,每月收益图

4个品种的累积收益率图

从上图中可以看出,红线(GE)和蓝线(SP500)的走势基本稳合,说明GE在从2010开始在跟着美国经济持续发展。绿线(YHOO)从2013初到2015年初大幅拉升,领先于SP500很多,说明这段时期YHOO所处的互联网行业,带来了非常大的市场红利;从2015年到2016年,又下跌很大,大起大落,受市场影响非常敏感。黑线(IBM)大部分时间都处于SP500的下方,说明美国经济这几年的高速发展,并没有给IBM带来很大的发展空间。如果从我们的目标来说,”比标普好的投资收益”那么我们只能选择GE或YHOO。

相关性分析

对4个品种进行相关性分析,发现GE和SP500相关系数为0.78,是3只股票中最相关的。而YHOO是与其他3个品种走势最不一样的。

最后,以SP500为市场组合,分别计算出3只股票的alpha和beta。


# 计算alpha
> CAPM.alpha(dat_ret[,1:3],dat_ret[,4],Rf=Rf)
                      IBM           GE         YHOO
Alpha: SP500 -0.000752943 0.0003502332 0.0003944279

# 计算beta
> CAPM.beta(dat_ret[,1:3],dat_ret[,4],Rf=Rf)
                  IBM       GE     YHOO
Beta: SP500 0.8218135 1.098877 1.064844

3只股票中,IBM的alpha是最小的,而且是负的,说明IBM落后于市场,买IBM不如直接SP500更好。GE的Beta是最大的,在上升时期beta越大,获得的市场收益也会越大。YHOO从Alpha和Beta上看,虽然与GE接近,但由于标准差,最大回撤等指标过大,会导致波动太大。

综上分析,我们如果配置部分GE和部分YHOO,就可以获得比标普好的收益,但由于GE和YHOO的beta都高于SP500,所以风险也会高于SP500,需要增加新的股票来分散风险,具体的定量分析,将在以后的文章中再进行介绍了。

5. Beta VS Alpha

最后,补充一些Alpha和Beta的说明。Alpha和Beta的认知最早是一个股市起源的概念,是一个关于投资组合的收益率分解的问题

  • Alpha:一般被认为是投资组合的超额收益,也既管理人的能力;
  • Beta:市场风险,最初主要指股票市场的系统性风险

Alpha是平均实际回报和平均预期回报的差额。

  • α>0,表示一基金或股票的价格可能被低估,建议买入。
  • α<0,表示一基金或股票的价格可能被高估,建议卖空。
  • α=0,表示一基金或股票的价格准确反映其内在价值,未被高估也未被低估。

Beta反映了单个证券与整体市场组合的联动性。

  • β>1,攻击性,市场上升时涨幅大。
  • β<1,防御性,市场下跌时跌幅小。
  • β=1,中立性,与市场波动一致。

从资本资产定价模型开始发展到现今,已经有很长的时间了。金融理论在一直发展,继资本资产定价模型之后又一重要的理论突破是套利定价理论,我将在下一篇文章中进行介绍。

本文中,我详细地介绍了资本资产定价模型的金融理论、推到过程、以及R语言实现,用我自己的理解进行阐述。希望能给走在量化道路上的朋友带来入门的指引和帮助,也希望找到像我一样,通过IT转金融的人,让我一起用IT技术+金融的思维在金融市场抢钱吧。

转载请注明出处:
http://blog.fens.me/finance-capm

打赏作者