options(repos = c(CRAN = “https://cloud.r-project.org”)) pkgs <- c(“tidyverse”,“tidyquant”,“quantmod”,“PerformanceAnalytics”, “quadprog”,“zoo”,“xts”,“TTR”,“scales”,“knitr”,“kableExtra”) to_install <- setdiff(pkgs, rownames(installed.packages())) if (length(to_install)) install.packages(to_install, dependencies = TRUE)
library(tidyverse) library(tidyquant) library(quantmod) library(PerformanceAnalytics) library(quadprog) library(zoo) library(scales) library(knitr) library(kableExtra)
tickers <- c(“AAPL”,“MSFT”,“NVDA”) # cámbialos si quieres start_date <- as.Date(“2023-10-01”) today <- Sys.Date() inv_usd <- 10000000 # 10 millones (sin guion bajo) set.seed(123)
suppressWarnings( try(getSymbols(“DGS10”, src=“FRED”, from=start_date, auto.assign=TRUE), silent=TRUE) ) rf_annual <- if (exists(“DGS10”)) tail(na.locf(DGS10),1)[[1]]/100 else 0.045
```markdown # Precios ajustados px <- tq_get(tickers, get = “stock.prices”, from = start_date) %>% group_by(symbol) %>% select(symbol, date, adjusted)
rets <- px %>% group_by(symbol) %>% arrange(date) %>% mutate(ret = log(adjusted/lag(adjusted))) %>% filter(!is.na(ret)) %>% ungroup()
divs <- map_df(tickers, ~{ d <- suppressWarnings( tryCatch(tq_get(.x, get=“dividends”, from=start_date), error=function(e) NULL) ) tibble(symbol=.x, div_12m = if (!is.null(d) && nrow(d)>0) sum(d$dividend, na.rm=TRUE) else 0) }) last_prices <- px %>% group_by(symbol) %>% summarise(last_px = last(adjusted), .groups=“drop”) div_yield <- divs %>% left_join(last_prices, by=“symbol”) %>% mutate(div_yield = ifelse(last_px>0, div_12m/last_px, 0)) %>% select(symbol, div_yield)
div_yield %>% mutate(div_yield = percent(div_yield)) %>% kable(“html”, caption=“Rendimiento por dividendo (aprox.)”) %>% kable_styling(full_width=FALSE)
R_wide <- rets %>% select(date, symbol, ret) %>% pivot_wider(names_from=symbol, values_from=ret) %>% drop_na()
mu_vec <- colMeans(R_wide[, tickers]) # medias diarias de cada acción Sigma <- cov(R_wide[, tickers]) # covarianza diaria ann_factor <- 252 mu_ann <- mu_vec * ann_factor vol_ann <- sqrt(diag(Sigma) * ann_factor)
tibble(Accion = tickers, mu_anual = mu_ann, sigma_anual = vol_ann) %>% mutate(mu_anual = percent(mu_anual), sigma_anual = percent(sigma_anual)) %>% kable(“html”, caption = “Medias y volatilidades anualizadas (históricas)”) %>% kable_styling(full_width = FALSE)
library(quadprog) library(ggplot2)
n <- length(mu_vec)
Dmat <- 2 * as.matrix(Sigma) + diag(1e-6, n) dvec <- rep(0, n)
Amat <- cbind(rep(1, n), diag(n)) # [n x (1+n)] bvec <- c(1, rep(0, n)) meq <- 1
target <- seq(min(mu_vec)0.9, max(mu_vec)1.1, length.out = 20)
rows <- lapply(target, function(m){ Amat2 <- cbind(Amat, mu_vec) # agrega t(mu) w >= m bvec2 <- c(bvec, m) sol <- tryCatch(solve.QP(Dmat, dvec, Amat2, bvec2, meq = meq), error = function(e) NULL) if (is.null(sol)) return(NULL) data.frame( mu_p = m, sd = sqrt(sol\(value), w1 = sol\)solution[1], w2 = sol\(solution[2], w3 = sol\)solution[3] ) })
rows <- rows[!sapply(rows, is.null)]
if (length(rows) == 0) { sol_mv <- solve.QP(Dmat, dvec, Amat, bvec, meq = meq) fe <- data.frame( mu_p = sum(mu_vec * sol_mv\(solution), sd = sqrt(sol_mv\)value), w1 = sol_mv\(solution[1], w2 = sol_mv\)solution[2], w3 = sol_mv$solution[3] ) } else { fe <- do.call(rbind, rows) }
fe\(mu_ann <- fe\)mu_p * 252 fe\(sd_ann <- fe\)sd * sqrt(252) fe\(sharpe <- (fe\)mu_ann - rf_annual) / fe$sd_ann
opt_idx <- which.max(fe$sharpe) opt_row <- fe[opt_idx, , drop = FALSE] w_opt <- as.numeric(opt_row[, c(“w1”,“w2”,“w3”)]) names(w_opt) <- tickers
print(data.frame(Acción = tickers, Peso = round(w_opt, 3)))
ggplot(fe, aes(x = sd_ann, y = mu_ann)) + geom_line() + geom_point(data = opt_row, aes(x = sd_ann, y = mu_ann), color = “red”, size = 3) + labs(title = “Frontera eficiente (anual)”, x = “Volatilidad anual”, y = “Retorno anual”)
mu_p_daily <- as.numeric(sum(w_opt * mu_vec)) sig_p_daily <- as.numeric(sqrt(t(w_opt) %% Sigma %% w_opt)) mu_p_ann <- mu_p_daily * ann_factor sig_p_ann <- sig_p_daily * sqrt(ann_factor)
S0_p <- sum(w_opt * last_prices$last_px)
quarterly_dates <- seq(from = as.Date(cut(today, “quarter”)), by = “quarter”, length.out = 8)
n_steps <- length(quarterly_dates) dt <- 0.25 # trimestre n_sims <- 5000
sim_paths <- matrix(NA, nrow = n_steps+1, ncol = n_sims) sim_paths[1,] <- S0_p for (j in 1:n_sims){ for (t in 2:(n_steps+1)){ z <- rnorm(1) sim_paths[t,j] <- sim_paths[t-1,j] * exp((mu_p_ann - 0.5sig_p_ann^2)dt + sig_p_annsqrt(dt)z) } }
sim_df <- as_tibble(sim_paths) %>% mutate(step = 0:n_steps, date = c(today, quarterly_dates)) %>% pivot_longer(-c(step,date), names_to=“path”, values_to=“S”)
sim_summary <- sim_df %>% group_by(step, date) %>% summarise(p05 = quantile(S,0.05), p50 = median(S), p95 = quantile(S,0.95), .groups=“drop”)
ggplot(sim_summary, aes(x=date)) + geom_ribbon(aes(ymin=p05, ymax=p95), alpha=0.2) + geom_line(aes(y=p50)) + labs(title=“Simulación GBM del portafolio (trimestral, 2 años)”, x=““, y=”Nivel simulado”)
exp_prices <- sim_df %>% group_by(step, date) %>% summarise(E_S = mean(S), SD_S = sd(S), .groups=“drop”) %>% filter(step>0) exp_prices %>% mutate(E_S = dollar(E_S), SD_S = dollar(SD_S)) %>% kable(“html”, caption=“Precios esperados trimestrales (cartera)”) %>% kable_styling(full_width=FALSE) # ========================= # SEGMENTO 6: Sharpe y VaR (1% y 5%) por trimestre (vía simulación) # ========================= sim_returns <- sim_df %>% group_by(path) %>% arrange(date) %>% mutate(ret = log(S/lag(S))) %>% filter(!is.na(ret))
var_by_q <- sim_returns %>% group_by(date) %>% summarise(VaR_1 = quantile(ret, 0.01), VaR_5 = quantile(ret, 0.05), .groups=“drop”)
var_by_q %>% mutate(VaR_1 = percent(VaR_1), VaR_5 = percent(VaR_5)) %>% kable(“html”, caption=“VaR trimestral (1% y 5%)”) %>% kable_styling(full_width=FALSE)
var_by_q %>% pivot_longer(starts_with(“VaR”), names_to=“nivel”, values_to=“valor”) %>% ggplot(aes(x=date, y=valor, color=nivel)) + geom_line(linewidth = 1) + # <- antes: size = 1 geom_point() + scale_y_continuous(labels=percent) + labs(title=“VaR trimestral del portafolio (simulado)”, x=““, y=”Retorno”)
pnorm_ <- pnorm; dnorm_ <- dnorm bsm_price <- function(S, K, r, q=0, sigma, T, type=c(“call”,“put”)){ type <- match.arg(type) if (T <= 0 || sigma <= 0) return(max((type==“call”)(S - K), (type==“put”)(K - S), 0)) d1 <- (log(S/K) + (r - q + 0.5sigma^2)T) / (sigmasqrt(T)) d2 <- d1 - sigmasqrt(T) if (type==“call”) Sexp(-qT)pnorm_(d1) - Kexp(-rT)pnorm_(d2) else Kexp(-rT)pnorm_(-d2) - Sexp(-qT)pnorm_(-d1) }
bsm_delta <- function(S, K, r, q=0, sigma, T, type=c(“call”,“put”)){ type <- match.arg(type) d1 <- (log(S/K) + (r - q + 0.5sigma^2)T) / (sigmasqrt(T)) if (type==“call”) exp(-qT)pnorm_(d1) else -exp(-qT)*pnorm_(-d1) }
iv_bsm <- function(S, K, r, q=0, T, price_mid, type=c(“call”,“put”), lo=0.01, hi=2, tol=1e-6, maxit=100){ f <- function(sig) bsm_price(S,K,r,q,sig,T,type) - price_mid if (f(lo)*f(hi) > 0) return(NA_real_) uniroot(f, c(lo,hi), tol=tol, maxiter=maxit)$root }
binom_american <- function(S, K, r, q=0, sigma, T, steps=200, type=c(“call”,“put”)){ type <- match.arg(type) dt <- T/steps u <- exp(sigmasqrt(dt)); d <- 1/u p <- (exp((r - q)dt) - d)/(u - d); disc <- exp(-rdt) # valores al vencimiento S_T <- S u^(0:steps) * d^(steps:0) if (type==“call”) V <- pmax(S_T - K, 0) else V <- pmax(K - S_T, 0) # hacia atrás con ejercicio temprano for (i in steps:1){ S_i <- S * u^(0:(i-1)) * d^((i-1):0) V <- disc(pV[2:(i+1)] + (1-p)V[1:i]) if (type==“call”) V <- pmax(V, S_i - K) else V <- pmax(V, K - S_i) } V[1] } # ========================= # SEGMENTO 8 (PLAN B): Valuación con ATM & 90 días usando σ histórica # ========================= # sigma anual histórica por activo (ya tienes Sigma en DIARIO) sigma_hist <- sqrt(diag(Sigma) 252) names(sigma_hist) <- tickers
q_yield <- if (exists(“div_yield”)) { setNames(pmax(div_yield\(div_yield, 0), div_yield\)symbol) } else { setNames(rep(0, length(tickers)), tickers) }
S0_vec <- setNames(last_prices\(last_px, last_prices\)symbol) T_yr <- 90/365
opt_tab <- tibble::tibble( sym = tickers, s0 = as.numeric(S0_vec[tickers]), k = as.numeric(S0_vec[tickers]), # ATM expiry= Sys.Date() + 90, t_yr = T_yr, sigma = as.numeric(sigma_hist[tickers]), q = as.numeric(q_yield[tickers]) )
opt_tab <- opt_tab %>% rowwise() %>% mutate( c_eur = bsm_price(s0, k, r = rf_annual, q = q, sigma = sigma, T = t_yr, type = “call”), p_eur = bsm_price(s0, k, r = rf_annual, q = q, sigma = sigma, T = t_yr, type = “put”), c_amer = binom_american(s0, k, r = rf_annual, q = q, sigma = sigma, T = t_yr, steps = 200, type = “call”), p_amer = binom_american(s0, k, r = rf_annual, q = q, sigma = sigma, T = t_yr, steps = 200, type = “put”) ) %>% ungroup()
opt_tab %>% transmute( Ticker = sym, Spot = round(s0,2), Strike =
round(k,2), Expira = expiry, σ hist (anual) =
round(sigma,3), IV (plan B) = round(sigma,3),
Call Eur = round(c_eur,2), Put Eur =
round(p_eur,2), Call Amer= round(c_amer,2),
Put Amer= round(p_amer,2) ) %>% kable(“html”, caption =
“Valuación de opciones (plan B: σ histórica, ATM, 90 días)”) %>%
kable_styling(full_width = FALSE)
L_total <- inv_usd * 0.85 # dinero apalancado alloc_L <- w_opt / sum(w_opt) * L_total # repartir según pesos óptimos
shares_i <- (inv_usd * w_opt) / last_prices\(last_px[match(names(w_opt), last_prices\)symbol)] names(shares_i) <- names(w_opt)
cover_frac <- 0.85 contracts_plan <- opt_tab %>% mutate( Peso = w_opt[match(sym, names(w_opt))], Shares = shares_i[match(sym, names(shares_i))], Ctr_target = ceiling(cover_frac * Shares / 100), # Costo de puts europeos (como referencia de prima; también puedes usar Americana) Premium_put = P_eur, Cash_need = Ctr_target * 100 * Premium_put, Budget = alloc_L[match(sym, names(alloc_L))], Ctr_final = pmax(0, floor(Budget / (100 * Premium_put))), # ajusta a presupuesto Cobertura_shares = Ctr_final * 100, Cobertura_pct = pmin(1, Cobertura_shares / Shares) )
cover_frac <- 0.85
L_total <- inv_usd * 0.85 alloc_L <- w_opt / sum(w_opt) * L_total
contracts_plan <- opt_tab %>% dplyr::mutate(+ Peso = w_opt[match(sym, names(w_opt))], Shares = if (exists(“shares_i”)) shares_i[match(sym, names(shares_i))] else (inv_usd * Peso) / s0, ctr_target = ceiling(cover_frac * Shares / 100), # 100 sh/contrato premium_put = p_eur, # <- minúsculas budget = alloc_L[match(sym, names(alloc_L))], ctr_final = pmax(0, floor(budget / (100 * premium_put))), covered_shares = ctr_final * 100, covered_pct = pmin(1, covered_shares / Shares) )
contracts_plan %>% dplyr::transmute( ticker = sym, expira = expiry, spot = round(s0, 2), strike = round(k, 2), put_price = round(premium_put, 2), acciones = round(Shares), contratos_target = ctr_target, contratos_comprados = ctr_final, cobertura_pct = scales::percent(covered_pct), presupuesto = scales::dollar(budget), costo = scales::dollar(ctr_final * 100 * premium_put) ) %>% kable(“html”, caption = “Plan de cobertura con PUTs (85% apalancado, rolling trimestral)”) %>% kable_styling(full_width = FALSE)
cat(” Estrategia de cobertura: se adquieren opciones PUT at-the-money con vencimiento ~90 días, una por cada activo del paquete, y se renuevan trimestralmente (rolling) hasta 2 años. Esto limita pérdidas del portafolio subyacente, especialmente en choques de mercado. La cobertura se financia con un apalancamiento del 85% a la tasa del bono 10Y (rf_annual). La asignación del dinero apalancado se hace proporcional a los pesos del portafolio óptimo. Para referencia, también se valúan las opciones americanas (binomial) y europeas (BSM). “)
cat(” ## Conclusiones
El portafolio de tres acciones fue optimizado con base en el modelo media-varianza, obteniendo un portafolio de máxima eficiencia bajo el criterio de Sharpe. Se simularon precios mediante un proceso de movimiento browniano geométrico (MGB), y se evaluaron las métricas de riesgo (VaR al 1% y 5%) junto al comportamiento trimestral esperado.
Para la cobertura, se valoraron opciones PUT tanto europeas como americanas utilizando volatilidades históricas y vencimiento a 90 días (rolling trimestral), cubriendo aproximadamente el 85% de la posición total mediante apalancamiento con la tasa del bono del tesoro (rf_annual). Esto permite mitigar las pérdidas potenciales del portafolio ante escenarios bajistas. “)