rm(list = ls())
#reading Data
price <- read.csv(file = "https://raw.githubusercontent.com/Neeraj2308/HM-Model/master/Data.csv")
#removing NA
price <- price[complete.cases(price), ]
#remove date
price <- price[, -1]
#number of securities
ns <- ncol(price)
#Getting return (we have taken log return)
ret <- apply(log(price), 2, diff)
sapply(as.data.frame(ret), summary)
## Infra Bank Auto Metal
## Min. -0.0966787371 -0.1348483039 -0.1101259230 -0.1257999876
## 1st Qu. -0.0089354252 -0.0090002003 -0.0072087712 -0.0081572173
## Median 0.0008721455 0.0007799968 0.0011342333 0.0010685098
## Mean 0.0005640156 0.0006519906 0.0007143575 0.0004013688
## 3rd Qu. 0.0105660259 0.0108416426 0.0092496360 0.0101112911
## Max. 0.1980336255 0.1754832373 0.1062658302 0.1261595169
#expected return and sd
er <- apply(ret, 2, mean) #daily expected return
std <- apply(ret, 2, sd) #daily sd
#covariance matrix
covm <- cov(ret)
#Global minimum variance (gmv) portfolio
Am <- rbind(2*covm, rep(1, ns))
Am <- cbind(Am, c(rep(1, ns), 0))
b <- c( rep(0, ns), 1)
w.gmv <- solve(Am) %*% b
w.gmv <- w.gmv[-(ns+1), ] #last value is lambda constraint. Hence not relevant
sum(w.gmv) #sum of weights are 1
## [1] 1
#variance of portfolio (Minimum variance)
(var.gmv <- t(w.gmv) %*% covm %*% w.gmv)
## [,1]
## [1,] 0.000128295
#expected return on minimum global variance portfolio
(ret.gmv <- t(w.gmv) %*% er)
## [,1]
## [1,] 0.0005783192
##alternate solution
uv <- rep(1, ns) #unit vector
w.gmv2 <- (solve(covm) %*% uv) * c(solve( (t(uv) %*% solve(covm) %*% uv)))
#we derived efficient for given return (ie. minimum risk for given return)
u0 <- er[3] #return that we want to achieved on our portfolio
if(u0 < ret.gmv) {
message("#u0 should be greater than return on minimum variance portfolio")
}
M <- cbind(er, uv)
B <- t(M) %*% solve(covm) %*% M
mu.tilde <- c(u0, 1)
#weight of efficient portfolio
(w.ep <- solve(covm) %*% M %*% solve(B) %*% mu.tilde)
## [,1]
## Infra -0.3860598
## Bank 0.0960876
## Auto 1.1236778
## Metal 0.1662944
#portfolio expected return
(ret.ep <- t(er) %*% w.ep) #equals to Auto expected return
## [,1]
## [1,] 0.0007143575
#portofolio variance
(var.ep <- t(w.ep) %*% covm %*% w.ep) # less than auto var but more than minimum variance portfolio
## [,1]
## [1,] 0.0001923228
w1 <- .4
w <- seq(from = -.5, to = 1.5, by = .01) #various combination of weights
#first method
eff <- function(w1) {
z <- w1 * w.gmv + (1 - w1) * w.ep
c(ret = t(z) %*% er * 248 , sd = sqrt(t(z) %*% covm %*% z * 248))
}
comb <- cbind(w, t(sapply(w, FUN = eff)))
efficient <- comb[, "ret"] > c(ret.gmv * 248)
xlim <- c(min(comb[, "sd"]), max(comb[, "sd"]))
ylim <- c(min(comb[, "ret"]), max(comb[, "ret"]))
col <- ifelse(efficient, "blue", "red")
plot(comb[ , "ret"] ~ comb[, "sd"], col = col, xlim = xlim, ylim = ylim,
xlab = "Portfolio RisK", ylab = "Portfolio Return", pch = 16, main = "Efficient Frontier",
cex = .7)
#second method (optional)
cov.gmv.ep <- t(w.gmv) %*% covm %*% w.ep
eff2 <- function(w1) {
ret.ef <- w1 * ret.gmv + (1 - w1) * ret.ep
var.ef <- w1^2 * var.gmv + (1 - w1)^2 * var.ep + 2 * w1 * (1-w1) * cov.gmv.ep
c(w = w1, ret = ret.ef * 248 , sd = sqrt(var.ef * 248))
}
comb2 <- t(sapply(w, eff2))
#SR = (mu - rf) / sigma
rfree <- .00 / 248 #daily risk free return
#assuming zero risk free rate of interest
(w.sr <- solve(covm) %*% (er - rfree) / c(t(uv) %*% solve(covm) %*% (er - rfree))) #weights of optimal portfolio
## [,1]
## Infra -0.16441270
## Bank 0.05490457
## Auto 0.81169580
## Metal 0.29781232
(ret.sr <- t(w.sr) %*% er ) #return
## [,1]
## [1,] 0.0006424395
(var.sr <- t(w.sr)%*% covm %*% w.sr) #portfolio variance
## [,1]
## [1,] 0.0001425195
#plotting results
comb <- rbind(comb, c(NA, ret.sr * 248, sqrt(var.sr * 248)))
col <- c(ifelse(efficient, "blue", "red"), "green")
xlim <- c(0.15, max(comb[, "sd"]))
ylim <- c(0.10, max(comb[, "ret"]))
cex <- c(ifelse(efficient, .6, .6), 1.2)
pch <- c(rep(1, nrow(comb) - 1), 16)
plot(comb[ , "ret"] ~ comb[, "sd"], col = col, xlim = xlim, ylim = ylim,
xlab = "Portfolio RisK", ylab = "Portfolio Return", pch = pch, main = "Efficient Frontier",
cex = cex)
abline(a = rfree, b = ret.sr * 248 / sqrt(var.sr * 248), lty = 2)
