Anexo. 4 Códigos

A continuación, se presenta el código utilizado para la realización del procedimiento descrito en el desarrollo del trabajo.

Datos

Obtención de Datos

Lo primero que se realizó fue cargar la tabla de las empresas.

empresas <- read_excel("data/000_empresas.xlsx")

Luego se extrajeron los ticks de las empresas.

ticks <- empresas |> 
  select(TICKERS) |> 
  pull()

Una vez almacenados los ticks de las empresas en la variable ticks se procedió a descargar los datos correspondientes a dichas empresas desde Yahoo Finance usando el paquete quantmod de Ryan y Ulrich (2023).

nombres_colum <- c("Date","Open","High","Low","Close","Volume","Adjusted")
qmd_data <- list()
for (i in 1:length(ticks)) {
  tick <- ticks[i]
  value <- getSymbols(
    tick,
    from = "2000-01-02",
    to = "2023-03-01",
    auto.assign = F,
    periodicity = "monthly") |>
    as.data.frame()
  dates <- row.names(value)
  row.names(value) <- NULL
  value <- cbind(dates,value)
  names(value) <- nombres_colum
  qmd_data[[tick]] <-  value
}

Con el objetivo de realizar un análisis exploratorio de los datos, se decidió realizar una evaluación visual de los datos históricos del precio ajustado para lo que se ejecutó:

lapply(qmd_data, function(x){
  x |>
    ggplot(aes(x=as.Date(Date), y=Adjusted))+
             geom_line(color="#065AD8")
})

Tras el análisis visual ejecutado con el fragmento de código anterior se persibió la existencia de precios constantes, así como cálculos erróneos en el precio ajustado correspondiente a los primeros años de algunas series. Con el objetivo de eliminar estas irregularidades se seleccionaron solo aquellas observaciones posteriores a enero del 2005.

returns_emps <- qmd_data |>
  lapply(function(x){
    emps <- x |>
      filter(Date >= "2005-01-31")
  })

Con el objetivo de determinar si los datos que habían sido importados contaban con valores faltantes se ejecutó el siguiente código:

na_values <- returns_emps |>
  sapply(function(x){
    na <- length(which(is.na(x)))
  })
emp_con_na <- which(na_values > 0)

Con el objetivo de solucionar el problema con respecto al incorrecto registro de los datos se decidió eliminar aquellas que no presentaran variaciones en los precios en más de 10 observaciones. Para lo que primero se computaron las rentabilidades ejecutando el siguiente código, mediante el cual además se eliminaron las series con valores faltantes.

returns_emps2 <- returns_emps[-emp_con_na] |>
  lapply(function(x){
    returns <- x |>
      select(Date, Adjusted) |>
      mutate(Return_Ad = Delt(Adjusted)[,1]) |>
      na.omit() |>
      select(Date, Return_Ad)
  })

Una vez computadas las rentabilidades se eliminaron aquellas series que presentaban en más de 10 observaciones rentabilidad 0, para lo que se ejecutó el siguiente código.

zero_values <- returns_emps2 |>
  sapply(function(x){
    zeros <- length(which(x[,2]==0))
  })
returns_emps3 <- returns_emps2[zero_values<10]

Indicadores

A continuación, se expone el código utilizado durante el proceso expuesto en el sub-epígrafe indicadores del capítulo 2.

Primero se descargaron los datos del IBEX, se computaron las rentabilidades del precio ajustado del mismo y se seleccionaron los valores posteriores a enero del 2005.

#Importando IBEX
IBEXsel <- getSymbols(
  "^IBEX",
  from = "1990-01-01",
  to = "2023-03-01",
  auto.assign = F,
  periodicity = "monthly") |>
  as.data.frame()
dates <- row.names(IBEXsel)
row.names(IBEXsel) <- NULL
IBEXsel <- cbind(dates,IBEXsel)
names(IBEXsel) <- nombres_colum
# Calculando rentabilidad y seleccionando observaciones posteriores a
# enero del 2005.
IBEXsel <- IBEXsel |>
  mutate(Return_I = Delt(Adjusted)[,1]) |>
  na.omit() |>
  filter(Date >= "2005-01-31") |>
  select(Date, Return_I)

