library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.2     ✔ tibble    3.2.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.0.4     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(tidyquant)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo 
## ── Attaching core tidyquant packages ─────────────────────── tidyquant 1.0.11 ──
## ✔ PerformanceAnalytics 2.0.8      ✔ TTR                  0.24.4
## ✔ quantmod             0.4.28     ✔ xts                  0.14.1── Conflicts ────────────────────────────────────────── tidyquant_conflicts() ──
## ✖ zoo::as.Date()                 masks base::as.Date()
## ✖ zoo::as.Date.numeric()         masks base::as.Date.numeric()
## ✖ dplyr::filter()                masks stats::filter()
## ✖ xts::first()                   masks dplyr::first()
## ✖ dplyr::lag()                   masks stats::lag()
## ✖ xts::last()                    masks dplyr::last()
## ✖ PerformanceAnalytics::legend() masks graphics::legend()
## ✖ quantmod::summary()            masks base::summary()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(lubridate)
tickers <- c("AAPL","MSFT","NVDA","TSLA","NESN.SW","SIEGY","BBAJIOO.MX","FEMSAUBD.MX")
start  <- as.Date("2022-01-01")
end    <- Sys.Date()
prices <- tq_get(tickers, from = start, to = end, get = "stock.prices") %>%
  select(symbol, date, adjusted) %>%
  arrange(symbol, date)
rets <- prices %>%
  group_by(symbol) %>%
  mutate(ret_log = log(adjusted/lag(adjusted))) %>%
  ungroup() %>%
  filter(!is.na(ret_log))
Tabla <- rets %>%
  group_by(symbol) %>%
  summarise(
    media = mean(ret_log, na.rm = TRUE),
    sd    = sd(ret_log,   na.rm = TRUE),
    var   = var(ret_log,  na.rm = TRUE),
    .groups = "drop"
  ) %>%
  as.data.frame()
rownames(Tabla) <- Tabla$symbol
Tabla <- Tabla[, c("media","sd","var")]

wide <- rets %>%
  select(date, symbol, ret_log) %>%
  pivot_wider(names_from = symbol, values_from = ret_log) %>%
  arrange(date)

mat_cov <- cov(wide[, tickers], use = "pairwise.complete.obs")
activos <- c("NESN.SW","MSFT","NVDA")

EY  <- as.numeric(Tabla[activos, 1])
S   <- as.matrix(mat_cov[activos, activos])
one <- rep(1, 3)

w_gmv_unc <- solve(S, one) / as.numeric(t(one) %*% solve(S, one))
names(w_gmv_unc) <- activos
ER_gmv_unc <- sum(w_gmv_unc * EY)
VR_gmv_unc <- as.numeric(t(w_gmv_unc) %*% S %*% w_gmv_unc)
SD_gmv_unc <- sqrt(VR_gmv_unc)

paso <- 0.01
grid <- expand.grid(a = seq(0,1,by=paso), b = seq(0,1,by=paso))
grid$c <- 1 - grid$a - grid$b
grid <- grid[grid$c >= 0, ]
W   <- as.matrix(grid[, c("a","b","c")])
ER  <- as.numeric(W %*% EY)
VR  <- rowSums((W %*% S) * W)
SD  <- sqrt(VR)
grid$ER <- ER; grid$VR <- VR; grid$SD <- SD

idx_gmv_lo <- which.min(grid$VR)
w_gmv_lo   <- as.numeric(W[idx_gmv_lo, ])
names(w_gmv_lo) <- activos
ER_gmv_lo  <- grid$ER[idx_gmv_lo]
SD_gmv_lo  <- grid$SD[idx_gmv_lo]

riesgo_obj_mult <- 1.25
SD_obj <- SD_gmv_lo * riesgo_obj_mult
cand <- grid[grid$SD <= SD_obj + 1e-12, ]
idx_maxER <- which.max(cand$ER)
w_risky   <- as.numeric(W[as.integer(rownames(cand))[idx_maxER], ])
names(w_risky) <- activos
ER_risky  <- cand$ER[idx_maxER]
SD_risky  <- cand$SD[idx_maxER]

pr <- function(w) sprintf("%s: %.2f%%", names(w), 100*w)
cat("\nGMV (sin restricción no-neg):\n", paste(pr(w_gmv_unc), collapse=" | "),
    "\nE[R]=", round(ER_gmv_unc,6), " SD=", round(SD_gmv_unc,6), "\n", sep="")
## 
## GMV (sin restricción no-neg):
## NESN.SW: 71.61% | MSFT: 28.88% | NVDA: -0.49%
## E[R]=-0.000209 SD=0.009306
cat("\nGMV long-only:\n", paste(pr(w_gmv_lo), collapse=" | "),
    "\nE[R]=", round(ER_gmv_lo,6), " SD=", round(SD_gmv_lo,6), "\n", sep="")
## 
## GMV long-only:
## NESN.SW: 72.00% | MSFT: 28.00% | NVDA: 0.00%
## E[R]=-0.000206 SD=0.009307
cat("\nCartera rendidora (SD <=", round(SD_obj,6), "):\n", paste(pr(w_risky), collapse=" | "),
    "\nE[R]=", round(ER_risky,6), " SD=", round(SD_risky,6), "\n", sep="")
## 
## Cartera rendidora (SD <=0.011634):
## NESN.SW: 31.00% | MSFT: 37.00% | NVDA: 32.00%
## E[R]=0.000233 SD=0.011633
ord <- order(SD)
sd_s <- SD[ord]
er_s <- ER[ord]

keep <- logical(length(er_s))
m <- -Inf
for (i in seq_along(er_s)) {
  if (er_s[i] > m) {
    keep[i] <- TRUE
    m <- er_s[i]
  }
}

ef_idx <- ord[keep]        
ef_ord <- order(SD[ef_idx])  
library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
SD_e  <- SD[ef_idx][ef_ord]
ER_e  <- ER[ef_idx][ef_ord]
W_e   <- W[ef_idx, ][ef_ord, ]  

fig <- plot_ly(
  x = SD_e,         
  y = ER_e,         
  z = W_e[,1],     
  type = "scatter3d",
  mode = "lines+markers",
  line = list(color = "red", width = 4),
  marker = list(size = 3, color = ER_e, colorscale = "Rainbow")
) %>%
  layout(
    scene = list(
      xaxis = list(title = "Riesgo"),
      yaxis = list(title = "Rendimiento"),
      zaxis = list(title = "Peso")
    ),
    title = "Frontera eficiente en 3D"
  )

fig