实证资产定价中的市场Beta估计
薛英杰 / 2022-11-22
在金融经济学中,单只股票对市场组合的暴露反映了该只股票承担的市场风险高低。根据资本资产定价模型(Sharpe 1964; Lintner 1965; Mossin 1966),资产期望收益的横截面变化应该是资产超额收益和市场组合超额收益协方差的函数,而市场超额收益对股票超额收益的回归系数被称为市场\(\beta\)。这儿我们不去讨论资本资产定价模型(CAPM)的理论基础,仅仅说明市场\(\beta\)在实证上如何估计。
固定\(\beta\)估计
为估计资本资产定价模型系数,我们通过以下回归模型来估计:
\[ r_{i,t}-r_{f,t}=\alpha_i+\beta_i(r_{m,t}-r_{f,t})+\epsilon_{i,t} \]
其中,\(r_{i,t}-r_{f,t}\)是股票超额收益,\(r_{m,t}-r_{f,t}\)是市场组合超额收益。
R 提供了简单的线性回归模型函数
lm()
来估计回归模型,具体代码如下:##加载程序包 pacman::p_load(data.table,stringr,dplyr,Tushare,reticulate,lubridate,mongolite,future.apply,zoo,tidyr, haven,readxl,foreign,scales,slider,furrr,tibble,ggplot2) ##开启Mongo数据库 system('net start MongoDB',intern = F,wait = T)
## [1] 2
##读入无风险收益数据 mogFret<-mongo(collection = "DRFRet",db="IntrestRate") Retf<-mogFret$find()%>% mutate(Tradedate=str_sub(TradingDate,1,7), MRet=(1+DRFRet)^(365/12)-1)%>% group_by(Tradedate)%>% summarise(rf=mean(MRet)) ##读入市场组合数据 mogfactor<-mongo(collection = "monthlyfactor",db="investfactor") Marketportfio<-mogfactor$find(query = '{"market":"全市场"}',fields ='{"_id":0,"Tradedate":1,"MKT":1}' ) ##读入股票数据并拼接所有数据 mogmogth<-mongo(collection = "qfqtradedatamonthly",db="stockdata") mothdata<-mogmogth$find(query='{}', fields='{"_id":0,"ts_code":1,"trade_date":1,"close":1,"pre_close":1}')%>% mutate(MRet=close/pre_close-1, trade_date=as.Date(trade_date,"%Y%m%d"), Tradedate=str_sub(trade_date,1,7))%>% merge(Marketportfio,by="Tradedate",all.x = T)%>% merge(Retf,by="Tradedate",all.x = T)%>% filter(Tradedate>="1995-01")%>% arrange(ts_code,Tradedate)%>% mutate(ret_excess=MRet-rf, mkt_excess=MKT-rf)%>% select(ts_code,Tradedate,trade_date,ret_excess,mkt_excess)%>%tibble::as_tibble() ##计算平安银行的股票Beta fit <- lm(ret_excess ~ mkt_excess, data = mothdata |> filter(ts_code == "000001.SZ") ) summary(fit)
## ## Call: ## lm(formula = ret_excess ~ mkt_excess, data = filter(mothdata, ## ts_code == "000001.SZ")) ## ## Residuals: ## Min 1Q Median 3Q Max ## -0.38942 -0.06252 -0.00644 0.04893 0.68285 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 0.002797 0.005839 0.479 0.632 ## mkt_excess 0.821452 0.062841 13.072 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 0.1069 on 334 degrees of freedom ## Multiple R-squared: 0.3384, Adjusted R-squared: 0.3365 ## F-statistic: 170.9 on 1 and 334 DF, p-value: < 2.2e-16
lm()
返回的对象包含所有我们平时关心的线性模型信息,summary()
返回估计系数的描述,以上结果表明平安银行的股票市场\(\hat\beta\)的估计值为1.15。资本资产定价模型表明风险资产的超额收益与投资该资产承担的风险成正比,也就是说,投资者购买平安银行获得的超额收益可以表示为单位风险的价格与风险大小的乘积,而市场超额收益就是单位风险的价格,市场\(\beta\)则反映了投资者承担的风险大小。
滚动窗口估计
投资者在购买风险后,其承担的风险并不是一成不变的,会大小随市场环境的变化而变化,而采用全样本估计的市场\(\beta\)只能反映投资者在样本期内承担风险的平均大小,并不能描绘其投资该股票承担风险的动态变化。为了更加精确地估计投资者承担市场风险的动态变化,我们采用滚动窗口估计法,在每一个小的子样本内来估计该时段的市场\(\beta\),具体如下:
##定义市场Beta的估计函数 estimate_capm <- function(data, min_obs = 1) { if (nrow(data) < min_obs) { beta <- as.numeric(NA) } else { fit <- lm(ret_excess ~ mkt_excess, data = data) beta <- as.numeric(coefficients(fit)[2]) } return(beta) } ## 定义滚动估计函数 roll_capm_estimation <- function(data, months, min_obs) { data <- data |> arrange(trade_date) betas <- slide_period_vec( .x = data, .i = data$trade_date, .period = "month", .f = ~ estimate_capm(., min_obs), .before = months - 1, .complete = FALSE ) return(tibble( month = unique(data$trade_date), beta = betas )) } beta_estimate = mothdata %>% filter(ts_code == "000001.SZ")%>% mutate(roll_capm_estimation(cur_data(),months = 60,min_obs = 48))%>%na.omit()
## Warning: There was 1 warning in `mutate()`. ## ℹ In argument: `roll_capm_estimation(cur_data(), months = 60, min_obs = 48)`. ## Caused by warning: ## ! `cur_data()` was deprecated in dplyr 1.1.0. ## ℹ Please use `pick()` instead.
tail(beta_estimate,60)
## # A tibble: 60 × 7 ## ts_code Tradedate trade_date ret_excess mkt_excess month beta ## <chr> <chr> <date> <dbl> <dbl> <date> <dbl> ## 1 000001.SZ 2018-03 2018-03-30 -0.0993 0.0374 2018-03-30 0.505 ## 2 000001.SZ 2018-04 2018-04-27 -0.00797 -0.0494 2018-04-27 0.499 ## 3 000001.SZ 2018-05 2018-05-31 -0.0651 -0.0358 2018-05-31 0.487 ## 4 000001.SZ 2018-06 2018-06-29 -0.111 -0.130 2018-06-29 0.340 ## 5 000001.SZ 2018-07 2018-07-31 0.0334 -0.00742 2018-07-31 0.345 ## 6 000001.SZ 2018-08 2018-08-31 0.0730 -0.0832 2018-08-31 0.316 ## 7 000001.SZ 2018-09 2018-09-28 0.0885 -0.00851 2018-09-28 0.305 ## 8 000001.SZ 2018-10 2018-10-31 -0.0150 -0.0983 2018-10-31 0.313 ## 9 000001.SZ 2018-11 2018-11-30 -0.0529 0.0519 2018-11-30 0.313 ## 10 000001.SZ 2018-12 2018-12-28 -0.0972 -0.0476 2018-12-28 0.314 ## # … with 50 more rows
需要强调的是,在上述代码中,滚动窗口估计使用了slider包, 其在各种数据上运行非常稳定,而且速度比较快,学习教程可以参考这儿。
我们可以多估计几只股票的市场\(\beta\),了解不同股票在市场风险上的暴露差异,下面选择平安银行(000001.SZ)、万科(000002.SZ)、中炬高新(600872.SH)、三一重工(600031.SH),分析市场\(\beta\)的动态变化。
examples <- tribble( ~ts_code, ~company, "000001.SZ", "平安银行", "000002.SZ", "万科", "600872.SH", "中炬高新", "600031.SH", "三一重工" ) beta_examples = mothdata |> inner_join(examples, by = "ts_code")|> group_by(ts_code) |> mutate(roll_capm_estimation(cur_data(), months = 60, min_obs = 48)) |> ungroup() |> select(ts_code, company, trade_date, beta) |> na.omit() beta_examples |> ggplot(aes( x = trade_date, y = beta, color = company, linetype = company)) + geom_line() + labs( x = NULL, y = NULL, color = NULL, linetype = NULL, title = "Monthly beta estimates for example stocks using 20 years of data" )
并行滚动估计
使用月度数据,计算量并不大,很快就能计算出来,但是当数据是日度频率时,分组计算会特别慢,因此,可以使用并行计算来做滚动估计,具体如下:
plan(multisession, workers = availableCores()) pdat<-mothdata|> select(ts_code,trade_date,ret_excess,mkt_excess)|> nest(data=c(trade_date,ret_excess,mkt_excess)) beta_daily <- pdat|> mutate(beta_daily = future_map( data, ~ roll_capm_estimation(., months = 60, min_obs = 50) )) |> unnest(c(beta_daily)) |> select(ts_code,month,beta_daily = beta) |> na.omit()
在不并行的代码中,我们使用group_by()
分组后进行滚动估计,而在并行计算中,我们需要使用nest()
函数按照股票代码ts_code
对数据进行嵌套,然后对嵌套数据在不同进程上进行运行来获得结果,嵌套前后的数据如下:
pdat<-mothdata|>
select(ts_code,trade_date,ret_excess,mkt_excess)
pdat|>head()
## # A tibble: 6 × 4
## ts_code trade_date ret_excess mkt_excess
## <chr> <date> <dbl> <dbl>
## 1 000001.SZ 1995-01-27 0.000175 -0.131
## 2 000001.SZ 1995-02-28 0.00299 -0.0208
## 3 000001.SZ 1995-03-31 0.00668 0.136
## 4 000001.SZ 1995-04-28 -0.105 -0.115
## 5 000001.SZ 1995-05-31 0.00580 0.150
## 6 000001.SZ 1995-06-30 -0.0698 -0.101
pdat|>
nest(data=c(trade_date,ret_excess,mkt_excess))%>%head(10)
## # A tibble: 10 × 2
## ts_code data
## <chr> <list>
## 1 000001.SZ <tibble [336 × 3]>
## 2 000002.SZ <tibble [332 × 3]>
## 3 000003.SZ <tibble [84 × 3]>
## 4 000004.SZ <tibble [329 × 3]>
## 5 000005.SZ <tibble [322 × 3]>
## 6 000006.SZ <tibble [333 × 3]>
## 7 000007.SZ <tibble [314 × 3]>
## 8 000008.SZ <tibble [332 × 3]>
## 9 000009.SZ <tibble [335 × 3]>
## 10 000010.SZ <tibble [291 × 3]>