# Caso de Política Pública

La **Dirección de Programas Sociales** debe presentar un presupuesto al Congreso para un programa de transferencias monetarias condicionadas. El problema: **múltiples fuentes de incertidumbre**:

-  **Costo unitario** por beneficiario (varía por logística, inflación)
-  **Número de beneficiarios** (demanda fluctuante)
-  **Costos logísticos** (transporte, operación regional)
-  **Tasa de deserción** (beneficiarios que abandonan el programa)

**Tu tarea:** Usar simulación Monte Carlo para:
- Estimar el presupuesto necesario con su rango de incertidumbre
- Identificar el percentil 90 (para no quedarse corto)
- Comunicar el riesgo presupuestario al Director

**Regla de oro:** Planificar con el promedio es planificar para el fracaso. Necesitamos conocer la distribución completa.

---

# Parte 1: Concepto de Simulación Monte Carlo

**¿Qué es Monte Carlo?**
- Técnica que usa muestreo aleatorio repetido
- "Tira los dados" miles de veces
- Genera distribución de resultados posibles
- Usado en finanzas, seguros, ingeniería, salud pública

**¿Por qué lo necesitamos?**
En política pública, raramente sabemos un valor exacto. Sabemos:
- "El costo unitario está entre $10-15 USD, pero más probable cerca de $12"
- "Esperamos 5,000 beneficiarios ± 400"
- "La logística cuesta entre $30k-50k USD"

Monte Carlo convierte esta incertidumbre en números útiles para decidir.

---

# Parte 2: Componentes del Modelo

```{webr-r}
#| warning: false

# Configuración de simulación
set.seed(2026)
B <- 10000  # Número de simulaciones (10 mil escenarios)

cat(" CONFIGURACIÓN DE LA SIMULACIÓN\n")
cat("═══════════════════════════════════════════════════════════\n")
cat("Corriendo", B, "escenarios posibles...\n\n")

# COMPONENTE 1: Costo unitario por beneficiario (distribución log-normal)
# Log-normal porque: no puede ser negativo, tiene cola derecha (inflación extrema)
# Media esperada: $12 USD, pero con variabilidad moderada

# Parámetros log-normal para media≈12, SD≈3
mu_log <- log(12)
sigma_log <- 0.30  # Variabilidad moderada

costo_unitario <- rlnorm(B, meanlog = mu_log, sdlog = sigma_log)

cat(" COSTO UNITARIO (por beneficiario)\n")
cat("   Distribución: Log-normal\n")
cat("   Justificación: Costos siempre positivos, cola derecha por inflación\n")
cat("   Media:      $", round(mean(costo_unitario), 2), " USD\n")
cat("   Mediana:    $", round(median(costo_unitario), 2), " USD\n")
cat("   DE:         $", round(sd(costo_unitario), 2), " USD\n")
cat("   Rango:      $", round(min(costo_unitario), 2), " - $", 
    round(max(costo_unitario), 2), " USD\n\n")

# COMPONENTE 2: Número de beneficiarios (distribución normal)
# Normal porque: suma de muchos factores independientes
# Media esperada: 5,000, SD: 600 (mayor incertidumbre)

n_beneficiarios <- round(rnorm(B, mean = 5000, sd = 600))
# Asegurar que no haya valores negativos (truncación realista)
n_beneficiarios <- pmax(n_beneficiarios, 3000)

cat(" NÚMERO DE BENEFICIARIOS\n")
cat("   Distribución: Normal (truncada en 3,000)\n")
cat("   Justificación: Suma de inscripciones de múltiples municipios\n")
cat("   Media:      ", round(mean(n_beneficiarios), 0), " personas\n")
cat("   Mediana:    ", round(median(n_beneficiarios), 0), " personas\n")
cat("   DE:         ", round(sd(n_beneficiarios), 0), " personas\n")
cat("   Rango:      ", min(n_beneficiarios), " - ", 
    max(n_beneficiarios), " personas\n\n")

# COMPONENTE 3: Costos logísticos con CORRELACIÓN FUERTE
# Gamma base PEQUEÑA, correlación GRANDE para enfatizar dependencia
# Justificación: Más beneficiarios → proporcionalmente más rutas y personal

# Primero generar gamma base con MENOS variabilidad independiente
costos_logistica_base <- rgamma(B, shape = 8, rate = 1/5000)  # Más concentrada

# Añadir componente correlacionado MUY FUERTE (driver principal)
n_beneficiarios_std <- (n_beneficiarios - mean(n_beneficiarios)) / sd(n_beneficiarios)
costos_logistica <- costos_logistica_base + n_beneficiarios_std * 10000

# Asegurar valores positivos
costos_logistica <- pmax(costos_logistica, 15000)

cat(" COSTOS LOGÍSTICOS (transporte, operación)\n")
cat("   Distribución: Gamma (con correlación FUERTE)\n")
cat("   Justificación: Costos logísticos dependen fuertemente de # beneficiarios\n")
cat("   Media:      $", round(mean(costos_logistica), 0), " USD\n")
cat("   Mediana:    $", round(median(costos_logistica), 0), " USD\n")
cat("   DE:         $", round(sd(costos_logistica), 0), " USD\n")
cat("   Rango:      $", round(min(costos_logistica), 0), " - $", 
    round(max(costos_logistica), 0), " USD\n\n")

# CORRELACIÓN ENTRE VARIABLES
cat(" CORRELACIONES ENTRE VARIABLES\n")
cat("   Corr(Beneficiarios ↔ Logística): ", 
    round(cor(n_beneficiarios, costos_logistica), 2), "\n")
cat("   Justificación: Más beneficiarios requieren proporcionalmente\n")
cat("                  más rutas, personal y coordinación logística\n\n")

cat("═══════════════════════════════════════════════════════════\n")
```
---

## INTERPRETACIÓN DE PARÁMETROS: COSTOS LOGÍSTICOS

Usamos distribución Gamma(shape = 8, rate = 1/5000) para el componente base.

**¿Cómo se traduce esto a dólares?**

Para una distribución Gamma:

- Media = shape / rate = 8 / (1/5000) = 40,000 USD
- Desviación estándar = √shape / rate = √8 / (1/5000) ≈ 14,142 USD

### COMPONENTE BASE: $40,000 ± $14,000

Representa costos logísticos "típicos" con variación aleatoria local (clima, geografía, eficiencia administrativa)