Luego se agregaron los valores de las rentabilidades del IBEX a las tablas de las rentabilidades de las acciones de las empresas seleccionadas, y se computaron y agregaron las variables listadas a continuación a cada una de las tablas:

• Volatilidad de la empresa
• Volatilidad del índice
• Correlación entre las rentabilidades de la empresa y el indice
• La Beta entre la empresa y el índice

returns_indc <- returns_emps3 |>
  lapply(function(x, ind = IBEXsel){
    emp <- x |>
      left_join(ind) |>
      mutate(
        VE = sqrt(cumsum((Return_Ad - cummean(Return_Ad))^2)/1:length(Return_Ad)),
        VI = sqrt(cumsum((Return_I - cummean(Return_I))^2)/1:length(Return_I)),
        Cor = cumsum((Return_Ad-cummean(Return_Ad))*(Return_I-cummean(Return_I)))/(sqrt(cumsum((Return_Ad-cummean(Return_Ad))^2))*sqrt(cumsum((Return_I-cummean(Return_I))^2)))
      )|>
      na.omit() |>
      mutate(
        Beta = (Cor*VE)/VI
      )
  })

Vectores

A continuación, se expone el código utilizado durante el proceso expuesto en el sub-epígrafe vectores del epígrafe modelado del Capítulo 2.

El primer paso llevado a cabo para la ejecución del proceso explicado en el sub-epígrafe en cuestión fue crear una función que permitió obtener las muestras consecutivas para cada serie utilizada. La función expuesta a continuación, como ya se mencionó, permite obtener las muestras consecutivas de una serie, para lo que se utilizan los parámetros mencionados en el sub-epígrafe, número de observaciones de entradas y número de observaciones de salida, así como un parámetro condicional con el que se indica si el vector a crear es de entrada o de salida.

vector2dmaker <- function(vec, ent, sal, eos=T){
  if(eos==T){
    emp <- 1
    term <- (length(vec) - (ent+sal-1))
    ob <- ent
  }else{
    emp <- ent + 1
    term <- (length(vec)-sal+1)
    ob <- sal
  }
  
  vec2d <- sapply(emp:term,
               function(x) vec[x:(x + ob-1)]) |>
    matrix(nrow = ob) |>
    t()
  
  return(vec2d)
}

A continuación, se muestra el código utilizado para la creación de los vectores de entrada de correspondiente a cada una de las series. Para lo que primero se crearon dos funciones una para las entradas y otra para las salidas.

# Función que se utlizará para crear las entradas tridimensionales
input3dmaker <- function(x,inp,out){
  empre <- x
  series <- 2:dim(x)[2]
  for (i in series) {
    if(i==series[1]){
      vec3d <- vector2dmaker(empre[[i]],ent=inp,sal=out)
    }else{
      vec3d <- abind(vec3d,vector2dmaker(empre[[i]],ent=inp,sal=out), along = 3)
    }
  }
  return(vec3d)
}

# Función que se utlizará para crear las salidas tridimensionales
output3dmaker <- function(x,inp,out){
  empre <- x[["Return_Ad"]]
  vec3d <- vector2dmaker(empre,ent=inp,sal=out,F)
  dim(vec3d) <- c(dim(vec3d),1)
  return(vec3d)
}

Luego se crearon las listas de vectores tridimensionales de entradas y salidas por empresa, ejecutándose el siguiente código otras dos veces con el objetivo de crear las listas vecs3d2e y vecs3d3e que corresponden a aquellos casos en los que se seleccionaron 2 y 3 entradas.

#Se define el horizonte temporal
ht <- 1

#Se definen las observaciones de entrada
oe <- 1

#Se crean los vectores de entrada 3d y el vector 2d de salida para tamaño de entrada 1
vecs3d1e <- list()
for(i in 1:length(returns_indc)){
  emp <- returns_indc[[i]]
  inps <- input3dmaker(emp, oe, ht)
  outs <- output3dmaker(emp, oe, ht)
  dates <- emp[(oe + ht):dim(emp)[1],1]
  id <- rep(names(returns_indc)[i],length(dates))
  tibblex <- tibble(
    Date = dates,
    ID = id,
    inputs = inps,
    outputs = outs
  )
  vecs3d1e[[names(returns_indc)[i]]] <- tibblex
}

