Optimización de portafolios en R

A continuación, entregamos un análisis de optimización de portafolio con programación R. Aquí se busca, en términos simples, minimizar la volatilidad de un portafolio dado un nivel de rentabilidad o en su defecto, dado un nivel de rentabilidad maximizamos el ratio de sharpe que corresponde a la rentabilidad por unidad de riesgo.

Mediante la simulación de portafolios y la obtención de la frontera eficiente lograrás obtener los gráficos que adjunto a continuación, y a su vez, puntos de interés como es el portafolio con máximo ratio de sharpe así como el portafolio de mínima varianza. Adicionalmente, podrás graficar los portafolios en función del ratio de sharpe y podrás visualizar el cómo afecta la incorporación y combinación de acciones a la eficiencia de tu portafolio.

Este modelo no solamente da luces en el ámbito de las inversiones si no que también de la gestión, a modo de ejemplo, si una empresa desea llevar a cabo proyectos que son independientes pero con baja probabilidad de éxito, le conviene llevar a cabo varios proyectos puesto que reduce la probabilidad de pérdida (diversificación), para ver más detalle te adjuntamos una nota al respecto en el siguiente link.

Simulación portafolios para obtener retornos y volatilidad.

Frontera eficiente mediante optimización.

Ponderación de acciones en función del ratio de sharpe

Gráfico comparativo normalizado.

De ante mano pedimos las disculpas por la no utilización de acentos en el código, lo cual puede llevar a confusión, en algunos casos, en su explicación.

Adjuntamos algunas consideraciones:

Obtenemos precios diarios desde Yahoo, en particular los precios ajustados, por lo que debes estar conectado a internet. El análisis es diario y se anualiza por 250 días que corresponderían a los días que se transa en un año. No se considera tasa de libre riesgo o en su defecto es 0%. Debes instalar las librerías que se encuentran descritas en el código mediante la función install.packages() y activarlas con library(). Finalmente, utilizamos la versión de R 4.3.1.

Adjuntamos el código con una breve explicación de la programación:

library(rvest)

library(dplyr)

library(quantmod)

library(purrr)

library(tidyverse)

library(timetk)

library(lubridate)

library(PerformanceAnalytics)

library(highcharter)

library(nloptr)

library(tidyr)

#Extraccion de datos desde Yahoo

read_html_tables <- function(url){

  webpage <- read_html(url)

  tables <- html_table(webpage, fill=T)

  return(tables)

}

#se define la url en donde se encuentra la data

url <-"https://es.wikipedia.org/wiki/%C3%8Dndice_de_Precio_Selectivo_de_Acciones"

#Se asignan tablas a la variable "tables"

tables <- read_html_tables(url)

#######Importate#####

#seleccinamos las acciones que deseamos analizar

#Se pueden añadir más acciones

#Ojo que las acciones deben estar en la url definidas puesto que de aquí se obtiene la data

# A modo de ejemplo, queremos analizar más acciones, implica que debemos cambiar esto

# c("COPEC","CHILE","CENCOSUD","FALABELLA", "ENTEL", "IAM", "PARAUCO") etc

# si cambias el tamaño de este vector, el programa lo detecta y no debes añadir nada mas al codigo

ipsa <- c("COPEC","CHILE","CENCOSUD", "IAM","PARAUCO")

# Fecha de inicio de extraccion de data

start_date <- "2020-01-01"

# Fecha de terrmino de extraccion de data

end_date <- "2024-08-31"

#Bajando data de Yahoo a traves de una funcion

download_data <- function(symbol){

  stock_data <- getSymbols(Symbols = paste0(symbol, ".SN"),

                           src = "yahoo",

                           from =start_date,

                           to = end_date,

                           auto.assign= FALSE) %>%

    data.frame()%>%

    rownames_to_column(var="Date")%>%

    mutate(Symbol = symbol)

  return(stock_data)

}

# la funcion map recoge un simbolo en el IPSA y lo aplica a la funcion download_data()

# actual como una "for" es decir, para cada elemento de IPSA aplica la funcion.

# es mas simple que realizar un for

# aqui se crea una lista con cada accion

data_list <- map(ipsa, download_data)

#tenemos un problema Houston puesto que la columa Date es un caracter y debemos pasarla a fecha 

#con la funcion ymd

for (h in 1:length(ipsa)) {

  data_list[[h]]$Date <- ymd(data_list[[h]]$Date)

}

#en este codigo convertimos la data en serie de tiempo, puesto que asi es mas simple gráficar en highcharts

for (h in 1:length(ipsa)) {

  data_list[[h]] <- tk_xts(data_list[[h]],data_var = data_list[[h]]$Date)

}

# eliminamos columnas en cada data frame de cada accion

# en este caso nos quedamos con el precio ajustado

# si quieres hacer un análisis con el precio de cierre debes utilizar esto

# data_list <- sapply(data_list, FUN = function(x) subset(x, select=-c(1,2,3,5,6)), simplify = FALSE)

# dado que estamos realizando un analisis de largo plazo, se debe trabajar con el precio ajustado