### COMPONENTE CORRELACIONADO: +$10,000 × (beneficiarios estandarizados)

Más beneficiarios → proporcionalmente más rutas, vehículos, personal

**Ejemplos:**

- Municipio pequeño (800 beneficiarios, -1 SD): -$10,000 ajuste
- Municipio promedio (1000 beneficiarios): $0 ajuste  
- Municipio grande (1200 beneficiarios, +1 SD): +$10,000 ajuste

### PISO MÍNIMO: $15,000

Costos fijos inevitables (coordinador, vehículo base)

### RESULTADO: 

- Costos logísticos entre $25,000 y $70,000 según tamaño del municipio.
- Es solo álgebra para "calibrar" la distribución a dólares reales.


---

**Preguntas de reflexión (discutir en clase):**

1. ¿Por qué usamos distribuciones diferentes para cada componente?
2. ¿Qué pasaría si usáramos Normal para todo?
3. ¿Por qué el costo unitario es Lognormal y no Normal?
   - *Pista:* ¿Puede ser negativo un costo?
4. ¿La truncación en 3,000 beneficiarios es realista? ¿Qué implica?
5. ¿Por qué es importante modelar la correlación entre beneficiarios y logística?

---

# Parte 2.5: La Magia Detrás del Telón

## ¿Cómo genera R números aleatorios con diferentes distribuciones?

```{webr-r}
# Demostración del método de transformación inversa
# PRINCIPIO: Todo empieza con Uniforme(0,1)

set.seed(2026)
n_demo <- 1000

# PASO 1: Generar números aleatorios uniformes entre 0 y 1
u <- runif(n_demo)

# PASO 2: Transformar a Normal(0,1) usando función cuantil inversa
x_normal <- qnorm(u)

# PASO 3: Transformar a Lognormal (como nuestros costos)
x_lognorm <- qlnorm(u, meanlog = log(12), sdlog = 0.30)

# PASO 4: Transformar a Gamma (como nuestra logística)
x_gamma <- qgamma(u, shape = 8, rate = 1/5000)

# Visualización del proceso
par(mfrow = c(2, 3), mar = c(4, 4, 3, 1))

# Fila 1: Las distribuciones
hist(u, breaks = 30, col = "lightgray", border = "white",
     main = "1. Uniforme(0,1)\n[Inicio]",
     xlab = "Valor", ylab = "Frecuencia", xlim = c(0, 1))

hist(x_normal, breaks = 30, col = "steelblue", border = "white",
     main = "2. Normal(0,1)\n[Transformada]",
     xlab = "Valor", ylab = "Frecuencia")

hist(x_lognorm, breaks = 30, col = "coral", border = "white",
     main = "3. Lognormal\n[Nuestros costos]",
     xlab = "Costo ($)", ylab = "Frecuencia")

# Fila 2: Las funciones de transformación
plot(u[1:200], x_normal[1:200], 
     pch = 16, col = rgb(0,0,1,0.3), cex = 0.8,
     main = "Transformación: U→Normal",
     xlab = "Uniforme(0,1)", ylab = "Normal")
curve(qnorm(x), from = 0.01, to = 0.99, add = TRUE, col = "darkblue", lwd = 2)

plot(u[1:200], x_lognorm[1:200], 
     pch = 16, col = rgb(1,0.5,0,0.3), cex = 0.8,
     main = "Transformación: U→Lognormal",
     xlab = "Uniforme(0,1)", ylab = "Lognormal")
curve(qlnorm(x, meanlog = log(12), sdlog = 0.30), 
      from = 0.01, to = 0.99, add = TRUE, col = "darkorange", lwd = 2)

plot(u[1:200], x_gamma[1:200], 
     pch = 16, col = rgb(0,0.5,0,0.3), cex = 0.8,
     main = "Transformación: U→Gamma",
     xlab = "Uniforme(0,1)", ylab = "Gamma")
curve(qgamma(x, shape = 8, rate = 1/5000), 
      from = 0.01, to = 0.99, add = TRUE, col = "darkgreen", lwd = 2)

par(mfrow = c(1, 1))

cat("\n PRINCIPIO FUNDAMENTAL DE SIMULACIÓN\n")
cat("═══════════════════════════════════════════════════════════\n\n")
cat("1. Todo generador de números aleatorios empieza con U(0,1)\n")
cat("2. Las funciones qnorm(), qlnorm(), qgamma() aplican la 'transformación inversa'\n")
cat("3. Esto convierte uniformes en CUALQUIER distribución\n\n")

cat("Ejemplo práctico:\n")
cat("  u = 0.75  (un número aleatorio uniforme)\n")
cat("  qnorm(0.75) = ", round(qnorm(0.75), 2), "  (percentil 75 de Normal)\n")
cat("  qlnorm(0.75, log(12), 0.30) = $", 
    round(qlnorm(0.75, log(12), 0.30), 2), "  (percentil 75 de nuestros costos)\n\n")

cat("Por eso decimos que las distribuciones son 'percentiles invertidos'!\n")
cat("═══════════════════════════════════════════════════════════\n")
```

**¿Por qué importa esto?**
- Entender que toda aleatoriedad es transformación matemática
- Saber que podemos crear *cualquier* distribución a partir de U(0,1)
- Reconocer que `rnorm()`, `rlnorm()` hacen esto automáticamente

---

# Parte 2.6: Test de Convergencia

## ¿10,000 simulaciones son suficientes?

