本节书摘来异步社区《量化金融R语言高级教程》一书中的第2章，第2.2节，作者： 【匈牙利】Edina Berlinger（艾迪娜•伯林格） , 等 译者： 高蓉 责编： 胡俊英，更多章节内容可以访问云栖社区“异步社区”公众号查看。
2.2 在R中建模
在接下来的部分中，我们将会学习在R的帮助下如何实现先前讲过的模型。
2.2.1 数据选择
在第4章大数据—高级分析中，我们将会详细讨论获取开源数据及对其进行高效处理的各种知识与方法。在这里，我们仅仅讲述股票价格时间序列和其他相关信息如何获取和如何用于因素模型估计。
我们使用quantmod包来收集数据库。
以下是在R中实现的代码：
library(quantmod)
stocks <- stockSymbols()
然后，我们需要等待几秒钟获取数据，接着可以看到输出：
Fetching AMEX symbols...
Fetching NASDAQ symbols...
Fetching NYSE symbols...
现在，我们得到一个R数据框对象，这个对象涵盖了在各个交易所（如AMEX，NASDAQ，或者NYSE）进行交易的大约6500只股票。为了查看数据集涵盖的变量，我们用str命令：
str(stocks)
'data.frame': 6551 obs. of 8 variables:
$ Symbol : chr "AA-P" "AAMC" "AAU" "ACU" ...
$ Name : chr "Alcoa Inc." "Altisource Asset Management Corp"...
$ LastSale : num 87 1089.9 1.45 16.58 16.26 ...
$ MarketCap: num 0.00 2.44e+09 9.35e+07 5.33e+07 2.51e+07 ...
$ IPOyear : int NA NA NA 1988 NA NA NA NA NA NA ...
$ Sector : chr "Capital Goods" "Finance" "Basic Industries"...
$ Industry : chr "Metal Fabrications" "Real Estate"...
$ Exchange : chr "AMEX" "AMEX" "AMEX" "AMEX" ...
我们可以删掉确实不需要的变量，也可以增加来自不同数据库的市场资本化信息和公司账面价值作为新变量，这些变量我们需要用来估计Fama-French模型：
stocks[1:5, c(1, 3:4, ncol(stocks))]
Symbol LastSale MarketCap BookValuePerShare
1 AA-P 87.30 0 0.03
2 AAMC 985.00 2207480545 -11.41
3 AAU 1.29 83209284 0.68
4 ACU 16.50 53003808 10.95
5 ACY 16.40 25309415 30.13
我们也会需要无风险收益率的时间序列。我们通过计算月度LIBOR市场上的美元利率确定这个序列：
library(Quandl)
Warning message:
package 'Quandl' was built under R version 3.1.0
LIBOR <- Quandl('FED/RILSPDEPM01_N_B',
start_date = '2010-06-01', end_date = '2014-06-01')
Warning message:
In Quandl("FED/RILSPDEPM01_N_B", start_date = "2010-06-01", end_date
= "2014-06-01") : It would appear you aren't using an authentication
token. Please visit http://www.quandl.com/help/r or your usage may be
limited.
因为数据依然分配给LIBOR变量，我们可以忽略这个警告信息。
Quandl包，tseries包和其他收集数据的包会在第4章（大数据—高级分析）中详细讨论。
这个方法也可以用来获取股票价格，并且标普500指数（S&P500指数）可以用来表示市场组合。
我们有一张表，记录了接近5000只股票价格在2010年6月1日到2014年6月1日之间的时间序列。其中第一列和最后几列如下所示：
d <- read.table("data.csv", header = TRUE, sep = ";")
d[1:7, c(1:5, (ncol(d) - 6):ncol(d))]
Date SP500 AAU ACU ACY ZMH ZNH ZOES ZQK ZTS ZX
1 2010.06.01 1070.71 0.96 11.30 20.64 54.17 21.55 NA 4.45 NA NA
2 2010.06.02 1098.38 0.95 11.70 20.85 55.10 21.79 NA 4.65 NA NA
3 2010.06.03 1102.83 0.97 11.86 20.90 55.23 21.63 NA 4.63 NA NA
4 2010.06.04 1064.88 0.93 11.65 18.95 53.18 20.88 NA 4.73 NA NA
5 2010.06.07 1050.47 0.97 11.45 19.03 52.66 20.24 NA 4.18 NA NA
6 2010.06.08 1062.00 0.98 11.35 18.25 52.99 20.96 NA 3.96 NA NA
7 2010.06.09 1055.69 0.98 11.90 18.35 53.22 20.45 NA 4.02 NA NA
如果我们的数据已经储存在硬盘中，那么可以使用read.table函数方便的读取。在第4章（大数据—高级分析）中，我们会讨论如何从互联网收集数据。
现在，我们得到了所需数据：市场组合（标普500）、股票价格、无风险利率（月度LIBOR）。
为了清洗数据库，我们删除了有缺失值、0价格或者负价格的变量。最简单的实现如下：
d <- d[, colSums(is.na(d)) == 0]
d <- d[, c(T, colMins(d[, 2:ncol(d)]) > 0)]
为了使用colMins函数，需要应用matrixStats包。现在，我们可以开始处理数据。
2.2.2 通过主成分分析估计APT
由于识别影响证券收益率的宏观变量很困难（Medvegyev-Szaz，2010，pp.42），我们在实践中使用因子分析相当不易。通常，我们通过主成分分析来寻找驱动收益率变动的潜因子。
在最初下载的6500只股票中，我们可以使用4500只股票数据。其他部分由于缺失值或者0价格的缘故删除了。现在，由于这里用不到日期，还有我们将标普500本身视为一个独立的因子，在主成分分析（principal component analysis，PCA）中不考虑，我们又删除了前两列数据。然后，计算对数收益率：
p <- d[, 3:ncol(d)]
r <- log(p[2:nrow(p), ] / p[1:(nrow(p) - 1), ])
也有其他方法可以计算给定资产的对数收益率，就是PerformanceAnalytics库中的return.calculate(data,method="log")函数。
因为我们的股票数目过于庞大，为了实施PCA，要么需要数据时间长度至少25年，要么需要减少股票数量。使因素模型在几十年间保持稳定，这是不可能的。因此，为了达到说明的目的，我们随机选择了百分之十的股票，并对这个样本计算PCA模型：
r <- r[, runif(nrow(r)) < 0.1]
runif(nrow(r)) < 0.1给出了一个4013维的0—1向量，选出了原表中几乎10%的列（在这个例子中，个数为393）。我们也可以使用以下样本函数，在网站http://stat.ethz.ch/R-manual/ R-devel/library/base/html/sample.html上你可以找到更多细节：
pca <- princomp(r)
结果，我们得到一个princomp类对象。这个对象有8个属性，其中最重要的属性是加载矩阵和sdev属性（包括组成部分的标准差）。第一主成分是数据集方差最大的向量。
我们确认一下主成分的标准差：
plot(pca$sdev)
结果如图2-1所示。
我们可以看出，前5个成分是独立的。因此，我们应该选择5个因子。但是，其他因子的标准差同样显著。所以，不能通过几个因子解释整个市场。
我们可以通过调用factanal函数确认结果。这个函数估计了五因子模型：
factanal(r, 5)
我们可以发现，实施这个计算花费了很多时间。因子分析与PCA相关，但从数学角度看稍复杂一些。结果，我们得到一个factanal类对象。它有很多属性。但是，此时我们仅对以下输出部分有兴趣：
Factor1 Factor2 Factor3 Factor4 Factor5
SS loadings 56.474 23.631 15.440 12.092 6.257
Proportion Var 0.144 0.060 0.039 0.031 0.016
Cumulative Var 0.144 0.204 0.243 0.274 0.290
Test of the hypothesis that 5 factors are sufficient.
The chi square statistic is 91756.72 on 75073 degrees of freedom.The
p-value is 0
结果显示，五因子模型适合数据。但是，可解释的方差仅仅接近30%。这意味着模型应该考虑扩展进其他因子。
2.2.3 Fama-French模型估计
我们有一个包含4015只股票5年价格的数据框架，和包含LIBOR时间序列的LIBOR数据框架。首先，我们需要计算收益率，再与LIBOR利率合并。
第一步，我们删掉数学计算不需要的数据。然后，对保留下来的每一列计算对数收益率：
d2 <- d[, 2:ncol(d)]
d2 <- log(tail(d1, -1)/head(d1, -1))
在计算了对数收益率之后，我们把日期放回到收益率中。然后，在最后一步合并两个数据集：
d <- cbind(d[2:nrow(d), 1], d2)
d <- merge(LIBOR, d, by = 1)
需要提醒读者的是，对数据框架实施merge函数，相当于SQL中的（内）联语句。
结果如下：
print(d[1:5, 1:5])]
Date LIBOR SP500 AAU ACU
2010.06.02 0.4 0.025514387 -0.01047130 0.034786116
2010.06.03 0.4 0.004043236 0.02083409 0.013582552
2010.06.04 0.4 -0.035017487 -0.04211149 -0.017865214
2010.06.07 0.4 -0.013624434 0.04211149 -0.017316450
2010.06.08 0.4 0.010916240 0.01025650 -0.008771986
我们将LIBOR利率转换为日度收益率：
d$LIBOR <- d$LIBOR / 36000
由于LIBOR利率的报价基于货币市场基差——以及（实际天数/360）的天数计算约定——而且时间序列包含使用百分比表示的利率，我们把LIBOR除以36000。现在，我们需要计算Fama-French模型的3个变量。正如在数据选择部分中所讲，我们有股票的数据框：
d[1:5, c(1,(ncol(d) - 3):ncol(d))]
Symbol LastSale MarketCap BookValuePerShare
1 AA-P 87.30 0 0.03
2 AAMC 985.00 2207480545 -11.41
3 AAU 1.29 83209284 0.68
4 ACU 16.50 53003808 10.95
5 ACY 16.40 25309415 30.13
我们得删掉那些没有价格数据的股票：
> stocks = stocks[stocks$Symbol %in% colnames(d),]
我们将市场上限作为一个变量。我们仍需对每只股票计算账面市值比：
stocks$BookToMarketRatio <-
stocks$BookValuePerShare / stocks$LastSale
str(stocks)
'data.frame': 3982 obs. of 5 variables:
$ Symbol : Factor w/ 6551 levels "A","AA","AA-P",..: 14
72...
$ LastSale : num 1.29 16.5 16.4 2.32 4.05 ...
$ MarketCap : num 8.32e+07 5.30e+07 2.53e+07 1.16e+08...
$ BookValuePerShare: num 0.68 10.95 30.13 0.19 0.7 ...
$ BookToMarketRatio: num 0.5271 0.6636 1.8372 0.0819 0.1728 ...
现在，我们需要计算SMB因素和HML因素。为了简化，我们将BIG公司定义为大于平均水平的公司。账面市值比实施同样的原则：
avg_size <- mean(stocks$MarketCap)
BIG <- as.character(stocks$Symbol[stocks$MarketCap > avg_size])
SMALL <- as.character(stocks[stocks$MarketCap < avg_size,1])
这些数组包括了BIG公司和SMALL公司。现在，我们可以定义SMB因素：
d$SMB <- rowMeans(d[,colnames(d) %in% SMALL]) –
rowMeans(d[,colnames(d) %in% BIG])
我们接着定义HML因素：
avg_btm <- mean(stocks$BookToMarketRatio)
HIGH <- as.character(
stocks[stocks$BookToMarketRatio > avg_btm, 1])
LOW <- as.character(
stocks[stocks$BookToMarketRatio < avg_btm, 1])
d$HML <- rowMeans(d[, colnames(d) %in% HIGH]) –
rowMeans(d[, colnames(d) %in% LOW])
第三个因素这样计算：
d$Market <- d$SP500 - d$LIBOR
定义完3个因素，我们在花旗集团（Citi）股票和伊克塞利克斯（EXEL）股票上试一试：
d$C <- d$C - d$LIBOR
model <- glm( formula = "C ～ Market + SMB + HML" , data = d)
GLM（general linear model，一般线性模型）函数作用如下：它将数据和公式作为参数读入。公式是一个形为响应～条件的字符串，其中响应是数据框中的一个变量名，条件指定了模型中的预测子，它包含在数据集中通过操作符“＋”分隔开的变量名。这个函数也可以用于Logistic回归，只是缺省状态设定为线性。
模型输出如下：
Call: glm(formula = "C～Market+SMB+HML", data = d)
Coefficients:
(Intercept) Market SMB HML
0.001476 1.879100 0.401547 -0.263599
Degrees of Freedom: 1001 Total (i.e. Null); 998 Residual
Null Deviance: 5.74
Residual Deviance: 5.364 AIC: -2387
模型概要输出如下：
summary(model)
Call:
glm(formula = "C～Market+SMB+HML", data = d)
Deviance Residuals:
Min 1Q Median 3Q Max
-0.09344 -0.01104 -0.00289 0.00604 2.26882
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.001476 0.002321 0.636 0.525
Market 1.879100 0.231595 8.114 1.43e-15 ***
SMB 0.401547 0.670443 0.599 0.549
HML -0.263599 0.480205 -0.549 0.583
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for gaussian family taken to be 0.005374535)
Null deviance: 5.7397 on 1001 degrees of freedom
Residual deviance: 5.3638 on 998 degrees of freedom
AIC: -2387
Number of Fisher Scoring iterations: 2
结果显示，唯一显著的因素是市场溢价。这表明花旗集团的股票收益率倾向于与整个市场本身共同变动。
使用以下命令可以画出结果：
estimation <- model$coefficients[1]+
model$coefficients[2] * d$Market +
model$coefficients[3]*d$SMB +
model$coefficients[4]*d$HML
plot(estimation, d$C, xlab = "estimated risk-premium",
ylab = "observed riks premium",
main = "Fama-French model for Citigroup")
lines(c(-1, 1), c(-1, 1), col = "red")
图2-2绘出了对花旗集团的Fama-French模型估计的风险溢价。
看图2-2可以发现，收益率中存在一个异常值。我们将这个异常值设为0，看看不考虑它之后的结果：
outlier <- which.max(d$C)
d$C[outlier] <- 0
如果我们运行相同的代码建立模型，并再次计算收益率的估计值和观测值，得到以下结果：
model_new <- glm( formula = "C ～ Market + SMB + HML" , data = d)
summary(model_new)
Call:
glm(formula = "C ～ Market + SMB + HML", data = d)
Deviance Residuals:
Min 1Q Median 3Q Max
-0.091733 -0.007827 -0.000633 0.007972 0.075853
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.0000864 0.0004498 -0.192 0.847703
Market 2.0726607 0.0526659 39.355 < 2e-16 ***
SMB 0.4275055 0.1252917 3.412 0.000671 ***
HML 1.7601956 0.2031631 8.664 < 2e-16 ***
---
Signif. codes: 0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘ ’1
(Dispersion parameter for gaussian family taken to be 0.0001955113)
Null deviance: 0.55073 on 1001 degrees of freedom
Residual deviance: 0.19512 on 998 degrees of freedom
AIC: -5707.4
Number of Fisher Scoring iterations: 2
根据以上结果，所有3个因素均显著。
GLM函数不返回{{R}^{2}}。lm函数对线性回归同样可以得到精确值。我们可以从模型总结中读出r.squared = 0.6446。
结果显示，变量可以解释花旗集团风险溢价中超过64%的变动。我们绘出新结果：
estimation_new <- model_new$coefficients[1]+
model_new$coefficients[2] * d$Market +
model_new$coefficients[3]*d$SMB +
model_new$coefficients[4]*d$HML
dev.new()
plot(estimation_new, d$C, xlab = "estimated risk-premium",ylab =
"observed riks premium",main = "Fama-French model for Citigroup")
lines(c(-1, 1), c(-1, 1), col = "red")
这个例子的输出如图2-3所示。
我们再检验另一只股票EXEL：
d$EXEL <- d$EXEL–d$LIBOR
model2 <- glm( formula = "EXEL～Market+SMB+HML" , data = d)
Call: glm(formula = "EXEL～Market+SMB+HML", data = d)
Coefficients:
(Intercept) Market SMB HML
-0.001048 2.038001 2.807804 -0.354592
Degrees of Freedom: 1001 Total (i.e. Null); 998 Residual
Null Deviance: 1.868
Residual Deviance: 1.364 AIC: -3759
模型小结输出如下：
summary(model2)
Call:
glm(formula = "EXEL～Market+SMB+HML", data = d)
Deviance Residuals:
Min 1Q Median 3Q Max
-0.47367 -0.01480 -0.00088 0.01500 0.25348
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.001773 0.001185 -1.495 0.13515
Market 1.843306 0.138801 13.280 < 2e-16 ***
SMB 2.939550 0.330207 8.902 < 2e-16 ***
HML -1.603046 0.535437 -2.994 0.00282 **
---
Signif. codes: 0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘ ’1
(Dispersion parameter for gaussian family taken to be 0.001357998)
Null deviance: 1.8681 on 1001 degrees of freedom
Residual deviance: 1.3553 on 998 degrees of freedom
AIC: -3765.4
Number of Fisher Scoring iterations: 2
根据以上结果，所有3个因子均显著。
GLM函数并不包括{{R}^{2}}。但lm函数对线性模型也可以得到精确的结果。我们从模型小结中得到r.squared = 0.2723。根据这个结果，我们认为变量可以解释EXEL风险溢价的超过27%的变动。
使用以下命令可以绘出图2-4：
estimation2 <- model2$coefficients[1] +
model2$coefficients[2] * d$Market +
model2$coefficients[3] * d$SMB + model2$coefficients[4] * d$HML
plot(estimation2, d$EXEL, xlab = "estimated risk-premium",
ylab = "observed riks premium",
main = "Fama-French model for EXEL")
lines(c(-1, 1), c(-1, 1), col = "red")
