← Catálogo
· Capítulo 1
T01_mECO_Script_VD_Limitada.R
# =============================================================================
# TEMA 1 — VARIABLE DEPENDIENTE LIMITADA: SCRIPT INTERACTIVO
# =============================================================================
# Autor: Carlos de Anta Puig
# Economista · Perito Financiero
# Miembro del Colegio de Economistas de Madrid
# Miembro del Instituto Español de Analistas Financieros (IEAF)
# Profesor de Econometría y Microeconometría
# carlos@cwconsultores.com
#
# Manual de Microeconometría - Digital Reasons (2026)
# https://github.com/carlanta/MicroEconometrics
#
# Versión: 1.0
# Fecha: 2026-03-27
#
# Ejecutar bloque a bloque en RStudio.
# =============================================================================
# Autor: Carlos de Anta Puig
# Economista · Perito Financiero
# Miembro del Colegio de Economistas de Madrid
# Miembro del Instituto Español de Analistas Financieros (IEAF)
# Profesor de Econometría y Microeconometría
# carlosmaria.deanta@unir.net
#
# Curso: Microeconometría — CdeA
# Objetivo: Explorar los distintos tipos de variables dependientes limitadas
# (binarias, multinomiales, censuradas, truncadas), entender por
# qué MCO falla y familiarizarse con las distribuciones asociadas.
#
# Version: 1.1
# Fecha: 2026-03-16
#
# INSTRUCCIONES:
# 1. Ejecuta este script línea a línea con Ctrl+Enter.
# 2. Los resultados aparecerán en la CONSOLA (abajo).
# 3. Los gráficos aparecerán en el panel PLOTS (abajo-derecha).
# 4. Las líneas que empiezan por # son COMENTARIOS (no se ejecutan).
# ============================================================================
# ============================================================================
# SECCION 1: INTRODUCCION
# ============================================================================
# Funcion de pausa interactiva
pausa <- function(msg = "\n>>> Pulsa ENTER para continuar...") {
readline(msg)
}
# Carga de paquetes
suppressPackageStartupMessages({
library(DT)
})
# Funcion auxiliar: tabla interactiva con exportacion
tabla_dt <- function(datos, titulo = "", decimales = 2, filas = TRUE) {
if (is.numeric(datos) && !is.null(names(datos))) {
datos <- data.frame(Nombre = names(datos), Valor = round(datos, decimales))
filas <- FALSE
}
DT::datatable(
datos,
caption = titulo,
extensions = 'Buttons',
options = list(
dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel'),
pageLength = 15,
scrollX = TRUE
),
rownames = filas
)
}
cat("\n")
cat("============================================================\n")
cat("║ ║\n")
cat("║ TEMA 1 — Datos con Variable Dependiente Limitada ║\n")
cat("║ Script compañero del manual de apuntes ║\n")
cat("║ ║\n")
cat("║ Profesor: Carlos de Anta Puig ║\n")
cat("║ Curso 2025-2026 | Q2 ║\n")
cat("║ ║\n")
cat("╚══════════════════════════════════════════════════════════════════╝\n")
cat("\n")
pausa()
# ============================================================================
# SECCION 2: VARIABLES BINARIAS Y DISTRIBUCION DE BERNOULLI
# ============================================================================
cat("═══════════════════════════════════════════════════════════════════\n")
cat(" SECCION 2: VARIABLES BINARIAS Y DISTRIBUCION DE BERNOULLI\n")
cat("═══════════════════════════════════════════════════════════════════\n\n")
# --- 2.1 Simulacion de una variable Bernoulli ---
cat("--- 2.1 Simulación de una variable Bernoulli ---\n\n")
# Fijamos la semilla para reproducibilidad
set.seed(2026)
# Simulamos 100 ensayos Bernoulli con p = 0.72
# (Ej: 72% de probabilidad de que una persona trabaje)
p <- 0.72
n <- 100
y_bernoulli <- rbinom(n, size = 1, prob = p)
# Veamos las primeras observaciones
cat("Primeras 20 observaciones (0 = no trabaja, 1 = trabaja):\n")
print(y_bernoulli[1:20])
# Tabla de frecuencias
cat("\nTabla de frecuencias:\n")
print(table(y_bernoulli))
# Proporcion muestral (debe estar cerca de p = 0.72)
cat("\nProporción muestral de y = 1:", mean(y_bernoulli), "\n")
cat("Valor teórico de p: ", p, "\n")
# --- 2.2 Gráfico de la distribución Bernoulli ---
cat("\n--- 2.2 Gráfico de la distribución Bernoulli ---\n\n")
par(mfrow = c(1, 3), mar = c(4, 4, 3, 1))
for (p_val in c(0.2, 0.5, 0.8)) {
bp <- barplot(c(1 - p_val, p_val),
names.arg = c("y = 0", "y = 1"),
col = c("#E74C3C", "#2ECC71"),
ylim = c(0, 1),
ylab = "Probabilidad",
main = paste0("Bernoulli(p = ", p_val, ")"),
border = "white")
abline(h = 0)
text(bp[1], 1 - p_val + 0.05, sprintf("%.1f", 1 - p_val), font = 2)
text(bp[2], p_val + 0.05, sprintf("%.1f", p_val), font = 2)
}
par(mfrow = c(1, 1))
# --- 2.3 Momentos de la Bernoulli ---
cat("--- 2.3 Momentos de la distribución Bernoulli ---\n\n")
cat("Para p = 0.72:\n")
cat(" Esperanza: E(y) = p =", p, "\n")
cat(" Varianza: Var(y) = p*(1-p) =", p * (1 - p), "\n")
cat(" Desv. típ.: SD(y) = sqrt(Var) =", sqrt(p * (1 - p)), "\n")
# Verificamos con los datos simulados
cat("\nVerificación con datos simulados (n =", n, "):\n")
cat(" Media muestral: ", round(mean(y_bernoulli), 4), "\n")
cat(" Varianza muestral: ", round(var(y_bernoulli), 4), "\n")
# NOTA: La varianza depende de p, que a su vez depende de x_i
# Esto VIOLA la homocedasticidad → MCO no funciona bien
pausa()
# ============================================================================
# SECCION 3: POR QUE MCO FALLA CON VARIABLES BINARIAS (MLP)
# ============================================================================
cat("\n═══════════════════════════════════════════════════════════════════\n")
cat(" SECCION 3: POR QUE MCO FALLA (MODELO LINEAL DE PROBABILIDAD)\n")
cat("============================================================\n")
cat("═══════════════════════════════════════════════════════════════════\n\n")
# --- 3.1 Generar datos con estructura Probit ---
cat("--- 3.1 Simulación del Modelo Lineal de Probabilidad ---\n\n")
set.seed(2026)
n <- 200
x <- runif(n, -3, 6)
# Probabilidad real de y=1 (modelo Probit)
p_real <- pnorm(0.5 + 0.8 * x)
# Generar variable binaria
y <- rbinom(n, 1, p_real)
# Estimar con MCO (Modelo Lineal de Probabilidad)
mlp <- lm(y ~ x)
cat("Resultados del Modelo Lineal de Probabilidad (MCO):\n\n")
print(summary(mlp))
# --- 3.2 Problemas del MLP ---
cat("\n--- 3.2 Problemas del MLP: predicciones fuera de [0, 1] ---\n\n")
predicciones <- predict(mlp)
cat("Rango de predicciones del MLP:\n")
cat(" Mínimo:", round(min(predicciones), 4), "\n")
cat(" Máximo:", round(max(predicciones), 4), "\n")
n_fuera <- sum(predicciones < 0 | predicciones > 1)
cat(" Predicciones fuera de [0,1]:", n_fuera, "de", n,
sprintf("(%.1f%%)\n", 100 * n_fuera / n))
# --- 3.3 Gráfico comparativo MLP vs Probit ---
cat("\n--- 3.3 Gráfico comparativo ---\n\n")
par(mar = c(4.5, 4.5, 3, 1))
plot(x, y, pch = 19, col = adjustcolor("gray40", alpha = 0.4),
xlab = "Variable explicativa (x)",
ylab = "y / Probabilidad predicha",
main = "MCO (MLP) vs Modelo correcto (Probit)",
ylim = c(-0.3, 1.3))
abline(mlp, col = "#E74C3C", lwd = 3)
abline(h = c(0, 1), lty = 2, col = "gray60")
x_ord <- sort(x)
lines(x_ord, pnorm(0.5 + 0.8 * x_ord), col = "#2ECC71", lwd = 3)
legend("bottomright",
legend = c("MCO (MLP)", "Probit", "Límites 0 y 1"),
col = c("#E74C3C", "#2ECC71", "gray60"),
lwd = c(3, 3, 1), lty = c(1, 1, 2), bty = "n")
# Zonas imposibles
rect(-4, -0.5, 8, 0,
col = adjustcolor("#E74C3C", alpha = 0.08), border = NA)
rect(-4, 1, 8, 1.5,
col = adjustcolor("#E74C3C", alpha = 0.08), border = NA)
# --- 3.4 Heterocedasticidad inherente ---
cat("--- 3.4 La heterocedasticidad es inherente al MLP ---\n\n")
cat("En el MLP, Var(u_i) = p_i * (1 - p_i)\n")
cat("Como p_i depende de x_i, la varianza NO es constante.\n\n")
# Mostramos cómo la varianza cambia
p_ejemplo <- seq(0.05, 0.95, by = 0.05)
var_ejemplo <- p_ejemplo * (1 - p_ejemplo)
par(mar = c(4.5, 4.5, 3, 1))
plot(p_ejemplo, var_ejemplo, type = "b", pch = 19,
col = "#3498DB", lwd = 2,
xlab = "p (probabilidad)", ylab = "Var(u) = p(1-p)",
main = "Varianza del error en el MLP (depende de p)")
abline(v = 0.5, lty = 2, col = "#E74C3C")
text(0.55, 0.24, "Máxima varianza\nen p = 0.5",
col = "#C0392B", cex = 0.8)
pausa()
# ============================================================================
# SECCION 4: DISTRIBUCION BINOMIAL
# ============================================================================
cat("\n═══════════════════════════════════════════════════════════════════\n")
cat(" SECCION 4: DISTRIBUCION BINOMIAL\n")
cat("═══════════════════════════════════════════════════════════════════\n\n")
# --- 4.1 La Binomial como suma de Bernoullis ---
cat("--- 4.1 La Binomial como suma de Bernoullis ---\n\n")
# Si n personas tienen cada una probabilidad p = 0.3 de comprar,
# ¿cuántas comprarán en total?
n_b <- 20
p_b <- 0.3
k <- 0:n_b
cat("Binomial(n =", n_b, ", p =", p_b, "):\n")
cat(" Esperanza: E(k) = n*p =", n_b * p_b, "\n")
cat(" Varianza: Var(k) = n*p*(1-p) =", n_b * p_b * (1 - p_b), "\n\n")
# Tabla de probabilidades (primeros valores)
prob_bin <- dbinom(k, n_b, p_b)
cat("Probabilidades P(k) para k = 0, 1, ..., 10:\n")
for (i in 0:10) {
cat(sprintf(" P(k = %2d) = %.4f\n", i, prob_bin[i + 1]))
}
# --- 4.2 Gráficos de la Binomial ---
cat("\n--- 4.2 Gráficos de la Binomial ---\n\n")
par(mfrow = c(1, 3), mar = c(4, 4, 3, 1))
# Caso 1: n=10, p=0.3
n1 <- 10; p1 <- 0.3
barplot(dbinom(0:n1, n1, p1), names.arg = 0:n1,
col = "#3498DB", border = "white",
ylab = "P(k)", main = paste0("Bin(n=", n1, ", p=", p1, ")"))
# Caso 2: n=20, p=0.5
n2 <- 20; p2 <- 0.5
barplot(dbinom(0:n2, n2, p2), names.arg = 0:n2,
col = "#9B59B6", border = "white",
ylab = "P(k)", main = paste0("Bin(n=", n2, ", p=", p2, ")"),
cex.names = 0.6)
# Caso 3: n=50, p=0.7
n3 <- 50; p3 <- 0.7
barplot(dbinom(0:n3, n3, p3), names.arg = 0:n3,
col = "#E67E22", border = "white",
ylab = "P(k)", main = paste0("Bin(n=", n3, ", p=", p3, ")"),
cex.names = 0.4)
par(mfrow = c(1, 1))
pausa()
# ============================================================================
# SECCION 5: VARIABLE LATENTE
# ============================================================================
cat("\n═══════════════════════════════════════════════════════════════════\n")
cat(" SECCION 5: EL CONCEPTO DE VARIABLE LATENTE\n")
cat("═══════════════════════════════════════════════════════════════════\n\n")
# --- 5.1 Simulación de la variable latente ---
cat("--- 5.1 Simulación: variable latente vs observada ---\n\n")
set.seed(42)
n <- 150
x_lat <- seq(-3, 3, length.out = n)
# Variable latente (utilidad neta de trabajar)
y_star <- 0.5 + 1.2 * x_lat + rnorm(n, 0, 1.5)
# Variable observada (decisión)
y_obs <- ifelse(y_star > 0, 1, 0)
cat("Resumen de la variable latente y*:\n")
print(summary(y_star))
cat("\nTabla de la variable observada y:\n")
print(table(y_obs))
cat("\nInterpretación:\n")
cat(" y* > 0 => la utilidad de trabajar > utilidad de no trabajar => y = 1\n")
cat(" y* ≤ 0 => la utilidad de trabajar ≤ utilidad de no trabajar => y = 0\n")
# --- 5.2 Gráfico doble panel ---
cat("\n--- 5.2 Gráfico: latente vs observada ---\n\n")
par(mfrow = c(1, 2), mar = c(4.5, 4.5, 3, 1))
# Panel izquierdo: variable latente
plot(x_lat, y_star, pch = 19,
col = ifelse(y_star > 0,
adjustcolor("#2ECC71", 0.6),
adjustcolor("#E74C3C", 0.6)),
xlab = "x (educación, experiencia...)",
ylab = "y* (utilidad latente)",
main = "Variable latente y*")
abline(h = 0, lty = 2, lwd = 2, col = "gray30")
abline(a = 0.5, b = 1.2, col = "#3498DB", lwd = 2)
legend("bottomright",
legend = c("y* > 0 (trabaja)", "y* ≤ 0 (no trabaja)"),
col = c("#2ECC71", "#E74C3C"), pch = 19, bty = "n", cex = 0.8)
# Panel derecho: variable observada
plot(x_lat, y_obs, pch = 19,
col = ifelse(y_obs == 1,
adjustcolor("#2ECC71", 0.5),
adjustcolor("#E74C3C", 0.5)),
xlab = "x", ylab = "y (observada)",
main = "Variable observada y", yaxt = "n", ylim = c(-0.1, 1.1))
axis(2, at = c(0, 1))
lines(sort(x_lat), pnorm(0.5 + 1.2 * sort(x_lat), 0, 1.5),
col = "#3498DB", lwd = 3)
legend("right", legend = "P(y=1|x)", col = "#3498DB",
lwd = 3, bty = "n")
par(mfrow = c(1, 1))
pausa()
# ============================================================================
# SECCION 6: VARIABLES MULTINOMIALES
# ============================================================================
cat("\n═══════════════════════════════════════════════════════════════════\n")
cat(" SECCION 6: VARIABLES MULTINOMIALES\n")
cat("═══════════════════════════════════════════════════════════════════\n\n")
# --- 6.1 Simulación de una variable multinomial no ordenada ---
cat("--- 6.1 Elección de medio de transporte (no ordenada) ---\n\n")
set.seed(2026)
probs_transp <- c(Coche = 0.45, Autobus = 0.25, Tren = 0.20, Avion = 0.10)
cat("Probabilidades teóricas:\n")
print(probs_transp)
# Simulamos 1000 elecciones
elecciones <- sample(names(probs_transp), 1000,
replace = TRUE, prob = probs_transp)
cat("\nFrecuencias observadas (n = 1000):\n")
freq <- table(factor(elecciones, levels = names(probs_transp)))
print(freq)
cat("\nProporciones observadas:\n")
print(round(prop.table(freq), 3))
# Gráfico de barras
par(mar = c(4.5, 4.5, 3, 1))
bp <- barplot(prop.table(freq),
col = c("#3498DB", "#2ECC71", "#E67E22", "#9B59B6"),
border = "white", ylim = c(0, 0.55),
ylab = "Proporción",
main = "Elección de medio de transporte (n = 1000)")
text(bp, as.numeric(prop.table(freq)) + 0.02,
sprintf("%.1f%%", as.numeric(prop.table(freq)) * 100), font = 2)
# --- 6.2 Variable multinomial ordenada ---
cat("\n--- 6.2 Nivel de satisfacción (ordenada) ---\n\n")
set.seed(2026)
# Simulamos satisfacción: 0=Muy insatisfecho, ..., 4=Muy satisfecho
probs_satis <- c(0.05, 0.10, 0.25, 0.35, 0.25)
nombres_satis <- c("Muy insatisf.", "Insatisfecho",
"Neutro", "Satisfecho", "Muy satisf.")
satisfaccion <- sample(0:4, 500, replace = TRUE, prob = probs_satis)
cat("Tabla de frecuencias:\n")
freq_s <- table(factor(satisfaccion, levels = 0:4, labels = nombres_satis))
print(freq_s)
# Gráfico
par(mar = c(6, 4.5, 3, 1))
bp2 <- barplot(prop.table(freq_s),
col = c("#E74C3C", "#E67E22", "#F1C40F", "#2ECC71", "#27AE60"),
border = "white", ylim = c(0, 0.45),
ylab = "Proporción",
main = "Nivel de satisfacción (variable ordenada)",
las = 2, cex.names = 0.8)
text(bp2, as.numeric(prop.table(freq_s)) + 0.02,
sprintf("%.0f%%", as.numeric(prop.table(freq_s)) * 100), font = 2)
# --- 6.3 Distribución Multinomial en R ---
cat("\n--- 6.3 La distribución Multinomial en R ---\n\n")
# La función dmultinom calcula P(k1, k2, ..., kJ)
# Ejemplo: P(exactamente 45 coches, 25 autobuses, 20 trenes, 10 aviones)
# en una muestra de 100
prob_exacta <- dmultinom(c(45, 25, 20, 10), prob = probs_transp)
cat("P(k = [45, 25, 20, 10] | n=100, p=[0.45, 0.25, 0.20, 0.10]):\n")
cat(" ", formatC(prob_exacta, format = "e", digits = 4), "\n")
# Simulación de la multinomial
cat("\nSimulación de 5 muestras multinomiales (n = 100):\n")
for (i in 1:5) {
muestra <- rmultinom(1, 100, probs_transp)
cat(sprintf(" Muestra %d: Coche=%d, Autobús=%d, Tren=%d, Avión=%d\n",
i, muestra[1], muestra[2], muestra[3], muestra[4]))
}
pausa()
# ============================================================================
# SECCION 7: DATOS CENSURADOS
# ============================================================================
cat("\n═══════════════════════════════════════════════════════════════════\n")
cat(" SECCION 7: DATOS CENSURADOS\n")
cat("═══════════════════════════════════════════════════════════════════\n\n")
# --- 7.1 Simulación de datos censurados ---
cat("--- 7.1 Simulación: gasto en seguros médicos (censurado en 0) ---\n\n")
set.seed(2026)
n <- 300
renta <- runif(n, 15, 80) # renta en miles de euros
# Variable latente: gasto potencial
gasto_latente <- -20 + 0.6 * renta + rnorm(n, 0, 8)
# Variable observada: censurada en 0
gasto_observado <- pmax(gasto_latente, 0)
cat("Resumen de la variable latente (gasto potencial):\n")
print(round(summary(gasto_latente), 2))
cat("\nResumen de la variable observada (censurada en 0):\n")
print(round(summary(gasto_observado), 2))
n_censurados <- sum(gasto_observado == 0)
cat("\nObservaciones censuradas (gasto = 0):", n_censurados,
sprintf("(%.1f%%)\n", 100 * n_censurados / n))
# --- 7.2 MCO sesgado con datos censurados ---
cat("\n--- 7.2 MCO está sesgado con datos censurados ---\n\n")
# MCO con datos censurados (SESGADO)
mco_censurado <- lm(gasto_observado ~ renta)
cat("MCO con datos censurados:\n")
cat(" Intercepto:", round(coef(mco_censurado)[1], 4), "\n")
cat(" Pendiente: ", round(coef(mco_censurado)[2], 4), "\n")
cat("\nValores verdaderos:\n")
cat(" Intercepto: -20\n")
cat(" Pendiente: 0.6\n")
cat("\n¡MCO subestima la pendiente y sobreestima el intercepto!\n")
cat("El modelo Tobit (Tema 3) corrige este sesgo.\n")
# --- 7.3 Gráfico de censura ---
cat("\n--- 7.3 Gráfico ---\n\n")
par(mfrow = c(1, 2), mar = c(4.5, 4.5, 3, 1))
# Panel 1: Variable latente
plot(renta, gasto_latente, pch = 19,
col = ifelse(gasto_latente > 0,
adjustcolor("#2ECC71", 0.5),
adjustcolor("#E74C3C", 0.5)),
xlab = "Renta (miles €)", ylab = "Gasto latente (miles €)",
main = "Variable latente y*")
abline(h = 0, lty = 2, lwd = 2)
abline(a = -20, b = 0.6, col = "#3498DB", lwd = 2)
# Panel 2: Variable censurada
plot(renta, gasto_observado, pch = 19,
col = ifelse(gasto_observado > 0,
adjustcolor("#2ECC71", 0.5),
adjustcolor("#E74C3C", 0.5)),
xlab = "Renta (miles €)", ylab = "Gasto observado (miles €)",
main = "Variable censurada y")
abline(h = 0, lty = 2, lwd = 2)
abline(mco_censurado, col = "#E74C3C", lwd = 2)
abline(a = -20, b = 0.6, col = "#2ECC71", lwd = 2, lty = 2)
legend("topleft",
legend = c("MCO (sesgado)", "Verdadera"),
col = c("#E74C3C", "#2ECC71"), lwd = 2, lty = c(1, 2),
bty = "n", cex = 0.8)
par(mfrow = c(1, 1))
# --- 7.4 Histograma de la variable censurada ---
cat("--- 7.4 Histograma de la variable censurada ---\n\n")
par(mar = c(4.5, 4.5, 3, 1))
hist(gasto_observado, breaks = 30,
col = adjustcolor("#3498DB", 0.6),
border = "white",
main = "Distribución de la variable censurada",
xlab = "Gasto observado (miles €)", ylab = "Frecuencia")
cat("La distribución tiene dos componentes:\n")
cat(" 1. Componente DISCRETO: masa de probabilidad en y = 0\n")
cat(" 2. Componente CONTINUO: para los valores y > 0\n")
cat("Esta mezcla es la razón por la que MCO no funciona.\n")
# --- 7.5 Tipos de censura ---
cat("\n--- 7.5 Tipos de censura (izquierda y derecha) ---\n\n")
set.seed(123)
y_star_c <- rnorm(500, 5, 3)
par(mfrow = c(1, 2), mar = c(4.5, 4.5, 3, 1))
# Censura por la izquierda
y_cens_izq <- pmax(y_star_c, 0)
hist(y_cens_izq, breaks = 30,
col = adjustcolor("#3498DB", 0.6), border = "white",
main = "Censura izquierda (c = 0)",
xlab = "y", ylab = "Frecuencia")
abline(v = 0, col = "#E74C3C", lwd = 2, lty = 2)
# Censura por la derecha
y_cens_der <- pmin(y_star_c, 10)
hist(y_cens_der, breaks = 30,
col = adjustcolor("#E67E22", 0.6), border = "white",
main = "Censura derecha (c = 10)",
xlab = "y", ylab = "Frecuencia")
abline(v = 10, col = "#E74C3C", lwd = 2, lty = 2)
par(mfrow = c(1, 1))
cat("Censura izquierda: y = max(y*, c) → no observamos valores < c\n")
cat("Censura derecha: y = min(y*, c) → no observamos valores > c\n")
pausa()
# ============================================================================
# SECCION 8: DATOS TRUNCADOS
# ============================================================================
cat("\n═══════════════════════════════════════════════════════════════════\n")
cat(" SECCION 8: DATOS TRUNCADOS\n")
cat("═══════════════════════════════════════════════════════════════════\n\n")
# --- 8.1 Simulación de datos truncados ---
cat("--- 8.1 Truncamiento vs. censura ---\n\n")
set.seed(2026)
n <- 300
x_tr <- runif(n, 0, 10)
y_star_tr <- -2 + 1.5 * x_tr + rnorm(n, 0, 2)
# Censura: observamos todos los x, pero y se fija en 0
y_censurada <- pmax(y_star_tr, 0)
n_cens <- sum(y_star_tr <= 0)
# Truncamiento: PERDEMOS las observaciones con y* <= 0
seleccion <- y_star_tr > 0
x_truncado <- x_tr[seleccion]
y_truncado <- y_star_tr[seleccion]
n_trunc_perdidos <- sum(!seleccion)
cat("Datos completos: n =", n, "\n")
cat("Censura: ", n_cens, "obs con y=0,",
n - n_cens, "obs con y>0 (total:", n, ")\n")
cat("Truncamiento: ", n_trunc_perdidos, "obs PERDIDAS,",
sum(seleccion), "obs restantes (total:", sum(seleccion), ")\n")
cat("\n¡En el truncamiento PERDEMOS datos e información!\n")
cat("En la censura al menos conservamos las X de todos.\n")
# --- 8.2 Gráfico comparativo ---
cat("\n--- 8.2 Gráfico comparativo ---\n\n")
par(mfrow = c(1, 2), mar = c(4.5, 4.5, 3, 1))
# Censura
plot(x_tr, y_censurada, pch = 19,
col = ifelse(y_star_tr > 0,
adjustcolor("#2ECC71", 0.5),
adjustcolor("#E74C3C", 0.5)),
xlab = "x", ylab = "y", main = "Censura en 0")
abline(h = 0, lty = 2, lwd = 2)
legend("topleft", legend = c("y > 0", "y = 0 (censurados)"),
col = c("#2ECC71", "#E74C3C"), pch = 19, bty = "n", cex = 0.8)
# Truncamiento
plot(x_truncado, y_truncado, pch = 19,
col = adjustcolor("#2ECC71", 0.5),
xlab = "x", ylab = "y", main = "Truncamiento en 0",
xlim = range(x_tr), ylim = range(y_star_tr))
abline(h = 0, lty = 2, lwd = 2)
rect(par("usr")[1], par("usr")[3], par("usr")[2], 0,
col = adjustcolor("#E74C3C", 0.08), border = NA)
text(8, -3, "Datos perdidos", col = "#C0392B", font = 3, cex = 0.85)
par(mfrow = c(1, 1))
pausa()
# ============================================================================
# SECCION 9: DISTRIBUCIONES CENSURADA Y TRUNCADA
# ============================================================================
cat("\n═══════════════════════════════════════════════════════════════════\n")
cat(" SECCION 9: DISTRIBUCIONES CENSURADA Y TRUNCADA\n")
cat("═══════════════════════════════════════════════════════════════════\n\n")
# --- 9.1 Normal censurada ---
cat("--- 9.1 Distribución Normal censurada ---\n\n")
mu <- 2; sigma <- 3
y_seq <- seq(-8, 12, length.out = 500)
dens <- dnorm(y_seq, mu, sigma)
# Masa en 0
p0 <- pnorm(0, mu, sigma)
cat("Para y* ~ N(mu=2, sigma=3), censurada en 0:\n")
cat(" P(y = 0) = Phi(-mu/sigma) = Phi(", round(-mu/sigma, 4), ") =",
round(p0, 4), "\n")
cat(" Es decir, el", sprintf("%.1f%%", p0 * 100),
"de las observaciones quedan censuradas en 0.\n")
par(mar = c(4.5, 4.5, 3, 1))
plot(y_seq[y_seq >= 0], dens[y_seq >= 0], type = "l",
lwd = 3, col = "#3498DB",
xlab = "y", ylab = "Densidad / Probabilidad",
main = "Normal censurada: mu=2, sigma=3",
xlim = c(-2, 12), ylim = c(0, max(dens) * 1.3))
polygon(c(0, y_seq[y_seq >= 0], max(y_seq)),
c(0, dens[y_seq >= 0], 0),
col = adjustcolor("#3498DB", 0.2), border = NA)
segments(0, 0, 0, p0, col = "#E74C3C", lwd = 4)
points(0, p0, pch = 19, col = "#E74C3C", cex = 2)
text(0.5, p0 + 0.01, sprintf("P(y=0) = %.3f", p0),
col = "#C0392B", font = 2, cex = 0.8, pos = 4)
lines(y_seq[y_seq < 0], dens[y_seq < 0],
lwd = 2, lty = 3, col = "gray60")
# --- 9.2 Normal truncada ---
cat("\n--- 9.2 Distribución Normal truncada ---\n\n")
y_pos <- seq(0.001, 12, length.out = 500)
dens_norm <- dnorm(y_pos, mu, sigma)
dens_trunc <- dens_norm / (1 - pnorm(0, mu, sigma))
cat("La densidad truncada se reescala dividiendo por P(y > 0):\n")
cat(" P(y > 0) = 1 - Phi(-mu/sigma) =", round(1 - p0, 4), "\n")
cat(" Factor de reescalado: 1 /", round(1 - p0, 4), "=",
round(1 / (1 - p0), 4), "\n")
par(mar = c(4.5, 4.5, 3, 1))
plot(y_pos, dens_trunc, type = "l", lwd = 3, col = "#E67E22",
xlab = "y", ylab = "Densidad",
main = "Normal truncada vs. normal original",
ylim = c(0, max(dens_trunc) * 1.1))
lines(y_pos, dens_norm, lwd = 2, lty = 2, col = "#3498DB")
polygon(c(y_pos, rev(y_pos)),
c(dens_trunc, rep(0, length(y_pos))),
col = adjustcolor("#E67E22", 0.15), border = NA)
legend("topright",
legend = c("Normal truncada", "Normal original"),
col = c("#E67E22", "#3498DB"), lwd = c(3, 2), lty = c(1, 2),
bty = "n")
abline(v = 0, lty = 3, col = "gray50")
pausa()
# ============================================================================
# SECCION 10: RESUMEN
# ============================================================================
cat("\n═══════════════════════════════════════════════════════════════════\n")
cat(" SECCION 10: RESUMEN\n")
cat("═══════════════════════════════════════════════════════════════════\n\n")
cat("Tipos de variables dependientes limitadas:\n\n")
cat(" 1. BINARIAS (y ∈ {0,1}):\n")
cat(" - Distribución: Bernoulli\n")
cat(" - Modelos: Probit, Logit (Tema 2)\n")
cat(" - Concepto clave: variable latente\n\n")
cat(" 2. MULTINOMIALES (y ∈ {0,1,...,J}):\n")
cat(" - Distribución: Multinomial\n")
cat(" - Tipos: ordenadas y no ordenadas\n")
cat(" - Modelos: Logit/Probit multinomial (avanzado)\n\n")
cat(" 3. CONTINUAS LIMITADAS:\n")
cat(" - CENSURADAS: observamos X de todos, pero y fijada en un punto\n")
cat(" Modelo: Tobit (Tema 3)\n")
cat(" - TRUNCADAS: perdemos TODA la información fuera del rango\n")
cat(" Modelo: Heckman (selección muestral)\n\n")
cat("PUNTO CLAVE: MCO no funciona con ninguno de estos tipos de datos.\n")
cat(" Necesitamos modelos específicos para cada caso.\n")
cat("\n╔══════════════════════════════════════════════════════════════════╗\n")
cat("║ FIN del Script — Tema 1 ║\n")
cat("║ Siguiente: Tema 2 — Modelos Probit y Logit ║\n")
cat("╚══════════════════════════════════════════════════════════════════╝\n")