Modelado y entrenamiento

A continuación, se presenta el código utilizado durante el proceso descrito en los distintos sub-epígrafes del epígrafe Modelado y entrenamiento.

Modelado

Para la creación de los modelos el primer paso a ejecutar es obtener la información de los vectores para los que se va a construir el modelo, lo que se hizo ejecutando el siguiente código:

data <- bind_rows(vecs3d1e)
data <- data  |>
  arrange(Date)
inputsinfo <- data|>
  select(inputs) |>
  pull() |>
  dim()
outputsinfo <- data|>
  select(outputs) |>
  pull() |>
  dim()

# Definir parámetros
n_ob_pas <- inputsinfo[2]
n_variables <- inputsinfo[3]
n_ob_fut <- outputsinfo[2]

Luego se constituyó la estructura de los modelos con los aspectos descritos en 2.5.1 Modelado.

# Capa de entrada
inp <- layer_input(
  shape = c(NULL,n_ob_pas,n_variables))

# Capas ocultas
# - CNN
cnn <- inp |>
  layer_conv_1d(
    filters = 64,
    kernel_size = 1,
    activation = layer_activation_leaky_relu())
# - LSTM
lstm <- cnn |>
  layer_lstm(64)

# Capa de Salida
out <- lstm |> 
  layer_dense(
    n_ob_fut*1)

# Juntar las capas para constituir el modelo 
model <- keras_model(inp, out)
# Estableciendo parámetros de aprendizaje
model |> 
  compile(loss = "mse", optimizer = optimizer_sgd(0.0005))

Nota

Puede encontrar modelos sin entrenar en la carpeta data del repositorio en el que se encuentra el presente trabajo. Los modelos fueron guardados usando la extensión hdf5 y bajo los nombres model1e, model2e y model3e.

Entrenamiento

El primer paso es definir la función a utilizar para el entrenamiento de los modelos. Esta función fue construida con el objetivo de emplear el método de entrenamiento descrito en 2.5.2 Entrenamiento. Como resultado esta función devuelve una lista que contiene las predicciones obtenidas y el modelo después de haber sido entrenado y tomará como entradas principales el tibble llamado data constituido en el primer paso que se expone en este Anexo en la sección Modelado y el modelo además de otros argumentos que permite la utilización de la función con unos inputs principales que no sean los utilizados en el presente trabajo.

wfv_train <- function(x, modelo, seq_var_name, inp_var_name = "inputs", out_var_name = "outputs", progress_bar=T){
  
  predictions <- c()
  seq_val <- unique(x[[seq_var_name]])
  
  if(progress_bar){
    pb <- txtProgressBar(min = 0, max = length(seq_val), initial = 0, style = 3)
  }
  
  
  # Iteración que se ejecutará para cada valor único en la variable que define la secuencia de los datos. Por ello es de vital importancia que los datos en el tibble x se encuentren ordenados por la variable de secuencia cuyo nombre se pasa a seq_var_name
  
  for (i in 1:length(seq_val)) {
    val_seq <- seq_val[i]
    #Extraer entradas y salidas correspondiente al periodo en la variable de secuencia actual
    inputs <- x |>
      filter(!!sym(seq_var_name) == val_seq) |>
      select(!!sym(inp_var_name)) |>
      pull()
    outputs <- x |>
      filter(!!sym(seq_var_name) == val_seq) |>
      select(!!sym(out_var_name)) |>
      pull()
    outputs <- outputs[,,1]
    
    #Usar entradas para obtener predicciones para los periodos en la variable secuencia a excepción del primero
    if(i > 1){
      pred <- modelo |>
        predict(inputs, verbose = 3)
      predictions <- rbind(predictions, pred)
    }
    
    # Entrenar el modelo
    modelo |>
      fit(
        inputs,
        outputs,
        epochs = 1,
        batch_size = 10,
        shuffle = F,
        verbose = 0)
    
    if(progress_bar){
      setTxtProgressBar(pb,i)
      }
    
  }
  
  if(progress_bar){
    close(pb)
  }
  
  results <- list()
  results[['predicciones']] <- predictions
  results[['modelo']] <- modelo
  return(results)
}