```{webr-r}
# Probar diferentes tamaños de muestra
tamanios <- c(100, 500, 1000, 5000, 10000, 50000)
p90_por_tamano <- numeric(length(tamanios))

set.seed(2026)

cat("\n TEST DE CONVERGENCIA\n")
cat("═══════════════════════════════════════════════════════════\n")
cat("Probando estabilidad del P90 con diferentes números de simulaciones...\n\n")

for (i in seq_along(tamanios)) {
  B_temp <- tamanios[i]
  cu_temp <- rlnorm(B_temp, meanlog = mu_log, sdlog = sigma_log)
  nb_temp <- pmax(round(rnorm(B_temp, mean = 5000, sd = 600)), 3000)
  cl_base <- rgamma(B_temp, shape = 8, rate = 1/5000)
  nb_std <- (nb_temp - mean(nb_temp)) / sd(nb_temp)
  cl_temp <- pmax(cl_base + nb_std * 10000, 15000)
  pt_temp <- (cu_temp * nb_temp) + cl_temp
  p90_por_tamano[i] <- quantile(pt_temp, 0.90)
}

# Gráfico de convergencia
plot(tamanios, p90_por_tamano / 1000,
     type = "b", pch = 16, col = "darkblue", lwd = 2,
     main = "Convergencia del P90: ¿Cuántas simulaciones necesito?",
     xlab = "Número de simulaciones",
     ylab = "P90 estimado (miles USD)",
     log = "x",
     ylim = range(p90_por_tamano / 1000) + c(-2, 2))

# Línea de referencia (10,000 sims)
abline(h = p90_por_tamano[5] / 1000, col = "red", lty = 2, lwd = 2)
text(200, p90_por_tamano[5]/1000 + 1, 
     "Valor con 10,000 sims", pos = 4, col = "red", cex = 0.9)

# Zona de estabilidad
rect(5000, min(p90_por_tamano/1000) - 2, 50000, max(p90_por_tamano/1000) + 2,
     col = rgb(0, 1, 0, 0.1), border = NA)
text(15000, max(p90_por_tamano/1000) + 1.5, 
     "Zona de estabilidad", col = "darkgreen", cex = 0.9)

cat("P90 estimado con diferentes números de simulaciones:\n")
for (i in seq_along(tamanios)) {
  diferencia <- abs(p90_por_tamano[i] - p90_por_tamano[5]) / p90_por_tamano[5] * 100
  cat(sprintf("  %6d sims: $%s  (diferencia vs 10k: %0.2f%%)\n", 
              tamanios[i], 
              format(round(p90_por_tamano[i], 0), big.mark = ","),
              diferencia))
}

cat("\n CONCLUSIÓN: Con 10,000+ simulaciones, el P90 se estabiliza.\n")
cat("   Variación < 1% entre 10,000 y 50,000 simulaciones.\n\n")
cat("═══════════════════════════════════════════════════════════\n")
```

**Lección importante:** Más simulaciones = mayor precisión, pero con rendimientos decrecientes.

---

# Parte 3: Cálculo del Presupuesto Total

```{webr-r}
# MODELO: Presupuesto Total = (Costo Unitario × Beneficiarios) + Logística
presupuesto_total <- (costo_unitario * n_beneficiarios) + costos_logistica

cat("\n💵 PRESUPUESTO TOTAL SIMULADO\n")
cat("═══════════════════════════════════════════════════════════\n\n")

# Estadísticas descriptivas
cat("Estadísticas centrales:\n")
cat("  Media:      $", format(round(mean(presupuesto_total), 0), 
                              big.mark = ","), " USD\n")
cat("  Mediana:    $", format(round(median(presupuesto_total), 0), 
                              big.mark = ","), " USD\n")
cat("  DE:         $", format(round(sd(presupuesto_total), 0), 
                              big.mark = ","), " USD\n\n")

# Percentiles críticos
percentiles <- quantile(presupuesto_total, 
                        probs = c(0.05, 0.10, 0.25, 0.50, 0.75, 0.90, 0.95))

cat("Percentiles (rango de escenarios posibles):\n")
cat("  P5  (optimista):     $", format(round(percentiles[1], 0), 
                                        big.mark = ","), "\n")
cat("  P10:                 $", format(round(percentiles[2], 0), 
                                        big.mark = ","), "\n")
cat("  P25 (cuartil inf):   $", format(round(percentiles[3], 0), 
                                        big.mark = ","), "\n")
cat("  P50 (mediana):       $", format(round(percentiles[4], 0), 
                                        big.mark = ","), "\n")
cat("  P75 (cuartil sup):   $", format(round(percentiles[5], 0), 
                                        big.mark = ","), "\n")
cat("  P90 (conservador):   $", format(round(percentiles[6], 0), 
                                        big.mark = ","), "\n")
cat("  P95 (muy conservador): $", format(round(percentiles[7], 0), 
                                          big.mark = ","), "\n\n")

# Rango intercuartil
iqr <- IQR(presupuesto_total)
cat("Rango Intercuartil (50% central): $", 
    format(round(iqr, 0), big.mark = ","), "\n\n")

# Coeficiente de variación
cv <- 100 * sd(presupuesto_total) / mean(presupuesto_total)
cat("Coeficiente de Variación: ", round(cv, 1), "%\n")
if (cv > 15) {
  cat("  Alta incertidumbre - considere planeación por escenarios\n")
} else if (cv > 10) {
  cat("⚡ Incertidumbre moderada - use percentiles conservadores\n")
} else {
  cat(" Incertidumbre baja - modelo relativamente predecible\n")
}

cat("═══════════════════════════════════════════════════════════\n")
```

---

# Parte 4: Visualización de la Distribución

```{webr-r}
# Histograma de presupuestos simulados
hist(presupuesto_total / 1000,  # Convertir a miles para mejor lectura
     breaks = 60,
     col = "steelblue",
     border = "white",
     main = "Distribución de Presupuesto Total Simulado (10,000 escenarios)",
     xlab = "Presupuesto Total (miles USD)",
     ylab = "Frecuencia de escenarios",
     xlim = c(40, 120))

# Líneas de referencia para percentiles clave
abline(v = percentiles[4] / 1000, col = "darkgreen", lwd = 3, lty = 1)  # P50
abline(v = percentiles[3] / 1000, col = "orange", lwd = 2, lty = 2)     # P25
abline(v = percentiles[5] / 1000, col = "orange", lwd = 2, lty = 2)     # P75
abline(v = percentiles[6] / 1000, col = "red", lwd = 3, lty = 3)        # P90

# Leyenda
legend("topright",
       legend = c("P50 (Mediana)", "P25/P75 (IQR)", "P90 (Conservador)"),
       col = c("darkgreen", "orange", "red"),
       lty = c(1, 2, 3),
       lwd = c(3, 2, 3),
       bty = "n",
       cex = 0.9)

# Anotación
text(x = percentiles[6] / 1000 + 5, 
     y = max(hist(presupuesto_total/1000, plot=FALSE, breaks=60)$counts) * 0.7,
     paste0("10% de escenarios\nexceden $", 
            format(round(percentiles[6]/1000, 0), big.mark=","), "k"),
     pos = 4, col = "red", cex = 0.9)
```

**Interpreta el histograma:**
- ¿La distribución es simétrica o asimétrica?
- ¿Hay cola larga a la derecha? (escenarios de alto costo)
- ¿Qué tan separados están P25 y P75?
- ¿Qué implica la forma de la distribución para tu planificación?