# si quieres analizar un escenario de corto plazo, debes trabajar con el precio de cierre

# de esta forma, rescatas el comportamiento del precio sin ajuste de dividendos

data_list <- sapply(data_list, FUN = function(x) subset(x, select=-c(1,2,3,4,5)), simplify = FALSE)

# fusionamos la lista en un dataframe llamado all_data

all_data <- do.call(cbind, data_list)

#Reasignamos nombres

names(all_data)<- ipsa

# obtenemos retornos simples (Xt/xt-1) omitiendo los NA que se produce al perder un dato

retornos <- Return.calculate(all_data, method = "simple")%>%na.omit()

################### Simulacion de portafolios

################### simularemos 1500 portafolios

numportfolio <- 1500

#generamos una lista de dimension de 1500 para recoger cada portafolio

matriz_port <- list(1:numportfolio)

#Creamos una funcion para realizar la simulacion

portfolio_simulacion <- function (x) 

{

  # creamos listas para recibir la data que se generara

  rets <- list(1)

  vols <- list(1)

  wts <- list(1)

  # lo que esta dentro del for se repitira 1500 veces

  for (i in 1: numportfolio) {

    #Se genera numeros aleatorios en funcion de la cantidad de acciones

    weights_port <- runif(ncol(x))

    #Se genera ponderadores que sumen 1

    weights_port <- weights_port/sum(weights_port)

    # el retorno esperado del portafolio es igual a las ponderaciones x retorno de cada accion

    # Sacamos retorno anualizado y multiplicamos por 250, ya que el retorno es diario

    # 250 serían los días de transacción promedio en un año, de esta forma anualizamos retornos y volatilidad

    # para llegar a 250 se debe considerar dias de transaccion de lunes a viernes, restando feriados

    # el retorno anualizado es simplemente 250 x retorno diario

    # la volatilidad anualizada es la desviacion estandar del retorno diario multiplicada por raiz de 250

    rets<- t(round(colMeans(x)*250, digits = 2))%*%weights_port

    #La volatilidad se obtiene de la matriz de covarianza por los ponderadores

    # se multiplica por raiz de 250, de esta forma obtenemos la volatilidad anual de cada portafolio 

    vols<- sqrt(t(weights_port)%*%cov(x)%*%weights_port)*sqrt(250)

    # Se asignan los ponderadores a la variable wts

    wts <- weights_port

    # se traspone y se conviere en un data frame

    wts <- as.data.frame(t(wts))

    names(wts) <- ipsa

    # se genera un data frame con el retorno y la volatilidad

    data_port <- data.frame(rets, vols)

    # se une en un data frame data_port y wts

    data_port <- cbind(data_port,wts)

    # cada repeticion en donde se obtiene data_port se asigna a matriz_port 1500 veces

    matriz_port[[i]]<- data_port

  }

  # la funcion portfolio_simulacion entrega una lista de 1500 datos con ponderadores 

  # retornos y volatilidad

  return(matriz_port)

}

# aplicamos la funcion a "retornos" y obtenemos matriz_port y se asigna a temp

temp <- portfolio_simulacion(retornos)

# unimos cada lista en un data frame

df_unido <- do.call(rbind, temp)

# calculamos ratio de sharpe y lo añaidmos a df_unido,considerando tasa de libre riesgo 0%

df_unido$sharp_r <- df_unido$rets/df_unido$vols

# dato de interes, cual es la cartera con mayor ratio de sharpe

ubicacion <- which(df_unido$sharp_r == max(df_unido$sharp_r))

# Extraemos los podneradores de la cartera optima

Cartera_optima <- df_unido[ubicacion,]

#Graficamos el resultado

# multiplicamos x 100 para dejar los valores en %

# se grafica en funcion del ratio de sharpe

# mientras mas grande el ratio se obtiene un color mas violeta

ggplot(df_unido, aes(x = vols*100, y = rets*100, color = sharp_r)) +

  geom_point() +

  scale_color_gradientn(colors = c("cyan", "violet"), 

                        limits = c(min(df_unido$sharp_r, na.rm = TRUE), 

                                   max(df_unido$sharp_r, na.rm = TRUE)),

                        name = "Sharp Ratio") +

  theme_minimal()+

  labs(x = "Volatilidad Anualizada (%)", y = "Retorno Anualizado(%)") 

###############Proceso de optimización################

###############Frontera eficiente################

# generamos una funcion para calcular retornos y volatilidad de un portafolio

# aqui obtenemos el ratio de sharpe y se coloca el signo negativo para minimizar

# al colocar el signo negativo estariamos minimizando la funcion pero en realidad es el maximo

portfolio_opt <- function(weights_port)

{

  #weights_port <- matrix(weights_port)

  p_rets<- t(round(colMeans(retornos)*250, digits = 2))%*%weights_port

  p_vols<- sqrt(t(weights_port)%*%cov(retornos)%*%weights_port)*sqrt(250)

  return(-p_rets/p_vols)

}

# Ponderadores iniciales

int_w <- rep(1, ncol(retornos))/ncol(retornos)