Una vez creada la función se obtuvieron las predicciones utilizando el siguiente código:

resultados <- wfv_train(data,model,'Date')
predicciones1e <- resultados$predicciones

Nota

Puede encontrar modelos entrenados en la carpeta data del repositorio en el que se encuentra el presente trabajo. Los modelos fueron guardados usando la extensión hdf5 y bajo los nombres model1etd, model2etd y model3etd.

Como se explica en 2.6.1 Predicciones además de las predicciones obtenidas por los modelos se computaron predicciones obtenidas a partir del uso de la media aritmética, para comparar con las obtenidas con los modelos. Para el computo de estas predicciones se creó la siguiente función:

wfv_means <- function(x, seq_var_name, inp_var_name = "inputs", out_var_name = "outputs", id_var_name, progress_bar=T){
  
  means <- c()
  seq_val <- unique(x[[seq_var_name]])
  
  if(progress_bar){
    pb <- txtProgressBar(min = 0, max = length(seq_val), initial = 0, style = 3)
  }
  
  for (i in 1:length(seq_val)) {
    val_seq <- seq_val[i]
    inputs <- x |>
      filter(!!sym(seq_var_name) == val_seq) |>
      select(!!sym(inp_var_name)) |>
      pull()
    inputspred <- x |>
      filter(!!sym(seq_var_name) == val_seq) |>
      select(!!sym(inp_var_name)) |>
      pull()
    outputs <- x |>
      filter(!!sym(seq_var_name) == val_seq) |>
      select(!!sym(out_var_name)) |>
      pull()
    outputs <- outputs[,,1]
    
    ids <- x |>
      filter(!!sym(seq_var_name) == val_seq) |>
      select(!!sym(id_var_name)) |>
      pull()
    
    if(i==1){
      dfmeans <- inputs[,,1] |>
        as.data.frame() |>
        cbind(ID = ids)
    }else{
      dfmeansupd <- inputs[,dim(inputs)[2],1] |>
        as.data.frame() |>
        cbind(ID = ids)
      names(dfmeansupd)[1] <- paste0("V",(dim(dfmeans)[2]))
      idsdf <- unique(c(ids, dfmeans[[id_var_name]]))
      idsdf <- data.frame(ID = idsdf)
      dfmeansupd <- dplyr::left_join(idsdf, dfmeansupd, by = "ID")
      ifelse(
        dim(dfmeansupd)[1] > dim(dfmeans)[1],
        dfmeans <- dplyr::left_join(dfmeansupd, dfmeans, by = "ID"),
        dfmeans <- dplyr::left_join(dfmeans, dfmeansupd, by = "ID")
        )
    }
    
    if(i > 1){
      MEANS <-  dfmeans |>
        rowwise() |>
        mutate(
          means = mean(c_across(-!!sym(id_var_name)), na.rm = T)) |>
        slice(match(ids,!!sym(id_var_name))) |>
        pull(means) |>
        as.matrix()
      means <- rbind(means, MEANS)
    }
    
    if(progress_bar){
      setTxtProgressBar(pb,i)
    }
    
  }
  
  if(progress_bar){
    close(pb)
  }
  
  return(means)
}

Una vez creada la función se obtuvieron las predicciones utilizando el siguiente código:

meanse1 <- wfv_means(data,'Date',id_var_name = "ID")

Nota

En adición a lo expuesto con anterioridad se crearon dos funciones getconfig y plot_modelk, en el archivo .Rprofile del repositorio en el que se encuentra este trabajo, que permiten graficar la estructura de los modelos mediante el uso del paquete Iannone (2023), como se ve en la Figura 17. El código a utilizar sería:

# Las funciones están creadas para graficar las estructuras utilizadas en el presente trabajo.
model |>
  getconfig() |>
  plot_modelk() |>
  grViz()

El procedimiento expuesto en las secciones Modelado y Entrenamiento del presente anexo fue repetido para construir los 10 modelos realizados a partir de cada grupo de vectores tridimensionales, sustituyendo en el primer código expuesto la llamada a vecs3d1e por vecs3d2e y vecs3d3e,según el grupo de vectores tridimensionales utilizado.

Resultado

