The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

Marco de Análisis de Sensibilidad para la Desagregación Económica Bayesiana

José Mauricio Gómez Julián

2025-10-09

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:

# 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\)

# 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:

# 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:

# 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:

Preservación de la estructura sectorial (correlación entre \(\bar p\) y \(\bar w\));
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}\).

# interpretability_score(P, W, use_q90 = TRUE)

5. API de punta a punta (`bayesian_disaggregate`)

El pipeline de conveniencia:

read_cpi() y read_weights_matrix() (Excel)
compute_L_from_P(P) y spread_likelihood(L, T, pattern)
regla posterior (weighted / multiplicative / dirichlet / adaptive)
métricas: coherencia, estabilidad (compuesta), interpretabilidad, eficiencia (heurística), puntaje compuesto
helpers de exportación: save_results() y un libro de Excel de un solo archivo para la configuración “mejor”

# 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.

# 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)

##     coherence stability interpretability efficiency composite
## 90%         1    0.7537           0.6887       0.65    0.7906

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.

# === 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

# 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)

##     coherence     stability interpret.90% 
##        1.0000        0.6459        0.6245

Apéndice B. Información de la sesión

sessionInfo()

## R version 4.4.2 (2024-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26100)
## 
## Matrix products: default
## 
## 
## locale:
## [1] LC_COLLATE=C                   LC_CTYPE=Spanish_Spain.utf8   
## [3] LC_MONETARY=Spanish_Spain.utf8 LC_NUMERIC=C                  
## [5] LC_TIME=Spanish_Spain.utf8    
## 
## time zone: America/Costa_Rica
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] openxlsx_4.2.8               readr_2.1.5                 
## [3] ggplot2_4.0.0                tidyr_1.3.1                 
## [5] dplyr_1.1.4                  BayesianDisaggregation_0.1.2
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.6       jsonlite_2.0.0     compiler_4.4.2     tidyselect_1.2.1  
##  [5] Rcpp_1.1.0         zip_2.3.3          jquerylib_0.1.4    scales_1.4.0      
##  [9] yaml_2.3.10        fastmap_1.2.0      R6_2.6.1           generics_0.1.4    
## [13] knitr_1.50         iterators_1.0.14   tibble_3.3.0       tzdb_0.5.0        
## [17] RColorBrewer_1.1-3 bslib_0.9.0        pillar_1.11.0      rlang_1.1.5       
## [21] cachem_1.1.0       stringi_1.8.7      xfun_0.53          S7_0.2.0          
## [25] sass_0.4.10        cli_3.6.5          withr_3.0.2        magrittr_2.0.4    
## [29] digest_0.6.37      foreach_1.5.2      grid_4.4.2         rstudioapi_0.17.1 
## [33] hms_1.1.3          lifecycle_1.0.4    vctrs_0.6.5        evaluate_1.0.5    
## [37] glue_1.8.0         farver_2.1.2       codetools_0.2-20   rmarkdown_2.29    
## [41] purrr_1.1.0        tools_4.4.2        pkgconfig_2.0.3    htmltools_0.5.8.1

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.

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.