#========================================
# FamaFrench_mon_69_98_3stocks
# one factor model
# ff three factor model
#========================================
library(tidyverse)
## Registered S3 methods overwritten by 'ggplot2':
## method from
## [.quosures rlang
## c.quosures rlang
## print.quosures rlang
## Registered S3 method overwritten by 'rvest':
## method from
## read_xml.response xml2
## ── Attaching packages ────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──
## ✔ ggplot2 3.1.1 ✔ purrr 0.3.2
## ✔ tibble 2.1.1 ✔ dplyr 0.8.5
## ✔ tidyr 1.0.2 ✔ stringr 1.4.0
## ✔ readr 1.3.1 ✔ forcats 0.4.0
## ── Conflicts ───────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
retdata = read_csv('FamaFrench_mon_69_98_3stocks.csv')
## Parsed with column specification:
## cols(
## date = col_double(),
## `Mkt-RF` = col_double(),
## SMB = col_double(),
## HML = col_double(),
## RF = col_double(),
## ge = col_double(),
## ibm = col_double(),
## mobil = col_double(),
## CRSP = col_double()
## )
colnames(retdata)[2]<- 'Mkt_RF'# Replace 'Mkt-RF' with 'Mkt_RF';
# attach(retdata)
#Below we use two different approaches to estimate covariance matrix
#===========================================================
# Single index model to compute covariance matrix
#===========================================================
# Method 1: by "lm" function
#===========================
stock.rets<-retdata %>% select(c(2,6,7,8))/100
glimpse(stock.rets)
## Observations: 360
## Variables: 4
## $ Mkt_RF <dbl> -0.0120, -0.0582, 0.0259, 0.0152, 0.0002, -0.0725, -0.070…
## $ ge <dbl> -0.011984, -0.060377, 0.066474, 0.059621, -0.035806, -0.0…
## $ ibm <dbl> -0.059524, -0.007004, 0.070303, 0.044586, -0.025000, 0.05…
## $ mobil <dbl> -0.014043, -0.078431, 0.215130, 0.029961, 0.026667, -0.12…
N <- dim(stock.rets)[1]
#Mkt.RF<-retdata %>% select(2)/100
fit = lm(formula = cbind(ge,ibm,mobil)~Mkt_RF, data = stock.rets)
Mkt_RF <- retdata[, 2]
sigF = as.numeric(var(Mkt_RF))
bbeta = as.matrix(fit$coefficients)
bbeta = as.matrix(bbeta[-1,])
bbeta
## [,1]
## ge 1.0580825
## ibm 0.8149949
## mobil 0.8158072
sigeps = crossprod(fit$residuals)/(N-2)
# sigeps = as.matrix(var(fit$residuals)) # you can use this way too
sigeps = diag(diag(sigeps))
sigeps
## [,1] [,2] [,3]
## [1,] 0.001702494 0.000000000 0.000000000
## [2,] 0.000000000 0.003225874 0.000000000
## [3,] 0.000000000 0.000000000 0.002913458
cov_1f = sigF*bbeta%*%t(bbeta)+sigeps
cov_1f
## ge ibm mobil
## ge 23.68078 18.23896 18.25714
## ibm 18.23896 14.05191 14.06268
## mobil 18.25714 14.06268 14.07961
#===================================
#Method 2: by formula "inv(X'X)*X'Y"
#===================================
ones = rep(1,N)
X = as.matrix(cbind(ones, stock.rets$Mkt_RF))
retdata1 = as.matrix(retdata[,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 23.68078 18.23896 18.25714
## ibm 18.23896 14.05191 14.06268
## mobil 18.25714 14.06268 14.07961
#===================================================================
# Using FF 3 factor model to compute covariance matrix
#===================================================================
# Method 1: by "lm" function
#============================
stock.rets<-retdata %>% select(c(2,3,4,6,7,8))/100
fit3 = lm(formula = cbind(ge, ibm, mobil)~Mkt_RF + SMB + HML, data=stock.rets)
sigF3 = as.matrix(var(cbind(stock.rets$Mkt_RF,
stock.rets$SMB,
stock.rets$HML)))
bbeta3 = as.matrix(fit3$coefficients)
bbeta3 = bbeta3[-1,]
bbeta3
## ge ibm mobil
## Mkt_RF 1.134331448 0.8050676 0.9803403
## SMB -0.369412445 -0.3099907 -0.3727822
## HML 0.009630701 -0.2981744 0.3726475
sigeps3 = crossprod(fit3$residuals)/(N-4)
sigeps3 = diag(diag(sigeps3))
cov_3f = t(bbeta3) * sigF3 * (bbeta3) + sigeps3
cov_3f
## Mkt_RF SMB HML
## ge 4.332517e-03 -1.258773e-04 -4.819330e-06
## ibm -1.258773e-04 3.199233e-03 -1.196042e-05
## mobil -4.819330e-06 -1.196042e-05 2.841930e-03
#===================================
#Method 2: by formula "inv(X'X)*X'Y"
#===================================
X.3 = cbind(ones, stock.rets$Mkt_RF, stock.rets$SMB, stock.rets$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 = t(b_hat.3)*sigF3*b_hat.3 + diag(diagD_hat.3)
cov_3f.3
##
## ge 4.332517e-03 -1.258773e-04 -4.819330e-06
## ibm -1.258773e-04 3.199233e-03 -1.196042e-05
## mobil -4.819330e-06 -1.196042e-05 2.841930e-03
cov_3f
## Mkt_RF SMB HML
## ge 4.332517e-03 -1.258773e-04 -4.819330e-06
## ibm -1.258773e-04 3.199233e-03 -1.196042e-05
## mobil -4.819330e-06 -1.196042e-05 2.841930e-03