---

# Parte 5: Análisis de Riesgo Presupuestario

```{webr-r}
# Definir diferentes escenarios de presupuesto más realistas
escenarios <- c(60000, 70000, 80000, 90000)
nombres_esc <- c("Austero", "Moderado", "Prudente", "Conservador")

cat("\n ANÁLISIS DE RIESGO: ¿Qué presupuesto solicitar?\n")
cat("═══════════════════════════════════════════════════════════\n\n")

for (i in 1:length(escenarios)) {
  prob_exceder <- mean(presupuesto_total > escenarios[i]) * 100
  prob_suficiente <- 100 - prob_exceder
  monto_faltante <- mean(pmax(presupuesto_total - escenarios[i], 0))
  
  cat("Escenario", nombres_esc[i], "($", 
      format(escenarios[i], big.mark = ","), "):\n")
  cat("   Probabilidad de ser SUFICIENTE:  ", sprintf("%4.1f%%", prob_suficiente), "\n")
  cat("  ❌ Probabilidad de QUEDARSE CORTO:  ", sprintf("%4.1f%%", prob_exceder), "\n")
  
  # Evaluar nivel de riesgo
  if (prob_exceder > 20) {
    cat("    RIESGO ALTO – Más del 20% de escenarios rebasan\n\n")
  } else if (prob_exceder > 10) {
    cat("  ⚡ RIESGO MODERADO – Entre 10-20% de escenarios rebasan\n\n")
  } else {
    cat("  ✓  RIESGO BAJO – Menos del 10% de escenarios rebasan\n\n")
  }
}

# Gráfico de probabilidades
probabilidades <- sapply(escenarios, function(e) mean(presupuesto_total > e) * 100)

barplot(probabilidades,
        names.arg = nombres_esc,
        col = c("darkred", "orange", "gold", "lightgreen"),
        main = "Riesgo de Exceder el Presupuesto",
        ylab = "Probabilidad de exceder (%)",
        ylim = c(0, 100),
        border = "white")
abline(h = c(10, 20, 50), lty = 2, col = "gray", lwd = 1.5)
text(0.5, 10, "10%", pos = 3, col = "gray30", cex = 0.8)
text(0.5, 20, "20%", pos = 3, col = "gray30", cex = 0.8)
text(0.5, 50, "50%", pos = 3, col = "gray30", cex = 0.8)

text(1:4, probabilidades + 5, 
     paste0(round(probabilidades, 1), "%"), 
     cex = 1.0, font = 2)

cat("═══════════════════════════════════════════════════════════\n\n")

cat("💬 Pregunta de política: ¿Qué nivel de riesgo es aceptable?\n")
cat("   • ¿Preferirías quedarte corto en 5% o 20% de los escenarios?\n")
cat("   • ¿Qué costo político tiene pedir presupuesto suplementario?\n")
cat("   • ¿Qué costo tiene sobre-estimar y que sobre presupuesto?\n\n")
```

**Pregunta clave:** ¿Cuánto riesgo está dispuesta a asumir tu institución?

---

# Parte 6: Análisis de Sensibilidad (Tornado Plot)

## ¿Qué variable genera más incertidumbre?

```{webr-r}
# Calcular contribución a la varianza de cada componente
# Método: Fijar cada variable en su media y ver cuánto disminuye la varianza

var_total <- var(presupuesto_total)

# Presupuesto con cada componente fijo en su media
presup_sin_costo <- (mean(costo_unitario) * n_beneficiarios) + costos_logistica
presup_sin_benef <- (costo_unitario * mean(n_beneficiarios)) + costos_logistica
presup_sin_logist <- (costo_unitario * n_beneficiarios) + mean(costos_logistica)

# Reducción de varianza al fijar cada componente
reduccion_costo <- (var_total - var(presup_sin_costo)) / var_total * 100
reduccion_benef <- (var_total - var(presup_sin_benef)) / var_total * 100
reduccion_logist <- (var_total - var(presup_sin_logist)) / var_total * 100

contribuciones <- c(reduccion_costo, reduccion_benef, reduccion_logist)
nombres <- c("Costo\nunitario", "Número\nbeneficiarios", "Costos\nlogística")

cat("\n ANÁLISIS DE SENSIBILIDAD: ¿Qué input genera más incertidumbre?\n")
cat("═══════════════════════════════════════════════════════════\n\n")

# Tornado plot
par(mar = c(5, 10, 4, 2))
barplot(sort(contribuciones, decreasing = TRUE),
        horiz = TRUE,
        names.arg = nombres[order(contribuciones, decreasing = TRUE)],
        col = c("darkred", "orange", "gold"),
        main = "Tornado Plot: Contribución a la Incertidumbre Total",
        xlab = "Reducción de varianza al eliminar incertidumbre (%)",
        las = 1,
        xlim = c(0, 100))
abline(v = 50, lty = 2, col = "gray", lwd = 2)
text(50, 1.8, "50%", pos = 3, col = "gray30", cex = 0.9)

# Agregar valores
text(sort(contribuciones, decreasing = TRUE) + 3, 
     seq(0.7, by = 1.2, length.out = 3),
     paste0(round(sort(contribuciones, decreasing = TRUE), 1), "%"),
     pos = 4, cex = 0.9, font = 2)

par(mar = c(5, 4, 4, 2))

cat("Contribución de cada componente a la incertidumbre total:\n\n")
orden <- order(contribuciones, decreasing = TRUE)
for (i in 1:3) {
  cat(sprintf("  %d. %s: %0.1f%%\n", 
              i, 
              gsub("\n", " ", nombres[orden[i]]), 
              contribuciones[orden[i]]))
}

cat("\n IMPLICACIÓN ESTRATÉGICA:\n")
cat("   Para reducir incertidumbre, enfócate en: ", 
    gsub("\n", " ", nombres[which.max(contribuciones)]), "\n")
cat("   Esta variable explica el ", round(max(contribuciones), 1), 
    "% de la variación total.\n\n")

# Recomendaciones específicas
if (which.max(contribuciones) == 1) {
  cat("   ESTRATEGIA RECOMENDADA:\n")
  cat("   • Negociar contratos a precio fijo con proveedores\n")
  cat("   • Compras anticipadas para mitigar inflación\n")
  cat("   • Hedging de tipo de cambio si hay insumos importados\n\n")
} else if (which.max(contribuciones) == 2) {
  cat("   ESTRATEGIA RECOMENDADA:\n")
  cat("   • Pre-registro obligatorio de beneficiarios\n")
  cat("   • Campañas de difusión para estimar demanda real\n")
  cat("   • Análisis de inscripciones históricas por región\n\n")
} else {
  cat("   ESTRATEGIA RECOMENDADA:\n")
  cat("   • Licitaciones de logística con cláusulas de precio máximo\n")
  cat("   • Consolidación de rutas de distribución\n")
  cat("   • Alianzas con distribuidores locales\n\n")
}

cat("═══════════════════════════════════════════════════════════\n")
```