A continuación, se presenta el código utilizado durante el proceso descrito en los distintos sub-epígrafes del epígrafe Resultado.

Predicciones

El análisis expuesto en 2.6.1 Predicciones fue realizado a partir de gráficas (vea Figura 8, Figura 9 y Figura 10), en las que se recogen los valores de los indicadores \(MSE\) y \(R^2\) para cada una de las estructuras probadas.

El primer paso para la obtención de estas gráficas fue el de computar los indicadores, para cada periodo de tiempo, para cada una de las predicciones obtenidas a partir de los distintos modelos construidos con cada estructura. Esto se realizó mediante el siguiente código.

#Extraer salidas reales
salidas <- data |>
  filter(
    Date > data$Date[1]
  ) |>
  select(outputs) |>
  pull()
salidas <- salidas[,,1]

#Computar indicadores MSE y R2
indicadores <- data |>
  filter(Date > data$Date[1]) |>
  cbind(predicciones = predicciones1e[,1]) |>
  cbind(means = meanse1) |>
  mutate(salidas = salidas) |>
  select(Date, predicciones, means, salidas) |>
  group_by(Date) |>
  summarise(
    r2 = 1 - (sum((salidas - predicciones)^2)/sum((salidas - means)^2)),
    mse = mse(predicciones, salidas),
  )

Los diferentes indicadores computados para cada uno de las 10 modelos entrenados con cada una de las estructuras fueron guardados en una lista llamada list_indicadores. Esto se realizo utilizando el siguiente código:

list_indicadores[["indicadores1"]] <-  indicadores

Una vez realizado esto se obtiene una lista que contiene 10 data frames (indicadores1,…,indicadores10), los cuales a su vez contienen los valores de los de \(MSE\) y \(R^2\) de las predicciones obtenidas por los modelos de RNA para cada una de las empresas agrupados por fecha. Por lo que luego se construyó la gráfica mediante el uso del siguiente código.

# Agrupar la información de las distintas construcciones en un solo data frame
indi_graf_data <- do.call(cbind,list_indicadores)

# Obtener los resultados medios, para cada periodo de tiempo, usando las distintas construcciones
indi_graf_data |>
  rowwise() |>
  mutate(
    Date = `indicadores1.Date`,
    meanmse = mean(c_across(contains("mse"))),
    meanr2 = mean(c_across(contains("r2")))
    ) |>
  select(
    Date, meanmse,meanr2
  )|>
  # Graficar
  mutate(
    Date = as.Date(Date)) |>
  ggplot(aes(x = Date, group = 1)) +
  geom_line(aes(y = meanmse, color = "MSE")) +
  geom_line(aes(y = meanr2, color = "R2")) +
  scale_color_manual(values = c("blue", "green")) +
  theme(axis.text.x = element_text(angle = 90)) +
  labs(x = "Fecha", y = "Indicadores", color = "Indicadores")

Además de las gráficas se utilizó también en el análisis de los resultados la Tabla 4, en la que se encuentran las empresas que obtuvieron mejores y peores indicadores para cada estructura, para la obtención de estos datos se uso el siguiente código:

indicadores_X_emp <- data |>
  filter(Date > data$Date[1]) |>
  cbind(predicciones = predicciones1e[,1]) |>
  cbind(means = meanse1) |>
  mutate(salidas = salidas) |>
  select(Date, predicciones, means, salidas, ID) |>
  group_by(ID) |>
  summarise(
    r2 = 1 - (sum((salidas - predicciones)^2)/sum((salidas - means)^2)),
    mse = mse(predicciones, salidas)
  ) |>
  select(ID, r2, mse)

Al igual que los indicadores computados por fecha para guardar los indicadores computados por empresa se creó una lista denominada list_indic_emp. Luego de haber almacenado los 10 data frames de indicadores por empresa en la lista se extrajeron las empresas con los mejores y peores resultados mediante el siguiente código:

# Agrupar la información de las distintas construcciones en un solo data frame
ind_emp_t <- do.call(rbind, list_indic_emp)

# Computar los R2 y MSE medios por empresa
ind_emp_t <- ind_emp_t |>
  group_by(ID) |>
  summarize(
    r2 = mean(r2),
    mse= mean(mse)) |>
  ungroup() |>
  arrange(desc(r2))

