library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.0.5
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v ggplot2 3.3.5 v purrr 0.3.4
## v tibble 3.1.5 v dplyr 1.0.7
## v tidyr 1.1.4 v stringr 1.4.0
## v readr 2.1.1 v forcats 0.5.1
## Warning: package 'ggplot2' was built under R version 4.0.5
## Warning: package 'tibble' was built under R version 4.0.5
## Warning: package 'tidyr' was built under R version 4.0.5
## Warning: package 'readr' was built under R version 4.0.5
## Warning: package 'purrr' was built under R version 4.0.5
## Warning: package 'dplyr' was built under R version 4.0.5
## Warning: package 'stringr' was built under R version 4.0.5
## Warning: package 'forcats' was built under R version 4.0.5
## -- 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
##
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.
retdata = FamaFrench_mon_69_98_3stocks_1_
glimpse(retdata)
## 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(retdata)[2]<- 'Mkt_RF'
stock.rets<-retdata %>% select(c(2,6,7,8))/100
glimpse(stock.rets)
## 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(stock.rets)[1]
Mkt.RF<-retdata %>% select(2)/100
fit = lm(formula = cbind(ge,ibm,mobil)~Mkt_RF, data = stock.rets)
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 = 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 0.004070402 0.001823896 0.001825714
## ibm 0.001823896 0.004630742 0.001406268
## mobil 0.001825714 0.001406268 0.004321127
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 0.004070402 0.001823896 0.001825714
## ibm 0.001823896 0.004630742 0.001406268
## mobil 0.001825714 0.001406268 0.004321127
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
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