---

# Parte 7: Análisis de Correlación

```{webr-r}
# Scatter plot mostrando la correlación entre beneficiarios y logística
plot(n_beneficiarios, costos_logistica / 1000,
     pch = 16, col = rgb(0, 0, 1, 0.2),
     main = "Relación entre Beneficiarios y Costos Logísticos",
     xlab = "Número de beneficiarios",
     ylab = "Costos logísticos (miles USD)",
     cex = 0.6)

# Agregar línea de tendencia
abline(lm(costos_logistica ~ n_beneficiarios), col = "red", lwd = 2)

# Calcular correlación
corr_value <- cor(n_beneficiarios, costos_logistica)
text(min(n_beneficiarios) + 500, max(costos_logistica)/1000 - 5,
     paste0("Correlación: ", round(corr_value, 2)),
     pos = 4, col = "red", cex = 1.2, font = 2)

cat("\n ANÁLISIS DE CORRELACIÓN\n")
cat("═══════════════════════════════════════════════════════════\n\n")
cat("Correlación (Beneficiarios ↔ Logística): ", round(corr_value, 3), "\n\n")

cat("¿Por qué importa la correlación?\n")
cat("  • Variables correlacionadas amplifican la incertidumbre\n")
cat("  • Ignorar correlaciones puede SUBESTIMAR el riesgo real\n")
cat("  • Modelar dependencias hace la simulación más realista\n\n")

cat("Comparación (simulación sin vs con correlación):\n\n")

# Simular escenario SIN correlación (gamma independiente)
set.seed(2026)
costos_log_indep <- pmax(rgamma(B, shape = 8, rate = 1/5000) + 40000, 15000)
presup_sin_corr <- (costo_unitario * n_beneficiarios) + costos_log_indep

cat("  SIN correlación:\n")
cat("    P90 =  $", format(round(quantile(presup_sin_corr, 0.90), 0), big.mark = ","), "\n")
cat("    IQR =  $", format(round(IQR(presup_sin_corr), 0), big.mark = ","), "\n\n")

cat("  CON correlación (modelo actual):\n")
cat("    P90 =  $", format(round(percentiles[6], 0), big.mark = ","), "\n")
cat("    IQR =  $", format(round(iqr, 0), big.mark = ","), "\n\n")

diferencia_p90 <- (percentiles[6] - quantile(presup_sin_corr, 0.90)) / 
                   quantile(presup_sin_corr, 0.90) * 100

cat("    Diferencia en P90: ", sprintf("%+.1f%%", diferencia_p90), "\n")
if (abs(diferencia_p90) > 5) {
  cat("     ¡Ignorar la correlación cambiaría significativamente tu presupuesto!\n")
  cat("     Subestimaríamos el riesgo en $", 
      format(round(percentiles[6] - quantile(presup_sin_corr, 0.90), 0), big.mark = ","), 
      " USD\n")
}

cat("\n═══════════════════════════════════════════════════════════\n")
```

---

# Parte 8: Simulación de Estrategias de Mitigación

```{webr-r}
# ESTRATEGIA A: Contratos fijos (reduce varianza del costo unitario)
set.seed(2026)
cu_estrategia_a <- rlnorm(B, meanlog = log(12), sdlog = 0.18)  # Menos volatilidad
presup_estrategia_a <- (cu_estrategia_a * n_beneficiarios) + costos_logistica

# ESTRATEGIA B: Pre-registro (reduce varianza de beneficiarios)
nb_estrategia_b <- pmax(round(rnorm(B, mean = 5000, sd = 350)), 3500)  # Menos volatilidad
nb_std_b <- (nb_estrategia_b - mean(nb_estrategia_b)) / sd(nb_estrategia_b)
cl_estrategia_b <- pmax(rgamma(B, shape = 8, rate = 1/5000) + nb_std_b * 10000, 15000)
presup_estrategia_b <- (costo_unitario * nb_estrategia_b) + cl_estrategia_b

# ESTRATEGIA C: Contratos + Pre-registro (combinada)
nb_c <- pmax(round(rnorm(B, mean = 5000, sd = 350)), 3500)
nb_std_c <- (nb_c - mean(nb_c)) / sd(nb_c)
cl_c <- pmax(rgamma(B, shape = 8, rate = 1/5000) + nb_std_c * 10000, 15000)
presup_estrategia_c <- (cu_estrategia_a * nb_c) + cl_c

# BASE (sin estrategias)
presup_base <- presupuesto_total

cat("\n🛡️ SIMULACIÓN DE ESTRATEGIAS DE MITIGACIÓN\n")
cat("═══════════════════════════════════════════════════════════\n\n")

estrategias <- list(
  "Base (sin mitigación)" = presup_base,
  "A: Contratos fijos" = presup_estrategia_a,
  "B: Pre-registro" = presup_estrategia_b,
  "C: Combinada (A+B)" = presup_estrategia_c
)

cat("Comparación de estrategias:\n\n")
cat(sprintf("%-25s | %10s | %10s | %10s | %5s\n", 
            "Estrategia", "P50", "P90", "IQR", "CV%"))
cat(strrep("─", 75), "\n")

for (nombre in names(estrategias)) {
  p50 <- quantile(estrategias[[nombre]], 0.50)
  p90 <- quantile(estrategias[[nombre]], 0.90)
  iqr_val <- IQR(estrategias[[nombre]])
  cv_val <- sd(estrategias[[nombre]]) / mean(estrategias[[nombre]]) * 100
  
  cat(sprintf("%-25s | $%9s | $%9s | $%9s | %4.1f%%\n",
              nombre,
              format(round(p50, 0), big.mark = ","),
              format(round(p90, 0), big.mark = ","),
              format(round(iqr_val, 0), big.mark = ","),
              cv_val))
}

cat("\n")

# Calcular reducciones
reduccion_iqr_a <- (IQR(presup_base) - IQR(presup_estrategia_a)) / IQR(presup_base) * 100
reduccion_iqr_b <- (IQR(presup_base) - IQR(presup_estrategia_b)) / IQR(presup_base) * 100
reduccion_iqr_c <- (IQR(presup_base) - IQR(presup_estrategia_c)) / IQR(presup_base) * 100

cat("Reducción de incertidumbre (IQR) vs baseline:\n")
cat("  Estrategia A: -", round(reduccion_iqr_a, 1), "%\n")
cat("  Estrategia B: -", round(reduccion_iqr_b, 1), "%\n")
cat("  Estrategia C: -", round(reduccion_iqr_c, 1), "%\n\n")

cat(" RECOMENDACIÓN:\n")
mejor_estrategia <- names(estrategias)[which.min(sapply(estrategias, IQR))]
cat("   La ", mejor_estrategia, " minimiza la incertidumbre.\n\n")

cat("═══════════════════════════════════════════════════════════\n")

# Visualización comparativa
par(mar = c(8, 4, 4, 2))
boxplot(presup_base / 1000,
        presup_estrategia_a / 1000,
        presup_estrategia_b / 1000,
        presup_estrategia_c / 1000,
        names = c("Base", "Contratos\nfijos", "Pre-\nregistro", "Combinada"),
        col = c("lightgray", "lightblue", "lightgreen", "gold"),
        main = "Comparación de Estrategias de Mitigación de Riesgo",
        ylab = "Presupuesto Total (miles USD)",
        xlab = "",
        las = 2,
        ylim = c(40, 120))

# Línea de referencia P90 base
abline(h = percentiles[6] / 1000, col = "red", lty = 2, lwd = 2)
text(1.5, percentiles[6]/1000 + 4, 
     "P90 sin mitigación", pos = 4, col = "red", cex = 0.9)

par(mar = c(5, 4, 4, 2))
```