# restricciones para cada accion, que va de 0 a 100% en ponderacion

# limite inferior de 0, es decir, no permitimos irnos corto en acciones o ponderadores negativos

lower_bounds <- rep(0,ncol(retornos))

# limite superior de 1 o 100%, es decir, una accion puede tener el 100%

upper_bounds <- rep(1,ncol(retornos))

# restricciones

constraint_function <- function(weights_port)

{

  #suma de ponderadores igual a 100%

  ponderador <-(sum(weights_port)-1)

  #ponderadores x retornos debe ser igual al retorno objetivo

  r_target <- (sum(t(weights_port) %*% colMeans(retornos)*250) - target_return)

  return(c(ponderador, r_target))

}

# generamos una lista de 30 para recoger la data que se generara

frontera_ef <- list(1:30)

# seteamos las rentabilidad entre 2% y 25%, generando 30 datos entre ellos, incluyendolos

# importante esta secuencia puesto que si cambias de activos, debes ver si dichas rentabilidades 

# se ajustan a los resultados

target <-seq(0.02,0.25, length= 30)

# como tenemos una lista de 30, realizamo un "for" para 30

# recogemos el retorno objetivo

# minimizamos el ratio de sharpe que dado que es negativo se maximiza

# del resultado sacamos el portafolio optimo y lo guardamos en mtx

# vola rescata la volatilidad del portafolio optimo anualziado

# generamos la frontera eficiente rescatando el retorno, la volatilidad y el ratio de sharpe

for (i in 1:30) {

  target_return <- target[i]

  #modelo para minimizar mediante el metodo SLSQP  

  opt_sharpe_sls <- slsqp(x0 = int_w,fn=portfolio_opt, lower = lower_bounds, upper=upper_bounds,

                          heq =constraint_function, control = list(xtol_rel = 1e-4, check_derivatives = TRUE))

  mtx <- as.matrix(opt_sharpe_sls[["par"]])

  vola <-   sqrt(t(mtx)%*%cov(retornos)%*%mtx)*sqrt(250)

  frontera_ef[[i]] <- c(target[i],vola, target[i]/vola, t(mtx))

}

#unimos la lista en un solo data frame

frontera_unido <- as.data.frame(do.call(rbind, frontera_ef) )

names(frontera_unido) <- c("retorno","volatilidad","sharpe",t(ipsa))

#grafico de portafolios simulados

ggplot(frontera_unido, aes(x = volatilidad*100, y = retorno*100, color = sharpe)) +

  geom_point() +

  scale_color_gradientn(colors = c("violet", "yellow"), 

                        limits = c(min(frontera_unido$sharpe, na.rm = TRUE), 

                                   max(frontera_unido$sharpe, na.rm = TRUE)),

                        name = "Sharp Ratio") +

  theme_minimal()+

  labs(x = "Volatilidad Anualizada (%)", y = "Retorno Anualizado(%)") 

#Cartera max ratio de sharpe

location <- which(frontera_unido$sharpe == max(frontera_unido$sharpe))

# Extraemos los podneradores de la cartera optima

Cartera_max_sharpe <- frontera_unido[location,]

#Cartera minima varianza

location_v <- which(frontera_unido$volatilidad == min(frontera_unido$volatilidad))

Cartera_min_var <- frontera_unido[location_v,]

for (i in 1:30) {

  target_return <- target[i]

  opt_sharpe <- nloptr(x0 = int_w, eval_f = portfolio_opt, lb =lower_bounds,

                       ub = upper_bounds, opts = options, eval_g_eq = constraint_function)

  mtx <- as.matrix(opt_sharpe[["solution"]])

  vola <-   sqrt(t(mtx)%*%cov(retornos)%*%mtx)*sqrt(250)

  frontera_ef[[i]] <- c(target[i],vola, target[i]/vola)

}

frontera_grap <- frontera_unido[,c(-1,-2)]

retornos_largo <- frontera_grap %>%

  pivot_longer(cols = -sharpe, names_to = "ACCIONES", values_to = "Ponderacion")

#graficamos cada portafolio con su ratio de sharpe

ggplot(retornos_largo, aes(x = sharpe, y = Ponderacion, fill = ACCIONES)) +

  geom_bar(stat = "identity", width= 0.01) +

  labs(title = "Ponderación de Acciones", x = "Sharp Ratio", y = "Ponderación") +

  theme_minimal()

#gráficamos los precios de cada accion normalizado en 0

names(data_list)<- ipsa

hchart_obj<- highchart(type = "stock") %>%

  hc_title(text = "Gráfico comparativo") %>%

  hc_xAxis((type="datatime"))%>%

  hc_tooltip(shared =T, crosshairs=T)%>%

  hc_add_theme(hc_theme_db())

for (i in seq_along(data_list)) {

  hchart_obj <- hchart_obj %>%

    hc_add_series(data_list[[i]], name = names(data_list)[i], compare="percent")

}

hchart_obj

Esperamos que haya sido útil.

Cualquier duda, no dudes en contactarnos.

Saludos

Equipo Principia…