library(fit.models)

library(robust)

library(readr)

FamaFrench_mon_69_98_3stocks <- read.csv("FamaFrench_mon_69_98_3stocks.csv")

retdata <- read.csv("famaFrench_mon_69_98_3stocks.csv")

attach(retdata)

fit = lm(formula = cbind(ge, ibm, mobil) ~ Mkt.RF + SMB + HML)

coef(fit)

##                       ge        ibm      mobil
## (Intercept)  0.894224681  0.6960039  0.7134499
## Mkt.RF       1.134331448  0.8050676  0.9803403
## SMB         -0.369412445 -0.3099907 -0.3727822
## HML          0.009630701 -0.2981744  0.3726475

covRob(data = fit$residuals, corr = T)

## Call:
## covRob(data = fit$residuals, corr = T)
## 
## Robust Estimate of Correlation: 
##             ge      ibm    mobil
## ge     1.00000  0.04436 -0.24996
## ibm    0.04436  1.00000 -0.07258
## mobil -0.24996 -0.07258  1.00000
## 
## Robust Estimate of Location: 
##       ge      ibm    mobil 
## -0.26939  0.03288 -0.12610

sigF = as.matrix(var(cbind(Mkt.RF, SMB, HML)))

sigF

##           Mkt.RF       SMB       HML
## Mkt.RF 21.150745  4.232563 -5.104485
## SMB     4.232563  8.181134 -1.076021
## HML    -5.104485 -1.076021  7.179737

bbeta = as.matrix(fit$coefficients)

bbeta = t(bbeta[-1, ])

bbeta

##          Mkt.RF        SMB          HML
## ge    1.1343314 -0.3694124  0.009630701
## ibm   0.8050676 -0.3099907 -0.298174424
## mobil 0.9803403 -0.3727822  0.372647515

n = dim(CRSP)

sigeps = (n-1)/(n-4)*var(as.matrix(fit$residuals))

sigeps = diag(diag(as.matrix(sigeps)))

cov_equities = as.vector(bbeta %*% sigF %*% t(bbeta)) + sigeps

cov(cbind(ge, ibm, mobil))

##             ge      ibm    mobil
## ge    40.65659 20.61864 14.01833
## ibm   20.61864 46.21756 11.26794
## mobil 14.01833 11.26794 43.13012