**Análisis visual:**
- ¿Qué estrategia reduce más la dispersión (altura de las cajas)?
- ¿Hay trade-offs entre reducir incertidumbre y costo esperado?
- ¿Cuál sería tu recomendación al Director?

---

# Parte 9: Tabla de Decisión para el Director

```{webr-r}
cat("\n╔═══════════════════════════════════════════════════════════╗\n")
cat("║  MEMO PARA DECISIÓN: PRESUPUESTO PROGRAMA TMC 2026       ║\n")
cat("╚═══════════════════════════════════════════════════════════╝\n\n")

cat("DE:       Dirección de Análisis Cuantitativo\n")
cat("PARA:     Director de Programas Sociales\n")
cat("FECHA:    ", format(Sys.Date(), "%d de %B, %Y"), "\n")
cat("ASUNTO:   Recomendación presupuestaria basada en simulación Monte Carlo\n\n")
cat("───────────────────────────────────────────────────────────\n\n")

cat("RESUMEN EJECUTIVO\n\n")

cat("Basado en 10,000 simulaciones Monte Carlo considerando:\n")
cat("  • Variabilidad en costo unitario (inflación, logística)\n")
cat("  • Incertidumbre en número de beneficiarios (demanda)\n")
cat("  • Fluctuación en costos logísticos regionales\n")
cat("  • Correlación FUERTE entre beneficiarios y costos operativos\n\n")

cat("┌─────────────────────────────────────────────────────────┐\n")
cat("│ OPCIONES DE PRESUPUESTO                                 │\n")
cat("└─────────────────────────────────────────────────────────┘\n\n")

cat("Opción 1: CONSERVADORA (P90)\n")
cat("  Monto:       $", format(round(percentiles[6], 0), big.mark = ","), 
    " USD\n")
cat("  Riesgo:      Solo 10% de escenarios rebasan este monto\n")
cat("  Beneficio:   Evita presupuesto suplementario\n")
cat("  Costo:       Podría sobrar presupuesto en escenarios favorables\n\n")

cat("Opción 2: PRUDENTE (P75)\n")
cat("  Monto:       $", format(round(percentiles[5], 0), big.mark = ","), 
    " USD\n")
cat("  Riesgo:      25% de escenarios rebasan este monto\n")
cat("  Beneficio:   Balance entre prudencia y eficiencia\n")
cat("  Costo:       Podría requerir ajuste presupuestal mid-year\n\n")

cat("Opción 3: MEDIANA (P50)\n")
cat("  Monto:       $", format(round(percentiles[4], 0), big.mark = ","), 
    " USD\n")
cat("  Riesgo:      50% de escenarios rebasan este monto\n")
cat("  Beneficio:   Uso eficiente de recursos\n")
cat("  Costo:       Alta probabilidad de necesitar suplemento\n\n")

cat("┌─────────────────────────────────────────────────────────┐\n")
cat("│ RECOMENDACIÓN FINAL                                     │\n")
cat("└─────────────────────────────────────────────────────────┘\n\n")

cat("Solicitar: $", format(round(percentiles[6], 0), big.mark = ","), 
    " USD (P90)\n\n")

cat("Justificación:\n")
cat("  1. Minimiza riesgo político de presupuesto insuficiente\n")
cat("  2. Costo de suplemento > Costo de sobre-estimar\n")
cat("  3. Excedentes pueden reasignarse dentro del sector\n")
cat("  4. Alta visibilidad del programa requiere conservadurismo\n\n")

cat("ADEMÁS, IMPLEMENTAR:\n\n")
cat("  Estrategia Combinada (Contratos fijos + Pre-registro):\n")
cat("  • Contratos a precio fijo con proveedores principales\n")
cat("  • Pre-registro obligatorio de beneficiarios\n")
cat("  → Reducción esperada de incertidumbre (IQR): ", 
    round(reduccion_iqr_c, 1), "%\n")
cat("  → Nuevo P90 estimado con estrategia: $",
    format(round(quantile(presup_estrategia_c, 0.90), 0), big.mark = ","), "\n\n")

cat("LIMITACIONES DEL MODELO:\n\n")
cat("    Este modelo NO considera:\n")
cat("      • Shocks macroeconómicos extremos (crisis financiera, pandemia)\n")
cat("      • Cambios regulatorios durante ejecución\n")
cat("      • Variación geográfica en costos (urbano vs rural)\n")
cat("      • Economías de escala con más beneficiarios\n")
cat("      • Eventos políticos que afecten demanda\n\n")

cat("  ⏰  Actualización recomendada:\n")
cat("      • Re-simular trimestralmente con datos reales observados\n")
cat("      • Ajustar distribuciones basándose en ejecución\n\n")

cat("═══════════════════════════════════════════════════════════\n")
```

