rm(list=ls())library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(tidyr)
library(tidyquant)## Loading required package: lubridate##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, union## Loading required package: PerformanceAnalytics## Loading required package: xts## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric##
## Attaching package: 'xts'## The following objects are masked from 'package:dplyr':
##
## first, last##
## Attaching package: 'PerformanceAnalytics'## The following object is masked from 'package:graphics':
##
## legend## Loading required package: quantmod## Loading required package: TTR## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo## == Need to Learn tidyquant? ====================================================
## Business Science offers a 1-hour course - Learning Lab #9: Performance Analysis & Portfolio Optimization with tidyquant!
## </> Learn more at: https://university.business-science.io/p/learning-labs-pro </>library(readr)#setwd(“D:/äºæ´²å¤§å¸ä¸èª²è³æ/Portfolio management 2017 Spring”)
dat <- read_csv("capm.csv") %>%
mutate(Date = as.character(Date) %>% as.Date(., "%Y/%m/%d")) %>%
filter(Date>=as.Date("1993-11-01") & Date<= as.Date("1998-11-30")) %>%
rename(rf = `Close-tbill`,
sp500 = `Close-sp500`,
msft = `Close-msft`,
ge = `Close-ge`,
ford = `Close-ford`) %>%
# convert risk-free rate into daily returns
mutate(rf = rf/(100*360)) ##
## -- Column specification --------------------------------------------------------
## cols(
## Date = col_character(),
## `Close-tbill` = col_double(),
## `Close-sp500` = col_double(),
## `Close-msft` = col_double(),
## `Close-ge` = col_double(),
## `Close-ford` = col_double()
## )#
glimpse(dat)## Rows: 1,278
## Columns: 6
## $ Date <date> 1993-11-01, 1993-11-02, 1993-11-03, 1993-11-04, 1993-11-05, 199~
## $ rf <dbl> 8.500000e-05, 8.666667e-05, 8.555556e-05, 8.527778e-05, 8.527778~
## $ sp500 <dbl> 469.10, 468.44, 463.02, 457.49, 459.57, 460.21, 460.33, 463.72, ~
## $ msft <dbl> 2.52, 2.50, 2.45, 2.38, 2.45, 2.45, 2.44, 2.53, 2.55, 2.51, 2.56~
## $ ge <dbl> 6.72, 6.69, 6.64, 6.46, 6.49, 6.46, 6.51, 6.52, 6.47, 6.50, 6.65~
## $ ford <dbl> 8.16, 8.31, 8.24, 7.98, 7.94, 8.02, 7.94, 8.05, 8.07, 8.10, 8.10~tail(dat)## # A tibble: 6 x 6
## Date rf sp500 msft ge ford
## <date> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1998-11-20 0.000121 1164. 28.3 28.6 25.9
## 2 1998-11-23 0.000124 1188. 29.7 29.2 26.1
## 3 1998-11-24 0.000125 1183. 30.3 29.1 25.9
## 4 1998-11-25 0.000124 1187. 31.0 28.8 25.7
## 5 1998-11-27 0.000123 1192. 31.9 28.5 25.8
## 6 1998-11-30 0.000123 1164. 30.4 27.9 25.8#ret4 <- dat %>% select(-rf) %>%
gather(key = stock, value = price, -Date) %>%
group_by(stock) %>%
tq_transmute( mutate_fun = periodReturn,
period = "daily",
type = "arithmetic",
col_rename = "daily.returns") %>%
ungroup() %>%
spread(stock, daily.returns) %>%
bind_cols(., rf = dat$rf) %>%
# subract each returns by risk-free rate
mutate(ford_rf = ford - rf,
ge_rf = ge - rf,
msft_rf = msft - rf,
sp500_rf = sp500 - rf) %>%
# Delete the first row with 0 data
slice(-1) %>%
select(Date, ends_with("_rf"))
#
head(ret4)## # A tibble: 6 x 5
## Date ford_rf ge_rf msft_rf sp500_rf
## <date> <dbl> <dbl> <dbl> <dbl>
## 1 1993-11-02 0.0183 -0.00455 -0.00802 -0.00149
## 2 1993-11-03 -0.00851 -0.00756 -0.0201 -0.0117
## 3 1993-11-04 -0.0316 -0.0272 -0.0287 -0.0120
## 4 1993-11-05 -0.00510 0.00456 0.0293 0.00446
## 5 1993-11-08 0.00999 -0.00471 -0.000085 0.00131
## 6 1993-11-09 -0.0101 0.00765 -0.00417 0.000174ret4.reg <- ret4 %>% lm(formula = cbind(msft_rf, ge_rf, ford_rf) ~ sp500_rf, data = .)
b_hat <- ret4.reg$coefficientsdiagD_hat <- ret4.reg$residuals %>% cov() %>% diag() %>% diag(nrow = length(.))
diagD_hat## [,1] [,2] [,3]
## [1,] 0.0002861862 0.000000e+00 0.0000000000
## [2,] 0.0000000000 9.148714e-05 0.0000000000
## [3,] 0.0000000000 0.000000e+00 0.0002484053sigm2 = var(ret4$sp500_rf)
omega = sigm2*t(b_hat)%*%b_hat + diagD_hat
omega## msft_rf ge_rf ford_rf
## msft_rf 0.0004170745 1.222959e-04 1.046283e-04
## ge_rf 0.0001222959 2.057549e-04 9.775984e-05
## ford_rf 0.0001046283 9.775984e-05 3.320421e-04#————————————————————————— # You can also use OLS formula: beta=inv(X’X)X’Y to get the estimated beta #————————————————————————— n = length(Y) ones = rep(1,1277) X = cbind(ones, Y) X = ret4\(sp500_rf X = as.matrix(X, nco = 1) X = cbind(ones, X) Y = cbind(ret4\)msft_rf, ret4\(ge_rf, ret4\)ford_rf) b_hat.1 = solve(t(X)%%X)%%t(X)%*%Y b_hat.1 b_hat
E_hat = Y - X%%b_hat.1 res_var.1 = diag(t(E_hat)%%E_hat)/(n-2) diagD_hat.1 = diag(res_var.1) diagD_hat.1
omega.1 = var(ret4$sp500_rf)t(b_hat.1)%%b_hat.1 + diagD_hat.1 omega.1 omega