library(tidyr)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(readr)
m_fac9003_3_ <- read_csv("m-fac9003(3).csv")
## Rows: 168 Columns: 14
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (14): AA, AGE, CAT, F, FDX, GM, HPQ, KMB, MEL, NYT, PG, TRB, TXN, SP500
## 
## ℹ 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.
#The monthly series of 3-month Treasury bill rates of the secondary market is used as the proxy for risk-free rate and has been deducted from the monthly returns to obtain the excess returns for stocks and market index.  
#Compute global minimum variance portfolio weights and return by using single index model.
t <- dim(m_fac9003_3_)[1]
market <- m_fac9003_3_[,14]
fac1 <- m_fac9003_3_[,c(-14)]
fac1 <- as.matrix(fac1)
num <- rep(1,t)
xmat <- cbind(num, market)
xmat <- as.matrix(xmat)
betta_hat <- solve(t(xmat)%*%xmat)%*%t(xmat)%*%fac1
E_hat <- fac1 - xmat%*%betta_hat
diagD_hat <- diag(t(E_hat)%*%E_hat)/(t-2)
ret <- apply(fac1, 2, var)
r2 <- 1- diag(t(E_hat)%*%E_hat)/((t-1)*ret)
rsq <- sqrt(diagD_hat)
cov_factor <- as.vector(var(market))*t(betta_hat)%*%betta_hat + diag(diagD_hat)
one.vec <- rep(1,13)
a <- solve(cov_factor)%*%one.vec
b <- t(one.vec)%*%a
gmv <- a / as.numeric(b)
gmv
##            [,1]
## AA   0.03501615
## AGE -0.03373670
## CAT  0.05314427
## F    0.05576422
## FDX  0.06451069
## GM   0.12369835
## HPQ -0.03290554
## KMB  0.28022778
## MEL  0.01901912
## NYT  0.19373120
## PG   0.19019732
## TRB  0.14314200
## TXN -0.09180887
barplot(as.vector(gmv), ylim = c(-0.1, 0.3), names.arg = rownames(gmv), cex.names = 0.5)