---

# Parte 10: Simulación de Escenarios Extremos

```{webr-r}
# ¿Qué pasa en el peor escenario razonable?
cat("\n ANÁLISIS DE ESCENARIOS EXTREMOS\n")
cat("═══════════════════════════════════════════════════════════\n\n")

# Peor 5% de escenarios
peor_5pct <- presupuesto_total >= quantile(presupuesto_total, 0.95)
presup_extremo <- presupuesto_total[peor_5pct]

cat("ESCENARIO PESIMISTA (Peor 5%):\n")
cat("  Presupuesto mínimo:  $", 
    format(round(min(presup_extremo), 0), big.mark = ","), "\n")
cat("  Presupuesto máximo:  $", 
    format(round(max(presup_extremo), 0), big.mark = ","), "\n")
cat("  Promedio en peor 5%: $", 
    format(round(mean(presup_extremo), 0), big.mark = ","), "\n\n")

# ¿Qué combinación lleva al peor escenario?
idx_peor <- which.max(presupuesto_total)
cat("EL PEOR ESCENARIO ÚNICO:\n")
cat("  Presupuesto total:   $", 
    format(round(presupuesto_total[idx_peor], 0), big.mark = ","), "\n")
cat("  Costo unitario:      $", round(costo_unitario[idx_peor], 2), "\n")
cat("  Beneficiarios:       ", n_beneficiarios[idx_peor], "\n")
cat("  Logística:           $", 
    format(round(costos_logistica[idx_peor], 0), big.mark = ","), "\n\n")

# Mejor escenario
idx_mejor <- which.min(presupuesto_total)
cat("EL MEJOR ESCENARIO:\n")
cat("  Presupuesto total:   $", 
    format(round(presupuesto_total[idx_mejor], 0), big.mark = ","), "\n")
cat("  Costo unitario:      $", round(costo_unitario[idx_mejor], 2), "\n")
cat("  Beneficiarios:       ", n_beneficiarios[idx_mejor], "\n")
cat("  Logística:           $", 
    format(round(costos_logistica[idx_mejor], 0), big.mark = ","), "\n\n")

cat("RANGO TOTAL (min-max):  $", 
    format(round(max(presupuesto_total) - min(presupuesto_total), 0), 
           big.mark = ","), "\n\n")

cat(" INSIGHT: Los escenarios extremos ocurren cuando MÚLTIPLES variables\n")
cat("   se alinean en la misma dirección (todas altas o todas bajas).\n")
cat("   La FUERTE correlación entre beneficiarios y logística amplifica esto.\n\n")

cat("═══════════════════════════════════════════════════════════\n")
```

---

# Parte 11: Errores Comunes en Simulación Monte Carlo

```{webr-r}
cat("\n🚫 TOP 5 ERRORES QUE EVITAR EN SIMULACIÓN MONTE CARLO\n")
cat("═══════════════════════════════════════════════════════════\n\n")

cat("1. USAR DISTRIBUCIONES SIN JUSTIFICACIÓN\n")
cat("   ❌ Malo:  'Usé Normal porque es la que conozco'\n")
cat("    Bueno: 'Usé Lognormal porque la variable es siempre\n")
cat("             positiva y tiene cola derecha por inflación'\n\n")

cat("2. IGNORAR CORRELACIONES ENTRE VARIABLES\n")
cat("   ❌ Malo:  Asumir que costo unitario y beneficiarios\n")
cat("             son independientes\n")
cat("    Bueno: Modelar que más beneficiarios → proporcionalmente\n")
cat("             más logística (diferencia ~15% en P90!)\n\n")

cat("3. CONFUNDIR INTERPRETACIÓN DE PERCENTILES\n")
cat("   ❌ Malo:  'P90 significa que 90% de los casos\n")
cat("             cuestan esto'\n")
cat("    Bueno: 'P90 significa que solo 10% de los casos\n")
cat("             cuestan MÁS que esto'\n\n")

cat("4. NO VALIDAR SUPUESTOS CON DATOS\n")
cat("   ❌ Malo:  Inventar parámetros sin respaldo\n")
cat("    Bueno: Usar datos históricos, benchmarks sectoriales,\n")
cat("             o análisis de sensibilidad robusto\n\n")

cat("5. REPORTAR SOLO LA MEDIA\n")
cat("   ❌ Malo:  'El costo esperado es $70,000'\n")
cat("    Bueno: 'El costo esperado es $70,000\n")
cat("             (P25-P75: $62k-$78k; P90: $85k)'\n\n")

cat("6. SUBESTIMAR LA IMPORTANCIA DE LA VARIABILIDAD (SD)\n")
cat("   ❌ Malo:  'La media es casi igual, entonces no importa'\n")
cat("    Bueno: 'Un cambio del 15% en SD puede cambiar P90 en 5-10%\n")
cat("             especialmente en distribuciones de cola pesada'\n\n")

cat("═══════════════════════════════════════════════════════════\n\n")

cat("  CUÁNDO NO USAR MONTE CARLO:\n\n")
cat("  1. Cuando tienes ecuaciones cerradas para calcular varianza\n")
cat("  2. Cuando no hay incertidumbre real (todo es determinístico)\n")
cat("  3. Cuando no tienes datos para parametrizar: 'garbage in, garbage out'\n")
cat("  4. Cuando las correlaciones son complejas y desconocidas\n")
cat("  5. Para eventos de 'cisne negro' (cola extremadamente gorda)\n\n")
```

---

# Ejercicio Individual (Entregable)

**Parte A: Modificar Supuestos**

1. **Re-simula el modelo** con beneficiarios más variables:
   - Cambia SD de 600 a 1000
   - ¿Cómo cambia el P90?
   - ¿Cuánto más presupuesto necesitarías?

2. **Escenario de inflación:**
   - Aumenta el costo unitario promedio de $12 a $15
   - Mantén todo lo demás igual
   - Calcula el nuevo P90 y compara

3. **Análisis con datos reales (NUEVO):**
   - Revisa los parámetros ajustados vs asumidos en la Parte 2.6
   - ¿En cuál variable hubo mayor diferencia?
   - Si tuvieras que recolectar datos históricos con presupuesto limitado, ¿en cuál variable te enfocarías primero?