# Obtener las 10 empresas con mejores y peores indicadores
mejores10 <- head(ind_emp_t,10)
peores10 <- tail(ind_emp_t,10)

Y mediante la utilización de las variables anteriores y el uso de las funciones rbind() y cbind fue como se creó la Tabla 4.

Composición de carteras

En esta sección se explica como se realizó el análisis del comportamiento de los resultados obtenidos por las distintas carteras (ver Figura 11, Figura 12 y Figura 13). Para ello es necesario primero obtener la composición de las carteras, por fecha, a partir de las predicciones obtenidas mediante el uso de las medias aritméticas y los modelos de RNA.

Para el cálculo de la composición de las carteras se usó el paquete de R llamado quadprog (Berwin A. Turlach R port by Andreas Weingessel <Andreas.Weingessel@ci.tuwien.ac.at> Fortran contributions from Cleve Moler dpodi/LINPACK) (2019)), a continuación, se muestra el código utilizado para hallar la composición de carteras a partir de las predicciones de la media:

# Se creo un data frame en el que se guardó toda la información:
#   - Valores del IBEX, como índice de referencia
#   - Valores de las predicciones, tanto las obtenidas por el modelo de RNA como por las medias aritméticas

DATA <- data |>
  left_join(IBEXsel, by ="Date") |>
  mutate(IBEX = Return_I) |>
  arrange(Date) |>
  filter(
    Date > data$Date[1]
  ) |>
  mutate(predicciones = predicciones1e[,1]) |>
  mutate(
    Real = salidas,
    RNA = predicciones,
    Means = meanse1
  ) |>
  select(Date, Real, IBEX, RNA, Means, ID)

# A partir del data frame DATA se crearon:
#    - Un data frame cuyas columnas son los datos reales de cada una de las empresas para cada uno de los periodos de tiempo para los que se obtuvieron predicciones.
#    - Un data frame cuyas columnas son los datos obtenidos mediante el uso de las medias aritméticas de cada una de las empresas para cada uno de los periodos de tiempo para los que se obtuvieron predicciones.

pvtReal <- DATA |>
  select(Date, Real, ID) |>
  pivot_wider(
    names_from = ID,
    values_from = Real
  )

pvtMeans <- DATA |>
  select(Date, Means, ID) |>
  pivot_wider(
    names_from = ID,
    values_from = Means
  )

# Se creó el data frame en el que se guardará la composición de las carteras para cada uno de los periodos para los que se obtuvó predicción
weightsm <- data.frame()

# Iteración mediante la cual se halla la composición de las carteras

pb <- txtProgressBar(min = 0, max = length(unique(data$Date)[-1]), initial = 0, style = 3)

for (i in 1:length(unique(data$Date)[-1])) {
  if(i>1){
    
    # Se crea el data frame que comprende los datos a utilizar para hallar la composición de la cartera, este está creado por los datos reales hasta la fecha y la previsión del siguiente periodo
    datamQP <- pvtReal |>
      filter(Date < unique(data$Date)[-1][i]) |>
      rbind(pvtMeans |>
              filter(Date == unique(data$Date)[-1][i])
      )
    
    # Elimina aquellas empresas que no tengan ni datos reales o de previsión
    nare <- which(is.na(datamQP[dim(datamQP)[1],]))
    naremo <- which(is.na(datamQP[(dim(datamQP)[1]-1),]))
    nare <- c(nare,naremo)
    nare <- unique(nare)
    if(length(nare) != 0){
      carteram <- datamQP[, - nare]
    }else{
      carteram <- datamQP
    }
    
    # Extrae las previsiones
    returnm <- carteram[dim(carteram)[1], -1] |>
      as.matrix() |>
      t()
    
    # Calcula la matriz de covarianza
    covmm <- cov(carteram[, -1], use = "complete.obs")
    npcovmm <- nearPD(covmm)$mat |> 
      as.matrix()
    # Extrae el número de empresas
    n <- ncol(npcovmm)
    
    # Halla la composición de la cartera
    qp_outm <- solve.QP(
      Dmat = 2*npcovmm,
      dvec = rep(0,n),
      Amat = cbind(-1, diag(n)),
      bvec = c(-1, rep(0,n)),
      meq = 1)
    qp_outm <- qp_outm$solution
    qp_outm <- floor(qp_outm*100)/100
    for(j in 1:length(qp_outm)){
      if(qp_outm[j] < 0.001){
        qp_outm[j] <- 0
      }else{}
    }
    
    # Guarda la composición de la cartera
    names(qp_outm) <- names(carteram[, -1])
    weightsm <- bind_rows(weightsm, qp_outm)
  }
  
  setTxtProgressBar(pb,i)
}

