HW6

retdata <- read.csv("m-fac9003.csv")
head(retdata)

stockret <- as.matrix(retdata[, 1:13])
head(stockret,6) # 168 x 13

retdata <- retdata / 100
head(retdata)

mktret <- as.matrix(retdata[, 14])
head(mktret,6) # 168x1

stockret <- as.matrix(retdata[, 1:13])
mktret <- as.matrix(retdata[, 14])
TT <- dim(stockret)[1]

fit <- lm(formula = stockret ~ mktret)
fit

ones <- rep(1, 168)
X <- as.matrix(cbind(ones, mktret))
Y <- stockret
b_hat = solve(t(X)%*%X)%*%t(X) %*% stockret
b_hat 

sigF = as.numeric(var(mktret))
bbeta = as.matrix(fit$coefficients)
bbeta = as.matrix(bbeta[-1,])
bbeta

sigeps = crossprod(fit$residuals)/(TT-2)
sigeps = diag(diag(sigeps))
sigeps

cov_1f = sigF*bbeta%*%t(bbeta)+sigeps
cov_1f

ones = rep(1, TT)
X = as.matrix(cbind(ones, mktret))
Y = stockret
b_hat = solve(t(X)%*%X)%*%t(X) %*% stockret
E_hat = Y - X %*% b_hat
b_hat = as.matrix(b_hat[-1,])
diagD_hat = diag(t(E_hat) %*% E_hat)/(TT-2)
cov_1f.1 = as.numeric(var(mktret))*b_hat%*%t(b_hat) + diag(diagD_hat); 
cov_1f.1