#Read data for 10 stocks:
a1 <- read.table("AMZN.csv", sep=",", header=TRUE)
a2 <- read.table("COST.csv", sep=",", header=TRUE)
a3 <- read.table("CVX.csv", sep=",", header=TRUE)
a4 <- read.table("DIS.csv", sep=",", header=TRUE)
a5 <- read.table("GILD.csv", sep=",", header=TRUE)
a6 <- read.table("GS.csv", sep=",", header=TRUE)
a7 <- read.table("LMT.csv", sep=",", header=TRUE)
a8 <- read.table("MMM.csv", sep=",", header=TRUE)
a9 <- read.table("PG.csv", sep=",", header=TRUE)
a10 <- read.table("TSLA.csv", sep=",", header=TRUE)

a11 <- read.table("^GSPC.csv", sep=",", header=TRUE)

cov_matrix <- cov(a)

#Convert adjusted close proces into returns:
r1 <- (a1$Adj.Close[-1]-a1$Adj.Close[-nrow(a1)])/a1$Adj.Close[-nrow(a1)]

r2 <- (a2$Adj.Close[-1]-a2$Adj.Close[-nrow(a2)])/a2$Adj.Close[-nrow(a2)]

r3 <- (a3$Adj.Close[-1]-a3$Adj.Close[-nrow(a3)])/a3$Adj.Close[-nrow(a3)]

r4 <- (a4$Adj.Close[-1]-a4$Adj.Close[-nrow(a4)])/a4$Adj.Close[-nrow(a4)]

r5 <- (a5$Adj.Close[-1]-a5$Adj.Close[-nrow(a5)])/a5$Adj.Close[-nrow(a5)]

r6 <- (a6$Adj.Close[-1]-a6$Adj.Close[-nrow(a6)])/a6$Adj.Close[-nrow(a6)]

r7 <- (a7$Adj.Close[-1]-a7$Adj.Close[-nrow(a7)])/a7$Adj.Close[-nrow(a7)]

r8 <- (a8$Adj.Close[-1]-a8$Adj.Close[-nrow(a8)])/a8$Adj.Close[-nrow(a8)]

r9 <- (a9$Adj.Close[-1]-a9$Adj.Close[-nrow(a9)])/a9$Adj.Close[-nrow(a9)]

r10 <- (a10$Adj.Close[-1]-a10$Adj.Close[-nrow(a10)])/a10$Adj.Close[-nrow(a10)]


r <- (a11$Adj.Close[-1]-a11$Adj.Close[-nrow(a11)])/a11$Adj.Close[-nrow(a11)]

rr <- as.data.frame(cbind(r1,r2,r3,r4,r5,r6,r7,r8,r9,r10))

#Compute the variance covariance matrix:
cov_matrix <- cov(rr)

#Risk free asset:
rf <- 0.001

#Compute the mean returns for each stock:
means <- colMeans(rr)

#Compute R vector:
R <- means - rf

#Compute z vector:
z <- solve(cov_matrix) %*% R

#Finally compute the portfolio weights (point of tangency):
x <- z/sum(z)

#Compute expected return and standard deviation of the point of tangency:
rrr <- t(x) %*% means

sigma_port <- (t(x) %*% cov_matrix %*% x)^.5

#Place all of them on the graph:
sigmas <- diag(cov_matrix)^.5

plot(sigmas, means, cex=3, xlim=c(0,.18), pch=19)

points(sigma_port, rrr, cex=3, pch=19, col="blue")

#20% Rf + 80% tangency point:
rfg <- 0.2*rf + 0.8* rrr
sdfg <- 0.8*sigma_port

points(sdfg, rfg, cex=3, pch=19, col="green")


#-20% Rf + 120% tangency point:
rfg1 <- -0.2*rf + 1.2* rrr
sdfg1 <- 1.2*sigma_port

points(sdfg1, rfg1, cex=3, pch=19, col="purple")



#An example using the excess return to beta ratio using 5 stocks.
#Run the regressions:
q1 <- lm(r1 ~ r)
q2 <- lm(r2 ~ r)
q3 <- lm(r3 ~ r)
q4 <- lm(r4 ~ r)
q5 <- lm(r5 ~ r)

#Get alphas:
alpha <- c(summary(q1)$coef[1],summary(q2)$coef[1],summary(q3)$coef[1],summary(q4)$coef[1],summary(q5)$coef[1])

#Get betas:
beta <- c(summary(q1)$coef[2],summary(q2)$coef[2],summary(q3)$coef[2],summary(q4)$coef[2],summary(q5)$coef[2])

#Compute mean return:
rbar <- alpha + beta*mean(r)

#Get sigma^2:
sigma2 <- c(summary(q1)$sigma^2,summary(q2)$sigma^2,summary(q3)$sigma^2,summary(q4)$sigma^2,summary(q5)$sigma^2)


#Compute ratio:
ratio <- (rbar-.001)/beta

#Put them together (see handout table):
xx <- cbind(1:5,alpha,beta,rbar, sigma2,ratio)

#Order the table based on the excess return to beta ratio:
aaa <- xx[order(-ratio),]


#Create the extra columns (the following is just an example):
q <- seq(1:5)

qq <- c(sum(q[1]),sum(q[1:2]),sum(q[1:3]),sum(q[1:4]),sum(q[1:5]))

q1 <- rep(0,5)
for(i in (1:5)){q1[i] <- sum(q[1:i])}