--- title: "Marco de Análisis de Sensibilidad para la Desagregación Económica Bayesiana" author: "José Mauricio Gómez Julián" date: "`r Sys.Date()`" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Marco de Análisis de Sensibilidad} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- > **Cómo leer este manual.** > Las secciones 1–3 desarrollan la **teoría** (con ecuaciones); las secciones 4–6 presentan **diagnósticos** y **métricas**; las secciones 7–8 proveen **código reproducible**: una demo sintética rápida (se evalúa al compilar con *knit*) y un **pipeline con datos reales** completo (deshabilitado por defecto por velocidad; habilítalo fijando `eval=TRUE`). Todo el código es consistente con las funciones exportadas por el paquete `BayesianDisaggregation`. ```{r setup, include=TRUE} # Valores predeterminados globales de los chunks knitr::opts_chunk$set( echo = TRUE, message = FALSE, warning = FALSE, fig.width = 9, fig.height = 6 ) # Librerías suppressPackageStartupMessages({ library(BayesianDisaggregation) library(dplyr) library(tidyr) library(ggplot2) library(readr) library(openxlsx) }) # Verbosidad del registro desde el paquete log_enable("INFO") set.seed(2024) ``` # 1. Planteamiento del problema Observamos un **índice agregado** (p. ej., IPC) por período $t=1,\dots,T$, y queremos una **desagregación sectorial** en $K$ componentes cuyas participaciones por período yacen en el **símplice unitario**: $$ W_t = (w_{t1},\dots,w_{tK}),\qquad w_{tk}\ge 0,\quad \sum_{k=1}^K w_{tk}=1. $$ Partimos de una **matriz de pesos previa** $P\in\mathbb{R}^{T\times K}$ (filas en el símplice), y construimos una **verosimilitud de sectores** $L\in\Delta^{K-1}$ (un vector no negativo que suma uno). Un **perfil temporal** entonces distribuye $L$ a $LT\in\mathbb{R}^{T\times K}$. Finalmente, una **regla de actualización determinista** combina $P$ y $LT$ para obtener el posterior $W$. # 2. Construcción de la verosimilitud sectorial $L$ ## 2.1 ACP/SVD de la matriz previa centrada Sea $P$ validada (finita, no negativa, filas $\approx 1$; pequeñas desviaciones se renormalizan). **Centramos** columnas en el tiempo: $$ X = P - \mathbf{1}\,\bar p^\top,\quad \bar p = \frac{1}{T}\sum_{t=1}^T P_{t\cdot}. $$ Calculamos la SVD $X = U\Sigma V^\top$. Sea $v$ el **primer vector singular derecho** (cargas de la CP1). Mapeamos a una saliencia no negativa vía valores absolutos y normalizamos: $$ \ell_k = |v_k|,\qquad L_k = \frac{\ell_k}{\sum_j \ell_j}. $$ Si la CP1 es **degenerada** (varianza casi nula o columnas idénticas), recurrimos a las **medias de columna** de $P$ (renormalizadas). Esto está implementado en: ```{r L_from_P, include=TRUE} # Llamada de ejemplo (los internos están en el paquete): # L <- compute_L_from_P(P) ``` **Diagnósticos adjuntos a `L`:** atributos `"pc1_loadings"`, `"explained_var"` y `"fallback"`. ## 2.2 Difusión temporal de $L$ Creamos una matriz dependiente del tiempo $LT$ aplicando un perfil de pesos no negativo $w_t$ y renormalizando por filas: $$ LT_{t,k} \propto w_t L_k,\qquad \sum_k LT_{t,k}=1. $$ **Patrones** incorporados: * `constant`: $w_t=1$ * `recent`: linealmente creciente en $t$ (más peso a los períodos recientes) * `linear`: rampa afín entre extremos * `bell`: bulto simétrico tipo gaussiano alrededor de $T/2$ ```{r spread-L, include=TRUE} # Llamada de ejemplo: # LT <- spread_likelihood(L, T_periods = nrow(P), pattern = "recent") ``` # 3. Reglas de actualización posterior (deterministas, sin MCMC) Dados $P$ y $LT$ (ambos por filas en el símplice), definimos cuatro actualizaciones deterministas: * **Promedio ponderado** (parámetro de mezcla $\lambda\in[0,1]$): $$ W = \mathsf{norm}_1\{\lambda P + (1-\lambda)LT\}. $$ * **Multiplicativa** (producto elemento a elemento con renormalización): $$ W = \mathsf{norm}_1\{P\odot LT\}. $$ * **Media Dirichlet** (conjugación analítica, $\gamma>0$, menor $\gamma\Rightarrow$ más aguda): $$ \alpha_{\text{post}} = \frac{P}{\gamma} + \frac{LT}{\gamma},\qquad W = \frac{\alpha_{\text{post}}}{\mathbf{1}^\top\alpha_{\text{post}}}. $$ * **Adaptativa** (mezcla por sector según la volatilidad previa): $$ \phi_k=\min\!\Big(\frac{\sigma_k}{\bar\sigma},\,0.8\Big),\quad W_{t\cdot}=\mathsf{norm}_1\{(1-\phi)\odot P_{t\cdot} + \phi\odot LT_{t\cdot}\}. $$ Todas están expuestas en el paquete: ```{r posteriors, include=TRUE} # posterior_weighted(P, LT, lambda = 0.7) # posterior_multiplicative(P, LT) # posterior_dirichlet(P, LT, gamma = 0.1) # posterior_adaptive(P, LT) ``` # 4. Coherencia, estabilidad e interpretabilidad ## 4.1 Coherencia respecto de $L$ Definimos las **medias temporales** previa/posterior: $$ \bar p = \frac{1}{T}\sum_t P_{t\cdot},\qquad \bar w = \frac{1}{T}\sum_t W_{t\cdot}. $$ Sea $\rho(\cdot,\cdot)$ una **correlación robusta** (máximo de |Pearson| y |Spearman|). La **coherencia** escala el **incremento** $\Delta\rho=\max(0,\rho(\bar w,L)-\rho(\bar p,L))$: $$ \text{coherence}=\min\{1,\ \text{const} + \text{mult}\cdot\Delta\rho\}. $$ ## 4.2 Estabilidad numérica y temporal * **Estabilidad numérica (penalización exponencial)** sobre desviación de suma por fila y negativos: $$ S_{\text{num}}=\exp\{-a\cdot\overline{|\sum_k W_{tk}-1|} - b\cdot \#(W<0)\}. $$ * **Estabilidad temporal** vía $|\Delta|$ promedio (menor variación $\Rightarrow$ mayor puntaje): $$ S_{\text{tmp}} = \frac{1}{1+\kappa\cdot \overline{|\Delta W|}},\quad \overline{|\Delta W|}=\frac{1}{K}\sum_k \frac{1}{T-1}\sum_{t}|W_{t+1,k}-W_{t,k}|. $$ * **Estabilidad compuesta** (pesos predeterminados 60% numérica, 40% temporal): $$ S_{\text{comp}} = 0.6\,S_{\text{num}} + 0.4\,S_{\text{tmp}}. $$ Funciones del paquete: ```{r metrics-fns, include=TRUE} # coherence_score(P, W, L, mult = 3.0, const = 0.5) # numerical_stability_exp(W, a = 1000, b = 10) # temporal_stability(W, kappa = 50) # stability_composite(W, a = 1000, b = 10, kappa = 50) ``` ## 4.3 Interpretabilidad Dos principios: 1. **Preservación** de la estructura sectorial (correlación entre $\bar p$ y $\bar w$); 2. **Plausibilidad** de cambios sectoriales promedio (penalizar desplazamientos relativos extremos). Implementación: $$ \text{pres}=\max\{0,\rho(\bar p,\bar w)\},\qquad r_k=\frac{|\,\bar w_k-\bar p_k\,|}{\bar p_k+\epsilon},\quad \text{plaus}= \frac{1}{1+2\cdot Q_{0.9}(r_k)}. $$ Luego $\text{interp}=0.6\,\text{pres}+0.4\,\text{plaus}$. ```{r interp-fn, include=TRUE} # interpretability_score(P, W, use_q90 = TRUE) ``` # 5. API de punta a punta (`bayesian_disaggregate`) El pipeline de conveniencia: 1. `read_cpi()` y `read_weights_matrix()` (Excel) 2. `compute_L_from_P(P)` y `spread_likelihood(L, T, pattern)` 3. regla posterior (`weighted` / `multiplicative` / `dirichlet` / `adaptive`) 4. métricas: coherencia, estabilidad (compuesta), interpretabilidad, eficiencia (heurística), puntaje compuesto 5. *helpers* de exportación: `save_results()` y un libro de Excel de un solo archivo para la configuración “mejor” ```{r api, include=TRUE} # Firma de ejemplo (ver Sección 8 para datos reales): # bayesian_disaggregate(path_cpi, path_weights, # method = c("weighted","multiplicative","dirichlet","adaptive"), # lambda = 0.7, gamma = 0.1, # coh_mult = 3.0, coh_const = 0.5, # stab_a = 1000, stab_b = 10, stab_kappa = 50, # likelihood_pattern = "recent") ``` # 6. Interpretación de visualizaciones clave * **Mapa de calor del posterior $W$**: cada **celda** es la participación de un sector en un año; las **filas** son años, las **columnas** son sectores. *Léelo así*: celdas más oscuras = mayor participación sectorial; la **suavidad horizontal** indica estabilidad temporal; los **patrones verticales** (bandas) muestran importancia sectorial persistente. * **Líneas de sectores principales**: para los sectores más relevantes por participación media, **líneas** que siguen la participación del sector en el tiempo. *Léelo así*: niveles consistentes = estabilidad; cambios de tendencia coinciden con cambios en la estructura macro. * **Hoja de IPC sectorial**: $\hat Y_{t,k} = \text{CPI}_t \times W_{t,k}$. *Léelo así*: descomposición dolarizada (o escalada en índice) del agregado. # 7. Demo sintética reproducible (se evalúa al compilar) Este *chunk* sintetiza un ejemplo pequeño que puedes compilar con seguridad. ```{r demo, include=TRUE} # Matriz previa sintética (filas en el símplice) T <- 10; K <- 6 set.seed(123) P <- matrix(rexp(T*K), nrow = T) P <- P / rowSums(P) # Vector de verosimilitud desde P (ACP/SVD; robusto con alternativa) L <- compute_L_from_P(P) # Difundir en el tiempo con patrón "recent" LT <- spread_likelihood(L, T_periods = T, pattern = "recent") # Probar un par de posteriores W_weighted <- posterior_weighted(P, LT, lambda = 0.7) W_adaptive <- posterior_adaptive(P, LT) # Métricas para el adaptativo coh <- coherence_score(P, W_adaptive, L) stab <- stability_composite(W_adaptive, a = 1000, b = 10, kappa = 50) intr <- interpretability_score(P, W_adaptive) eff <- 0.65 comp <- 0.30*coh + 0.25*stab + 0.25*intr + 0.20*eff data.frame(coherence = coh, stability = stab, interpretability = intr, efficiency = eff, composite = comp) %>% round(4) ``` # 8. Pipeline completo con datos reales (habilitar/deshabilitar evaluación) > **Cambia a `eval=TRUE` tras fijar tus rutas**. Por defecto mantenemos este *chunk* apagado para renderizar rápido en cualquier máquina. ```{r real-pipeline, eval=FALSE} # === Crear datos sintéticos para demo compatible con CRAN === demo_dir <- tempdir() # Crear datos sintéticos de IPC (imitando tu estructura) set.seed(123) cpi_demo <- data.frame( Year = 2000:2010, CPI = cumsum(c(100, rnorm(10, 0.5, 2))) # Caminata aleatoria iniciando en 100 ) cpi_file <- file.path(demo_dir, "synthetic_cpi.xlsx") openxlsx::write.xlsx(cpi_demo, cpi_file) # Crear matriz de pesos sintética (imitando estructura de pesos de VAB) set.seed(456) years <- 2000:2010 sectors <- c("Agriculture", "Manufacturing", "Services", "Construction", "Mining") weights_demo <- data.frame(Year = years) for(sector in sectors) { weights_demo[[sector]] <- runif(length(years), 0.05, 0.35) } # Normalizar filas para sumar 1 (restricción de símplice) weights_demo[, -1] <- weights_demo[, -1] / rowSums(weights_demo[, -1]) weights_file <- file.path(demo_dir, "synthetic_weights.xlsx") openxlsx::write.xlsx(weights_demo, weights_file) # Usar rutas de datos sintéticos path_cpi <- cpi_file path_w <- weights_file out_dir <- demo_dir cat("Usando datos sintéticos para la demo de CRAN:\n") cat("Archivo CPI:", path_cpi, "\n") cat("Archivo de pesos:", path_w, "\n") cat("Directorio de salida:", out_dir, "\n") # --- Ejecución base (predeterminados robustos) --- base_res <- bayesian_disaggregate( path_cpi = path_cpi, path_weights = path_w, method = "adaptive", lambda = 0.7, # registrado en métricas; no usado por "adaptive" gamma = 0.1, coh_mult = 3.0, coh_const = 0.5, stab_a = 1000, stab_b = 10, stab_kappa = 60, likelihood_pattern = "recent" ) xlsx_base <- save_results(base_res, out_dir = file.path(out_dir, "base")) print(base_res$metrics) # --- Búsqueda de cuadrícula mínima para la demo (tamaño reducido) --- n_cores <- 1 # Un solo núcleo para cumplimiento CRAN grid_df <- expand.grid( method = c("weighted", "adaptive"), # Métodos reducidos lambda = c(0.5, 0.7), # Opciones reducidas gamma = 0.1, # Opción única coh_mult = 3.0, # Opción única coh_const = 0.5, # Opción única stab_a = 1000, stab_b = 10, stab_kappa = 60, # Opción única likelihood_pattern = "recent", # Opción única KEEP.OUT.ATTRS = FALSE, stringsAsFactors = FALSE ) grid_res <- run_grid_search( path_cpi = path_cpi, path_weights = path_w, grid_df = grid_df, n_cores = n_cores ) write.csv(grid_res, file.path(out_dir, "grid_results.csv"), row.names = FALSE) best_row <- grid_res %>% arrange(desc(composite)) %>% slice(1) print("Mejor configuración de la búsqueda en cuadrícula:") print(best_row) # --- Re-ejecutar la mejor configuración para exportación limpia --- best_res <- bayesian_disaggregate( path_cpi = path_cpi, path_weights = path_w, method = best_row$method, lambda = if (!is.na(best_row$lambda)) best_row$lambda else 0.7, gamma = if (!is.na(best_row$gamma)) best_row$gamma else 0.1, coh_mult = best_row$coh_mult, coh_const = best_row$coh_const, stab_a = best_row$stab_a, stab_b = best_row$stab_b, stab_kappa = best_row$stab_kappa, likelihood_pattern = best_row$likelihood_pattern ) xlsx_best <- save_results(best_res, out_dir = file.path(out_dir, "best")) # --- Un Excel con todo (incluyendo hiperparámetros) --- sector_summary <- tibble( Sector = colnames(best_res$posterior)[-1], prior_mean = colMeans(as.matrix(best_res$prior[, -1])), posterior_mean = colMeans(as.matrix(best_res$posterior[, -1])) ) wb <- createWorkbook() addWorksheet(wb, "Hyperparameters"); writeData(wb, "Hyperparameters", best_row) addWorksheet(wb, "Metrics"); writeData(wb, "Metrics", best_res$metrics) addWorksheet(wb, "Prior_P"); writeData(wb, "Prior_P", best_res$prior) addWorksheet(wb, "Posterior_W"); writeData(wb, "Posterior_W", best_res$posterior) addWorksheet(wb, "Likelihood_t"); writeData(wb, "Likelihood_t", best_res$likelihood_t) addWorksheet(wb, "Likelihood_L"); writeData(wb, "Likelihood_L", best_res$likelihood) addWorksheet(wb, "Sector_Summary"); writeData(wb, "Sector_Summary", sector_summary) for (sh in c("Hyperparameters","Metrics","Prior_P","Posterior_W", "Likelihood_t","Likelihood_L","Sector_Summary")) { freezePane(wb, sh, firstRow = TRUE) addFilter(wb, sh, rows = 1, cols = 1:ncol(readWorkbook(wb, sh))) setColWidths(wb, sh, cols = 1:200, widths = "auto") } # --- Añadir IPC sectorial (agregado por pesos posteriores) --- W_post <- best_res$posterior # Year + sectores cpi_df <- read_cpi(path_cpi) # Year, CPI sector_cpi <- dplyr::left_join(W_post, cpi_df, by = "Year") %>% dplyr::mutate(dplyr::across(-c(Year, CPI), ~ .x * CPI)) # Verificación de calidad: suma de sectores vs CPI check_sum <- sector_cpi %>% dplyr::mutate(row_sum = rowSums(dplyr::across(-c(Year, CPI))), diff = CPI - row_sum) cat("Verificación de calidad (primeras 5 filas):\n") print(head(check_sum, 5)) addWorksheet(wb, "Sector_CPI") writeData(wb, "Sector_CPI", sector_cpi) freezePane(wb, "Sector_CPI", firstRow = TRUE) addFilter(wb, "Sector_CPI", rows = 1, cols = 1:ncol(sector_cpi)) setColWidths(wb, "Sector_CPI", cols = 1:200, widths = "auto") excel_onefile <- file.path(out_dir, "best", "Best_Full_Output_withSectorCPI.xlsx") saveWorkbook(wb, excel_onefile, overwrite = TRUE) cat("Resultados completos guardados en:", excel_onefile, "\n") # --- Gráficos rápidos (guardados como PNGs) --- dir_plots <- file.path(out_dir, "best", "plots") if (!dir.exists(dir_plots)) dir.create(dir_plots, recursive = TRUE) W_long <- best_res$posterior %>% pivot_longer(-Year, names_to = "Sector", values_to = "Weight") p_heat <- ggplot(W_long, aes(Year, Sector, fill = Weight)) + geom_tile() + scale_fill_viridis_c() + labs(title = "Pesos posteriores (W): mapa de calor", x = "Año", y = "Sector", fill = "Participación") + theme_minimal(base_size = 11) + theme(axis.text.y = element_text(size = 6)) ggsave(file.path(dir_plots, "posterior_heatmap.png"), p_heat, width = 12, height = 9, dpi = 220) top_sectors <- best_res$posterior %>% summarise(across(-Year, mean)) %>% pivot_longer(everything(), names_to = "Sector", values_to = "MeanShare") %>% arrange(desc(MeanShare)) %>% slice(1:3) %>% pull(Sector) # Reducido a top 3 para la demo p_lines <- best_res$posterior %>% select(Year, all_of(top_sectors)) %>% pivot_longer(-Year, names_to = "Sector", values_to = "Weight") %>% ggplot(aes(Year, Weight, color = Sector)) + geom_line(linewidth = 0.9) + labs(title = "Top 3 sectores por participación media (posterior W)", y = "Participación", x = "Año") + theme_minimal(base_size = 11) ggsave(file.path(dir_plots, "posterior_topSectors.png"), p_lines, width = 11, height = 6, dpi = 220) cat("Demo completada exitosamente. Todos los archivos escritos al directorio temporal.\n") ``` # 9. Guía práctica y valores predeterminados * Prefiere `method="adaptive"` cuando las volatilidades sectoriales previas son heterogéneas; de lo contrario, `weighted` con $\lambda\in[0.7,0.9]$ es sólida y suele encabezar la cuadrícula. * Los parámetros predeterminados de **coherencia** `(mult=3.0, const=0.5)` producen un puntaje acotado e interpretable 0–1 que enfatiza la **mejora** sobre la previa. * La penalización **exponencial** numérica es intencionalmente pronunciada: mantiene a raya las desviaciones de suma por fila y los negativos en corridas automatizadas y búsquedas en cuadrícula. * Para informes, exporta **Sector_CPI** para ilustrar la descomposición económica $\hat Y_{t,k}$. # Apéndice A. Invariantes y comprobaciones rápidas ```{r invariants, include=TRUE} # Ejemplo: invariantes en una corrida sintética fresca T <- 6; K <- 5 set.seed(7) P <- matrix(rexp(T*K), nrow = T); P <- P / rowSums(P) L <- compute_L_from_P(P) LT <- spread_likelihood(L, T, "recent") W <- posterior_multiplicative(P, LT) # Invariantes stopifnot(all(abs(rowSums(P) - 1) < 1e-12)) stopifnot(all(abs(rowSums(LT) - 1) < 1e-12)) stopifnot(all(abs(rowSums(W) - 1) < 1e-12)) c( coherence = coherence_score(P, W, L), stability = stability_composite(W), interpret = interpretability_score(P, W) ) %>% round(4) ``` # Apéndice B. Información de la sesión ```{r session, include=TRUE} sessionInfo() ``` **Notas** - El *chunk* de *datos reales* está con `eval=FALSE` para que el manual se renderice en cualquier lugar. Cámbialo a `TRUE` en tu máquina para ejecutarlo completamente contra tus archivos de Excel. - La exportación “de un solo archivo, la mejor” incluye **Hyperparameters, Metrics, Prior_P, Posterior_W, Likelihood_t, Likelihood_L, Sector_Summary, Sector_CPI**, con encabezados congelados y filtros para análisis rápido. - Los gráficos se escriben en `.../best/plots/` y coinciden con la guía de interpretación de la Sección 6.