library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.5 ✓ purrr 0.3.4
## ✓ tibble 3.1.6 ✓ dplyr 1.0.8
## ✓ tidyr 1.2.0 ✓ stringr 1.4.0
## ✓ readr 2.1.2 ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(readr)
FamaFrench_mon_69_98_3stocks_1_ <- read_csv("FamaFrench_mon_69_98_3stocks(1).csv")
## Rows: 360 Columns: 9
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (9): date, Mkt-RF, SMB, HML, RF, ge, ibm, mobil, CRSP
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
glimpse(FamaFrench_mon_69_98_3stocks_1_)
## Rows: 360
## Columns: 9
## $ date <dbl> 196901, 196902, 196903, 196904, 196905, 196906, 196907, 19690…
## $ `Mkt-RF` <dbl> -1.20, -5.82, 2.59, 1.52, 0.02, -7.25, -7.05, 4.65, -2.88, 4.…
## $ SMB <dbl> -0.80, -3.90, -0.28, -0.85, -0.27, -5.31, -3.27, 0.89, 1.20, …
## $ HML <dbl> 1.57, 0.93, -0.45, 0.06, 0.74, -1.15, 1.36, -3.83, -3.24, -3.…
## $ RF <dbl> 0.53, 0.46, 0.46, 0.53, 0.48, 0.51, 0.53, 0.50, 0.62, 0.60, 0…
## $ ge <dbl> -1.1984, -6.0377, 6.6474, 5.9621, -3.5806, -3.8196, -4.3056, …
## $ ibm <dbl> -5.9524, -0.7004, 7.0303, 4.4586, -2.5000, 5.8777, -3.9230, 6…
## $ mobil <dbl> -1.4043, -7.8431, 21.5130, 2.9961, 2.6667, -12.9870, -6.0981,…
## $ CRSP <dbl> -0.6714, -5.3641, 3.0505, 2.0528, 0.5038, -6.7388, -6.5173, 5…
colnames(FamaFrench_mon_69_98_3stocks_1_)[2]<- 'Mkt_RF'
stocks<-FamaFrench_mon_69_98_3stocks_1_ %>% select(c(2,6,7,8))/100
glimpse(stocks)
## Rows: 360
## Columns: 4
## $ Mkt_RF <dbl> -0.0120, -0.0582, 0.0259, 0.0152, 0.0002, -0.0725, -0.0705, 0.0…
## $ ge <dbl> -0.011984, -0.060377, 0.066474, 0.059621, -0.035806, -0.038196,…
## $ ibm <dbl> -0.059524, -0.007004, 0.070303, 0.044586, -0.025000, 0.058777, …
## $ mobil <dbl> -0.014043, -0.078431, 0.215130, 0.029961, 0.026667, -0.129870, …
N <- dim(stocks)[1]
Mkt.RF<-FamaFrench_mon_69_98_3stocks_1_ %>% select(2)/100
ones <- rep(1,N)
X <- as.matrix(cbind(ones, stocks$Mkt_RF))
retdata1 <- as.matrix(FamaFrench_mon_69_98_3stocks_1_[,c(6,7,8)]/100)
b_hat <- solve(t(X)%*%X)%*%t(X)%*%retdata1
E_hat <- retdata1 - X%*%b_hat
b_hat <- as.matrix(b_hat[-1,])
diagD_hat <- diag(t(E_hat)%*%E_hat)/(N-2)
cov_1f.1 <- as.numeric(var(Mkt.RF))*b_hat%*%t(b_hat) + diag(diagD_hat);
cov_1f.1
## ge ibm mobil
## ge 0.004070402 0.001823896 0.001825714
## ibm 0.001823896 0.004630742 0.001406268
## mobil 0.001825714 0.001406268 0.004321127
#compute covariance matrix for three stocks based on factor models
stocks1 <- FamaFrench_mon_69_98_3stocks_1_ %>% select(c(2,3,4,6,7,8))/100
sigF3 <- as.matrix(var(cbind(stocks1$Mkt_RF,
stocks1$SMB,
stocks1$HML)))
X.3 <- cbind(ones, stocks$Mkt_RF, stocks$SMB, stocks$HML)
b_hat.3 <- solve(t(X.3)%*%(X.3))%*%t(X.3)%*%retdata1
E_hat.3 <- retdata1 - X.3%*%b_hat.3
b_hat.3 <- as.matrix(b_hat.3[-1,])
diagD_hat.3 <- diag(t(E_hat.3)%*%E_hat.3)/(N-4)
cov_3f.3 <- as.vector(sigF3)*t(b_hat.3)%*%b_hat.3 + diag(diagD_hat.3)
## Warning in as.vector(sigF3) * t(b_hat.3) %*% b_hat.3: Recycling array of length 1 in vector-array arithmetic is deprecated.
## Use c() or as.vector() instead.
cov_3f.3
## [,1] [,2] [,3]
## [1,] 0.006892504 0.0010366803 -0.0012502397
## [2,] 0.001036680 0.0052477995 -0.0002635496
## [3,] -0.001250240 -0.0002635496 0.0046883562