close(pb)

# Sustituir los pesos y observaciones reales con valores faltantes con cero
pvtReal[is.na(pvtReal)] <- 0
weightsm[is.na(weightsm)] <- 0

Luego para hallar la rentabilidad de la cartera se multiplicaron las composiciones por las rentabilidades reales, se asumió que se invertía uno en el primer periodo y se realizó una sumatoria acumulativa a los largo de los valores para obtener el comportamiento de la rentabilidad a lo largo del tiempo.

# Hallando las rentabilidades de las carteras conformadas a partir de las predicciones de la media aritmética

return_CM <-  weightsm * pvtReal[-1,-1]
return_CM <- rowSums(return_CM)
return_CM <- c(1,return_CM)
return_CM <- data.frame(
  Date = pvtReal[,1],
  Mean = return_CM
)

Los mismos pasos que se realizaron para hallar el comportamiento de la rentabilidad de las carteras a partir de las medias aritméticas se realizaron para hallar el comportamiento a partir de las predicciones obtenidas por el modelo de RNA como se ve en el código a continuación.

# A partir del data frame DATA se creó un data frame cuyas columnas son los datos obtenidos mediante el uso del model de RNA de cada una de las empresas para cada uno de los periodos de tiempo para los que se obtuvieron predicciones.

pvtRNA <- DATA |>
  select(Date, RNA, ID) |>
  pivot_wider(
    names_from = ID,
    values_from = RNA
  )

# Se creó el data frame en el que se guardara la composición de las carteras para cada uno de los periodos para los que se obtuvó predicción
weightse <- data.frame()

# Iteración mediante la cual se halla la composición de las carteras

pb <- txtProgressBar(min = 0, max = length(unique(data$Date)[-1]), initial = 0, style = 3)

for (i in 1:length(unique(data$Date)[-1])) {
  if(i>1){
    # Se crea el data frame que comprende los datos a utilizar para hallar la composición de la cartera, este estó creado por los datos reales hasta la fecha y la previsión del siguiente periodo
    dataeQP <- pvtReal |>
      filter(Date < unique(data$Date)[-1][i]) |>
      rbind(pvtRNA |>
              filter(Date == unique(data$Date)[-1][-1][i])
            )
    # Elimina aquellas empresas que no tengan ni datos reales o de previsión
    nare <- which(is.na(dataeQP[dim(dataeQP)[1],]))
    naremo <- which(is.na(dataeQP[(dim(dataeQP)[1]-1),]))
    nare <- c(nare,naremo)
    nare <- unique(nare)
    if(length(nare) != 0){
      carterae <- dataeQP[, - nare]
    }else{
      carterae <- dataeQP
    }
    
    # Extrae las previsiones
    returne <- carterae[dim(carterae)[1], -1] |>
      as.matrix() |>
      t()
    
    # Calcula la matriz de covarianza
    covme <- cov(carterae[, -1], use = "complete.obs")
    npcovme <- nearPD(covme)$mat |> 
      as.matrix()
    # Extrae el número de empresas
    n <- ncol(npcovme)
    
    # Halla la composición de la cartera
    qp_oute <- solve.QP(
      Dmat = 2*npcovme,
      dvec = rep(0,n),
      Amat = cbind(-1, diag(n)),
      bvec = c(-1, rep(0,n)),
      meq = 1)
    qp_oute <- qp_oute$solution
    qp_oute <- floor(qp_oute*100)/100
    for(j in 1:length(qp_oute)){
      if(qp_oute[j] < 0.001){
        qp_oute[j] <- 0
      }else{}
    }
    
    # Guarda la composición de la cartera
    names(qp_oute) <- names(carterae[, -1])
    weightse <- bind_rows(weightse, qp_oute)
  }
  
  setTxtProgressBar(pb,i)
}

