rm(list = ls()) # Elimina todos los objetos del entorno de trabajo en R
# Leyendo los datos
price <- read.csv(file = "https://raw.githubusercontent.com/Neeraj2308/HM-Model/master/Data.csv")
# Eliminando valores NA
price <- price[complete.cases(price), ]
# Eliminar la columna de fecha
price <- price[, -1]
# Número de valores en la serie de datos
ns <- ncol(price)
# Obtener los rendimientos (se usa el logaritmo de los retornos)
ret <- apply(log(price), 2, diff)
# Mostrar un resumen estadístico de los retornos
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
# Retorno esperado y desviación estándar
er <- apply(ret, 2, mean) # Retorno esperado diario
std <- apply(ret, 2, sd) # Desviación estándar diaria
# Matriz de covarianza
covm <- cov(ret)
# Portafolio de varianza mínima global (GMV)
Am <- rbind(2*covm, rep(1, ns)) # Construcción de la matriz del sistema de ecuaciones
Am <- cbind(Am, c(rep(1, ns), 0)) # Añadir restricción de pesos que sumen 1
b <- c(rep(0, ns), 1) # Vector de restricciones
w.gmv <- solve(Am) %*% b # Resolver el sistema para obtener los pesos óptimos
w.gmv <- w.gmv[-(ns+1), ] # Eliminar el último valor (restricción de lambda)
sum(w.gmv) # Verificar que la suma de los pesos sea 1
## [1] 1
# Retorno esperado y desviación estándar
er <- apply(ret, 2, mean) # Retorno esperado diario
std <- apply(ret, 2, sd) # Desviación estándar diaria
# Matriz de covarianza
covm <- cov(ret)
# Portafolio de varianza mínima global (GMV)
Am <- rbind(2*covm, rep(1, ns)) # Construcción de la matriz del sistema de ecuaciones
Am <- cbind(Am, c(rep(1, ns), 0)) # Agregar restricción de pesos que sumen 1
b <- c(rep(0, ns), 1) # Vector de restricciones
w.gmv <- solve(Am) %*% b # Resolver el sistema para obtener los pesos óptimos
w.gmv <- w.gmv[-(ns+1), ] # Eliminar el último valor (restricción de lambda)
sum(w.gmv) # Verificar que la suma de los pesos sea 1
## [1] 1
# Varianza del portafolio (Varianza mínima)
(var.gmv <- t(w.gmv) %*% covm %*% w.gmv)
## [,1]
## [1,] 0.000128295
# Retorno esperado del portafolio de varianza mínima global
(ret.gmv <- t(w.gmv) %*% er)
## [,1]
## [1,] 0.0005783192
## Solución alternativa
uv <- rep(1, ns) # Vector unitario
w.gmv2 <- (solve(covm) %*% uv) * c(solve( (t(uv) %*% solve(covm) %*% uv)))
# Derivamos la frontera eficiente para un retorno dado (es decir, el menor riesgo para un retorno dado)
u0 <- er[3] # Retorno que queremos lograr en nuestro portafolio
if(u0 < ret.gmv) {
message("# u0 debe ser mayor que el retorno del portafolio de varianza mínima")
}
M <- cbind(er, uv) # Matriz con los retornos esperados y un vector de unos
B <- t(M) %*% solve(covm) %*% M # Matriz de coeficientes para la optimización
mu.tilde <- c(u0, 1) # Vector de restricciones
# Pesos del portafolio eficiente
(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
# Varianza del portafolio
(var.ep <- t(w.ep) %*% covm %*% w.ep) # Menor que la varianza del portafolio de mercado, pero mayor que la del portafolio de varianza mínima
## [,1]
## [1,] 0.0001923228
w1 <- .4
w <- seq(from = -.5, to = 1.5, by = .01) # Varias combinaciones de pesos
# Primer método
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 = "Riesgo del Portafolio", ylab = "Retorno del Portafolio", pch = 16,
main = "Frontera Eficiente")
# Segundo método (opcional)
cov.gmv.ep <- t(w.gmv) %*% covm %*% w.ep # Covarianza entre el portafolio de mínima varianza y el portafolio eficiente
eff2 <- function(w1) {
ret.ef <- w1 * ret.gmv + (1 - w1) * ret.ep # Retorno esperado del portafolio eficiente
var.ef <- w1^2 * var.gmv + (1 - w1)^2 * var.ep + 2 * w1 * (1-w1) * cov.gmv.ep # Varianza del portafolio eficiente
c(w = w1, ret = ret.ef * 248 , sd = sqrt(var.ef * 248)) # Resultados anualizados
}
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 ) # Retorno del portafolio de máxima razón de Sharpe
## [,1]
## [1,] 0.0006424395
(var.sr <- t(w.sr) %*% covm %*% w.sr) # Varianza del portafolio de máxima razón de Sharpe
## [,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)
