Documento auto-contenido para RStudio. Knit a HTML (y luego podés publicarlo en RPubs).
Todos los gráficos usan base R (sin paquetes externos).
Comparamos dos formas de agregar crecimientos sectoriales:
El “verdadero” para evaluar el error será la
suma de cantidades \(\sum_i
q_{i,t}\) (normalizada con base 100).
Medimos desempeño con RMSE de \(\Delta\log\) (crecimientos
logarítmicos).
Hechos clave: 1. Sin drift (precios relativos constantes): fijo = encadenado (mismos pesos y mismos crecimientos), aun con quiebres en cantidades. 2. Con drift relativo: el encadenado se re-pondera con \(p_{t-1}\). Según cómo ese drift se alinee con la evolución de cantidades, puede ganar o perder frente al fijo. 3. Inflación común alta (\(\mu\)) diluye shocks de precios relativos: el ratio \(\frac{1+\mu+\Delta}{1+\mu}=1+\frac{\Delta}{1+\mu}\) se acerca a 1 cuando \(\mu\) crece.
# ---------- Pesos (ponderadores de crecimiento) ----------
weights_fixed_growth <- function(Q, P){
Tt <- ncol(Q); S <- nrow(Q)
W <- matrix(NA_real_, S, Tt-1)
rownames(W) <- rownames(Q); colnames(W) <- colnames(Q)[-1]
p0 <- P[,1]
for(t in 2:Tt){
V <- Q[,t-1] * p0
W[,t-1] <- V / sum(V)
}
W
}
weights_chain_growth <- function(Q, P){
Tt <- ncol(Q); S <- nrow(Q)
W <- matrix(NA_real_, S, Tt-1)
rownames(W) <- rownames(Q); colnames(W) <- colnames(Q)[-1]
for(t in 2:Tt){
V <- Q[,t-1] * P[,t-1]
W[,t-1] <- V / sum(V)
}
W
}
# ---------- Índices y error (vs “verdadero” = suma de q) ----------
idx_volumen_verdadero <- function(Q){ base <- sum(Q[,1]); 100 * colSums(Q) / base }
idx_laspeyres_fijo <- function(Q,P){
p0 <- P[,1]; base <- sum(Q[,1]*p0); 100 * colSums(Q*p0) / base
}
idx_laspeyres_encadenado <- function(Q,P){
Tt <- ncol(Q); idx <- numeric(Tt); idx[1] <- 100
for(t in 2:Tt){
num <- sum(Q[,t] * P[,t-1])
den <- sum(Q[,t-1] * P[,t-1])
idx[t] <- idx[t-1] * (num/den)
}
idx
}
rmse_diff_log_growth <- function(idxA, idxB){
gA <- diff(log(idxA)); gB <- diff(log(idxB)); sqrt(mean((gA - gB)^2))
}
# ---------- Barras apiladas (varios años) ----------
barras_anios <- function(W, anios, titulo="Ponderadores (t-1)",
dec=0, umbral=0.05){
stopifnot(!is.null(colnames(W)))
if (is.null(rownames(W))) rownames(W) <- paste0("Sector", seq_len(nrow(W)))
cols_chr <- trimws(as.character(colnames(W)))
anios_chr <- trimws(as.character(anios))
idx <- match(anios_chr, cols_chr)
if (all(is.na(idx))) stop("Ninguno de los años solicitados existe.")
if (any(is.na(idx))) warning("Años inexistentes: ",
paste(anios_chr[is.na(idx)], collapse=", "),
". Se grafican los disponibles.")
idx <- idx[!is.na(idx)]
M <- W[, idx, drop=FALSE]
pal <- setNames(gray.colors(nrow(M)), rownames(M))
bp <- barplot(M, beside=FALSE, ylim=c(0,1), col=pal[rownames(M)],
main=titulo, xlab="Año t", ylab="Ponderador (suma=1)",
names.arg=cols_chr[idx], las=1)
legend("topright", legend=rownames(M), fill=pal[rownames(M)], bty="n")
for (j in seq_len(ncol(M))){
cum <- cumsum(M[,j]); base <- cum - M[,j]; mid <- base + M[,j]/2
show <- which(M[,j] >= umbral)
if (length(show)){
lbl <- paste0(round(100*M[show,j], dec), "%")
text(x=bp[j], y=mid[show], labels=lbl, cex=0.9)
}
}
invisible(bp)
}
# ---------- Barras apiladas (un solo año) ----------
barras_ponderadores <- function(W, anio, titulo=""){
cols <- trimws(colnames(W)); anio <- as.character(anio)
j <- match(anio, cols); if(is.na(j)) stop("Año no encontrado en W")
M <- matrix(W[,j], ncol=1); rownames(M) <- rownames(W)
pal <- gray.colors(nrow(M))
bp <- barplot(M, beside=FALSE, ylim=c(0,1), col=pal,
main=titulo, names.arg=anio, xlab="Año t",
ylab="Ponderador (suma=1)", las=1)
legend("topright", legend=rownames(M), fill=pal, bty="n")
cum <- cumsum(M[,1]); base <- cum - M[,1]; mid <- base + M[,1]/2
lbl <- paste0(round(100*M[,1],0), "%"); text(x=bp, y=mid, labels=lbl, cex=0.9)
}
# ---------- Simulador Monte Carlo con ruido proporcional a mu ----------
simular_QP_mc <- function(
S = 3,
years = 2022:2031,
qty0 = c(100,200,300),
price0 = c(1000,50,30),
mu_qty = c(0.03,0.02,0.025),
sd_qty = c(0.03,0.03,0.03),
mu_inflacion = 1.50, # 150% anual
rel_price_drift = c(0,0,0), # "sin drift" esperado
sd_factor = 0.25, # sd_price = sd_factor * mu_inflacion
break_from = NA # NA = sin quiebre
){
Tt <- length(years)
Q <- matrix(NA_real_, S, Tt); P <- matrix(NA_real_, S, Tt)
rownames(Q) <- rownames(P) <- c("Comp","Mesas","Sillas")
colnames(Q) <- colnames(P) <- years
Q[,1] <- qty0; P[,1] <- price0
for(t in 2:Tt){
gq <- rnorm(S, mu_qty, sd_qty)
if (!is.na(break_from) && years[t] >= break_from) gq[1] <- gq[1] + 0.12
Q[,t] <- pmax(Q[,t-1]*(1+gq), .Machine$double.eps)
sd_price <- rep(sd_factor*mu_inflacion, S)
gp <- rnorm(S, mu_inflacion + rel_price_drift, sd_price)
P[,t] <- pmax(P[,t-1]*(1+gp), .Machine$double.eps)
}
list(Q=Q, P=P)
}
# ---------- Promedio de ponderadores y una corrida 'ganadora' del encadenado ----------
mc_ponderadores <- function(R=300, years=2022:2031, mu_inf=1.50, sd_factor=0.25,
rel_price_drift=c(0,0,0), break_from=NA, seed=123){
set.seed(seed)
S <- 3; Tt <- length(years)
Wfix_acc <- array(0, dim=c(S, Tt-1)); rownames(Wfix_acc) <- c("Comp","Mesas","Sillas"); colnames(Wfix_acc) <- years[-1]
Wch_acc <- Wfix_acc
win_best <- +Inf; best_run <- NULL
for(r in 1:R){
sim <- simular_QP_mc(S=S, years=years, mu_inflacion=mu_inf,
rel_price_drift=rel_price_drift, sd_factor=sd_factor,
break_from=break_from)
Wfix <- weights_fixed_growth(sim$Q, sim$P)
Wch <- weights_chain_growth (sim$Q, sim$P)
Wfix_acc <- Wfix_acc + Wfix; Wch_acc <- Wch_acc + Wch
V <- idx_volumen_verdadero(sim$Q)
Lf <- idx_laspeyres_fijo(sim$Q, sim$P)
Lc <- idx_laspeyres_encadenado(sim$Q, sim$P)
diff_rmse <- rmse_diff_log_growth(Lc,V) - rmse_diff_log_growth(Lf,V) # negativo => gana encadenado
if (diff_rmse < win_best){
win_best <- diff_rmse; best_run <- list(Wfix=Wfix, Wch=Wch, Q=sim$Q, P=sim$P)
}
}
Wfix_mean <- Wfix_acc / R; Wch_mean <- Wch_acc / R
list(Wfix_mean=Wfix_mean, Wch_mean=Wch_mean, best_run=best_run,
mu_inf=mu_inf, sd_factor=sd_factor, years=years)
}
# --------- Correr Monte Carlo: caso A = sin quiebre; caso B = con quiebre ----------
years <- 2022:2031
resA <- mc_ponderadores(R=400, years=years, mu_inf=1.50, sd_factor=0.25,
rel_price_drift=c(0,0,0), break_from=NA, seed=123)
resB <- mc_ponderadores(R=400, years=years, mu_inf=1.50, sd_factor=0.25,
rel_price_drift=c(0,0,0), break_from=2027, seed=456)
anios_mostrar <- c(2026, 2028, 2031)
par(mfrow=c(2,2), mar=c(4,4,3,1))
barras_anios(resA$Wfix_mean, anios_mostrar,
titulo=sprintf("MC PROMEDIO — FIJO (p0) μ=%.0f%%, sd=%.0f%% de μ",
100*resA$mu_inf, 100*resA$sd_factor))
barras_anios(resA$Wch_mean, anios_mostrar,
titulo=sprintf("MC PROMEDIO — ENCADENADO (p_{t-1}) μ=%.0f%%, sd=%.0f%% de μ",
100*resA$mu_inf, 100*resA$sd_factor))
barras_anios(resA$best_run$Wfix, anios_mostrar, titulo="CORRIDA 'GANADORA' — FIJO (p0)")
barras_anios(resA$best_run$Wch, anios_mostrar, titulo="CORRIDA 'GANADORA' — ENCADENADO (p_{t-1})")
par(mfrow=c(1,1))
par(mfrow=c(2,2), mar=c(4,4,3,1))
barras_anios(resB$Wfix_mean, anios_mostrar,
titulo=sprintf("MC PROMEDIO — FIJO (p0) con QUIEBRE μ=%.0f%%, sd=%.0f%% de μ",
100*resB$mu_inf, 100*resB$sd_factor))
barras_anios(resB$Wch_mean, anios_mostrar,
titulo=sprintf("MC PROMEDIO — ENCADENADO (p_{t-1}) con QUIEBRE μ=%.0f%%, sd=%.0f%% de μ",
100*resB$mu_inf, 100*resB$sd_factor))
barras_anios(resB$best_run$Wfix, anios_mostrar, titulo="CORRIDA 'GANADORA' — FIJO (p0) con QUIEBRE")
barras_anios(resB$best_run$Wch, anios_mostrar, titulo="CORRIDA 'GANADORA' — ENCADENADO (p_{t-1}) con QUIEBRE")
par(mfrow=c(1,1))
hacer_escenario_drift <- function(mu=0.10, delta_rel=0.25,
years=2022:2031, t_shock=2026, t_break=2027){
S <- 3; Tt <- length(years)
Q <- matrix(NA_real_, S, Tt); P <- matrix(NA_real_, S, Tt)
rownames(Q) <- rownames(P) <- c("Comp","Mesas","Sillas")
colnames(Q) <- colnames(P) <- years
Q[,1] <- c(100,200,300); P[,1] <- c(1000,50,30)
for(t in 2:Tt){
g <- c(0.02, 0.01, 0.015)
if(years[t] >= t_break) g[1] <- 0.15
Q[,t] <- Q[,t-1]*(1+g)
gp <- rep(mu, S)
if (years[t] == t_shock) gp[1] <- mu + delta_rel
P[,t] <- P[,t-1]*(1+gp)
}
list(Q=Q,P=P,years=years,t_shock=t_shock,t_break=t_break,mu=mu,delta_rel=delta_rel)
}
experimento_mu <- function(mus = c(0.10,0.50,1.00,2.00,3.00),
delta_rel=0.25, years=2022:2031, t_shock=2026, t_break=2027){
out <- data.frame(mu=mus, rmse_fix=NA_real_, rmse_ch=NA_real_)
ejemplos <- list()
for(i in seq_along(mus)){
esc <- hacer_escenario_drift(mu=mus[i], delta_rel=delta_rel, years=years,
t_shock=t_shock, t_break=t_break)
V <- idx_volumen_verdadero(esc$Q)
Lf <- idx_laspeyres_fijo(esc$Q, esc$P)
Lc <- idx_laspeyres_encadenado(esc$Q, esc$P)
out$rmse_fix[i] <- rmse_diff_log_growth(Lf, V)
out$rmse_ch[i] <- rmse_diff_log_growth(Lc, V)
ejemplos[[i]] <- esc
}
list(tabla=out, ejemplos=ejemplos)
}
years <- 2022:2031; t_shock <- 2026; t_break <- 2027; delta_rel <- 0.25
res <- experimento_mu(mus = c(0.10,0.50,1.00,2.00,3.00),
delta_rel=delta_rel, years=years,
t_shock=t_shock, t_break=t_break)
res$tabla
## mu rmse_fix rmse_ch
## 1 0.1 0.06228775 0.06417837
## 2 0.5 0.06228775 0.06373843
## 3 1.0 0.06228775 0.06341155
## 4 2.0 0.06228775 0.06306244
## 5 3.0 0.06228775 0.06287882
op <- par(mfrow=c(1,2), mar=c(4,4,3,1))
plot(res$tabla$mu*100, res$tabla$rmse_fix, type="b", pch=17,
xlab="Inflación común μ (%)", ylab="RMSE Δlog (vs suma q)",
main=sprintf("Shock relativo Δ=%.0f p.p. en %d", 100*delta_rel, t_shock))
lines(res$tabla$mu*100, res$tabla$rmse_ch, type="b", pch=19)
legend("topright", c("Laspeyres fijo (p0)","Encadenado (p_{t-1})"),
pch=c(17,19), lty=1, bty="n")
plot(res$tabla$mu*100, res$tabla$rmse_ch - res$tabla$rmse_fix, type="b", pch=19,
xlab="Inflación común μ (%)", ylab="ΔRMSE (encadenado - fijo)",
main="Diferencia de error (negativo = mejor encadenado)")
abline(h=0, lty=2)
par(op)
mu_baja <- 0.10; mu_alta <- 2.00
esc_low <- res$ejemplos[[which(res$tabla$mu==mu_baja)]]
esc_high <- res$ejemplos[[which(res$tabla$mu==mu_alta)]]
Wc_low <- weights_chain_growth(esc_low$Q, esc_low$P)
Wc_high <- weights_chain_growth(esc_high$Q, esc_high$P)
par(mfrow=c(1,2), mar=c(4,4,3,1))
barras_ponderadores(Wc_low, esc_low$t_shock+1,
titulo=sprintf("Ponderadores encadenados — μ=%.0f%% — t=%d",
100*esc_low$mu, esc_low$t_shock+1))
barras_ponderadores(Wc_high, esc_high$t_shock+1,
titulo=sprintf("Ponderadores encadenados — μ=%.0f%% — t=%d",
100*esc_high$mu, esc_high$t_shock+1))
par(mfrow=c(1,1))
comparar_series <- function(esc, titulo=""){
yrs <- esc$years
V <- idx_volumen_verdadero(esc$Q)
Lf <- idx_laspeyres_fijo(esc$Q, esc$P)
Lc <- idx_laspeyres_encadenado(esc$Q, esc$P)
e_fix <- diff(log(Lf)) - diff(log(V))
e_ch <- diff(log(Lc)) - diff(log(V))
plot(yrs[-1], e_fix, type="b", pch=17, xlab="Año", ylab="Desvío en Δlog",
main=titulo); lines(yrs[-1], e_ch, type="b", pch=19)
abline(h=0, lty=2); legend("topright", c("Fijo","Encadenado"),
pch=c(17,19), lty=1, bty="n")
}
par(mfrow=c(1,2), mar=c(4,4,3,1))
comparar_series(esc_low, sprintf("Desvío de crecimiento — μ=%.0f%%", 100*esc_low$mu))
comparar_series(esc_high, sprintf("Desvío de crecimiento — μ=%.0f%%", 100*esc_high$mu))
par(mfrow=c(1,1))