close(pb)

# Sustituir los pesos con valores faltantes con cero
weightse[is.na(weightse)] <- 0

Luego para hallar la rentabilidad de la cartera se multiplicaron las composiciones por las rentabilidades reales, se asumió que se invertía uno en el primer periodo y se realizó una sumatoria acumulativa a lo largo de los valores para obtener el comportamiento de la rentabilidad a lo largo del tiempo.

# Hallando las rentabilidades de las carteras conformadas a partir de las predicciones del modelo RNA

return_CRNA <-  weightse * pvtReal[-1,-1]
return_CRNA <- rowSums(return_CRNA)
return_CRNA <- c(1,return_CRNA)
return_CRNA <- data.frame(
  Date = pvtReal[,1],
  RNA = return_CRNA
)

Luego, al igual que con los indicadores se creó una lista list_ret_RNA en la que se guardaron los data frames de los distintos modelos construidos con cada una de las estructuras. Después se ejecutó el siguiente código para obtener la gráfica.

# Hallando el comportamiento de las rentabilidades del IBEX para el periodo

IBEXvals <- IBEXsel |>
    filter(Date > unique(data$Date)[2]) |>
    select(2) |>
    pull()
IBEXvals <- c(1, IBEXvals)

data_rent_RNA <- do.call(cbind,list_ret_RNA)
data_rent_RNA <- data_rent_RNA |>
  mutate(
    Date = RNA1.Date,
    IBEX = IBEXvals,
    Means = return_CM$Mean) |>
  mutate_at(vars(contains(".RNA")), ~ cumsum(.)) |>
  mutate(
    IBEX = cumsum(IBEX),
    Means = cumsum(Means)) |>
  group_by(Date) |>
  summarize(
    meanRNA = mean(c_across(contains(".RNA"))),
    max_y = max(c_across(contains(".RNA"))),
    min_y = min(c_across(contains(".RNA"))),
    min_5 = unname(quantile(c_across(contains(".RNA")),0.05)),
    max_95 = unname(quantile(c_across(contains(".RNA")),0.95)),
    IBEX = IBEX,
    Means = Means)
data_rent_RNA |>
  mutate(
    Date = as.Date(Date)) |>
ggplot(aes(x = Date)) +
  geom_ribbon(aes(ymin = min_y, ymax = min_5), fill = "blue", alpha = 0.3) +
  geom_ribbon(aes(ymin = max_y, ymax = max_95), fill = "blue", alpha = 0.3) +
  geom_ribbon(aes(ymin = min_5, ymax = max_95), fill = "blue", alpha = 0.6) +
  geom_line(
    aes(y = meanRNA, color = "Media RNA1"),
    linetype = "dashed") +
  geom_line(aes(y = max_y), color = "blue") +
  geom_line(aes(y = min_y), color = "blue") +
  geom_line(aes(y = max_95), color = "blue") +
  geom_line(aes(y = min_5, color = "RNA1")) +
  geom_line(aes(y = IBEX, color = "IBEX")) +
  geom_line(aes(y = Means, color = "Medias")) +
  scale_color_manual(
    values = c(
      "Media RNA1"="blue",
      "RNA1" = "blue",
      "IBEX" = "red",
      "Medias" = "green")) +
  guides(
    color = guide_legend(
      override.aes = list(
        linetype = c("solid","dashed","solid","solid"))))+
  labs(x = "Fecha",
       y = "Rentabilidades",
       color = "Leyenda")+
  theme_minimal()

Berwin A. Turlach R port by Andreas Weingessel <Andreas.Weingessel@ci.tuwien.ac.at> Fortran contributions from Cleve Moler dpodi/LINPACK), S original by. 2019. quadprog: Functions to Solve Quadratic Programming Problems. https://CRAN.R-project.org/package=quadprog.

Iannone, Richard. 2023. DiagrammeR: Graph/Network Visualization. https://CRAN.R-project.org/package=DiagrammeR.

Ryan, Jeffrey A., y Joshua M. Ulrich. 2023. quantmod: Quantitative Financial Modelling Framework. https://CRAN.R-project.org/package=quantmod.