**Parte B: Análisis Escrito**

4. **Tabla comparativa** de las 4 estrategias:
   - Columnas: Estrategia, P50, P90, IQR, Costo de implementación
   - ¿Cuál recomendarías y por qué?

5. **Escribe un párrafo** (6 líneas) explicando:
   - ¿Por qué NO usar solo el promedio?
   - ¿Qué percentil usarías para presupuestar?
   - ¿Qué riesgo implica cada opción (P50 vs P90)?

6. **Ventajas de datos reales** (4-5 líneas):
   - ¿Qué ventajas viste al parametrizar con datos históricos vs supuestos?
   - ¿En qué situación justificarías el costo de recolectar datos primarios?
   - ¿Qué limitaciones tienen los datos históricos para predecir el futuro?

**Parte C: Comunicación**

7. **Crea un memo de 1 página** para el Director que incluya:
   - Recomendación de monto (con percentil)
   - 1 gráfico (histograma o boxplot)
   - 2-3 supuestos clave
   - 2 limitaciones del modelo
   - **NUEVO**: Una nota sobre la calidad de los parámetros (datos reales vs supuestos)

---

# Preguntas de Discusión

1. **Planificación Conservadora vs Eficiente:**
   - ¿Qué costo tiene planificar con P90 vs P50?
   - ¿En qué contextos políticos usarías cada uno?
   - ¿Cómo cambiaría tu decisión si fuera un programa piloto vs un servicio de salud crítico?

2. **Comunicación de Incertidumbre:**
   - ¿Cómo le explicas al Congreso "90% de confianza"?
   - ¿Preferirías dar un rango o un número único?
   - ¿Qué reacción política esperas si pides presupuesto conservador?

3. **Mejora del Modelo:**
   - ¿Qué otras fuentes de incertidumbre agregarías?
   - ¿Cómo modelarías economías de escala?
   - ¿Qué distribución usarías para eventos de "cisne negro"?

4. **Ética de la Incertidumbre:**
   - ¿Es ético pedir presupuesto conservador sabiendo que podría sobrar?
   - ¿Preferirías quedarte corto o sobre-estimar?
   - ¿Qué responsabilidad tiene el analista de comunicar riesgo honestamente?

5. **Importancia de las Correlaciones:**
   - Vimos que ignorar correlaciones subestima el P90 en ~15%
   - ¿En qué otros contextos de política pública son críticas las correlaciones?
   - ¿Cómo identificarías correlaciones importantes en tu propio trabajo?

6. **Datos Reales vs Supuestos Expertos:**
   - ¿Cuándo confiarías más en el juicio de expertos que en datos históricos?
   - Si los datos muestran algo contra-intuitivo, ¿cómo resolverías el conflicto?
   - ¿Qué harías si solo tienes datos de 5 municipios en lugar de 50?
   - ¿Cómo balanceas el costo de recolectar datos vs el riesgo de usar supuestos incorrectos?

---

# Recursos Adicionales

**Conceptos clave:**

- **Simulación Monte Carlo:** Muestreo aleatorio repetido para generar distribución de resultados
- **Distribuciones:**
  - Normal: Suma de factores independientes
  - Log-normal: Variables siempre positivas con cola derecha
  - Gamma: Tiempo hasta eventos, costos agregados
  - Triangular: Estimaciones de expertos (min/modo/max)
  - Beta: Variables acotadas entre 0 y 1
- **Parametrización empírica:** Ajustar distribuciones a datos históricos usando:
  - Método de momentos (igualar media y varianza observadas)
  - Máxima verosimilitud (MLE - maximizar probabilidad de observar los datos)
  - Validación visual (Q-Q plots, histogramas superpuestos)
- **Percentiles de riesgo:**
  - P90: Conservador (10% rebasa)
  - P75: Prudente (25% rebasa)
  - P50: Neutral (50% rebasa)
- **Análisis de sensibilidad:** Identificar qué input genera más varianza en output
- **Correlaciones:** Modelar dependencias entre variables (CRÍTICO - puede cambiar resultados 15%+)

**Paquetes de R útiles para ajuste de distribuciones:**
- `fitdistrplus`: Ajuste de distribuciones con múltiples métodos
- `MASS::fitdistr()`: Ajuste básico con MLE
- `EnvStats`: Distribuciones especializadas y tests de bondad de ajuste
- `actuar`: Distribuciones para ciencias actuariales y seguros

**Usos en política pública:**

-  Planificación presupuestaria
-  Análisis costo-beneficio bajo incertidumbre
- 🏥 Proyecciones de demanda en salud
- 🌾 Evaluación de riesgo en programas agrícolas
- ⚡ Proyecciones de energía y cambio climático
- 🎓 Planificación de infraestructura educativa
-  Preparación para emergencias y desastres

**Para profundizar:**
- Busca "The Black Swan" de Nassim Taleb (eventos extremos)
- Investiga "Value at Risk" (VaR) en finanzas
- Lee sobre "Precautionary Principle" en política ambiental
- Consulta "Thinking in Bets" de Annie Duke (decisiones bajo incertidumbre)

---

**Última reflexión:** 

*"En un mundo incierto, dar un solo número es mentir con confianza. La distribución completa es la verdad que necesitan quienes toman decisiones."*

El valor de Monte Carlo no está en predecir el futuro exacto, sino en cuantificar nuestra ignorancia de manera honesta y útil para la toma de decisiones.

**Lección especial sobre correlaciones:** Como vimos, ignorar una correlación importante (beneficiarios ↔ logística) puede subestimar el presupuesto necesario en 15-20%. En política pública, estas dependencias ocultas son la diferencia entre planificar bien y fallarle a los ciudadanos.

**Lección especial sobre parametrización:** La diferencia entre usar supuestos expertos vs datos reales puede parecer pequeña en la media (0.1%), pero crecer significativamente en los percentiles altos (5-10%). Si presupuestas con P90 o P95, **los datos no son un lujo, son una necesidad**. El costo de recolectar 50 observaciones históricas es típicamente menor que el costo político de pedir presupuesto suplementario.

**El verdadero poder de Monte Carlo** está en convertir la incertidumbre en una conversación informada sobre riesgo, donde los tomadores de decisiones pueden elegir conscientemente su nivel de tolerancia al riesgo en lugar de simplemente esperar que "todo salga bien".
