Del modelo económico al código Python: endogeneidad, aleatorización y variables instrumentales
Un error frecuente en el análisis de datos es ir directamente a los datos sin haber pensado en el modelo que los genera. Python puede calcular cualquier correlación en milisegundos — pero no puede decirte si esa correlación tiene interpretación causal. Eso requiere teoría.
En este bloque formalizamos el problema antes de escribir una sola línea de código. El objetivo es que cuando llegues al Bloque II, cada línea de Python te resulte una traducción directa de algo que ya entiendes.
El Servicio Público de Empleo ofrece un curso de formación profesional de tres meses a desempleados. Queremos saber si el curso aumenta la probabilidad de encontrar empleo. Disponemos de datos observacionales: sabemos quién asistió y quién encontró trabajo, pero la participación fue voluntaria.
El modelo estructural que queremos estimar es:
$Y_i = 1$ si el desempleado $i$ encuentra empleo en los 6 meses siguientes · $D_i = 1$ si asiste al curso · $U_i$ = todo lo que afecta al reempleo y el analista no observa
La participación voluntaria genera selección positiva: se apuntan los desempleados con mayor motivación y capital humano latente. Dos variables en $U_i$:
OLS requiere $\Cov(D_i, U_i) = 0$. Aquí se viola porque los más motivados (mayor $U_i$) tienen más probabilidad de apuntarse ($D_i = 1$), por lo que $\Cov(D_i, U_i) > 0$. Aplicando la fórmula del sesgo:
OLS sobreestima el efecto causal del curso. Parte del mayor reempleo observado entre los participantes se debe a su mayor motivación de base, no al curso.
La diferencia de medias que calcula el gestor es:
El segundo término es positivo: los participantes ya eran más empleables antes del curso. El gestor atribuye al programa algo que habría ocurrido de todas formas.
Para resolver el problema de endogeneidad, el servicio de empleo decide sortear las plazas del curso entre los solicitantes. Solo entran quienes salen elegidos en el sorteo. Llamamos $Z_i = 1$ si el desempleado $i$ es sorteado, $Z_i = 0$ si no.
En presencia de incumplimiento parcial — algunos sorteados no asisten, algunos no sorteados consiguen plaza por otras vías — el modelo completo tiene tres ecuaciones:
| Perfil | Si es sorteado ($Z=1$) | Si no es sorteado ($Z=0$) |
|---|---|---|
| Asiste solo si le toca el sorteo | $D=1$ | $D=0$ |
| Siempre busca plaza, con o sin sorteo | $D=1$ | $D=1$ |
| No asiste nunca, aunque le toque | $D=0$ | $D=0$ |
El sorteo asigna $Z_i$ de forma completamente aleatoria, sin usar ninguna información sobre el desempleado. Por tanto $Z_i$ es estadísticamente independiente de cualquier característica — observada o no — del individuo:
Si las plazas se asignaran por orden de solicitud, los más motivados se apuntarían antes y tendrían más probabilidad de obtener plaza. La asignación volvería a depender de $U_i$ y la exogeneidad se rompería: los grupos $Z=1$ y $Z=0$ ya no serían comparables.
El ITT mide el efecto de ser sorteado sobre el reempleo — el efecto de la asignación, no del tratamiento efectivo. Con incumplimiento, algunos sorteados no asisten, por lo que el ITT diluye el efecto real:
El estimador de Wald corrige esa dilución dividiendo el ITT por la primera etapa $\pi_1$ — la fracción de individuos que cambia su comportamiento al cambiar $Z_i$. De ese modo recupera el efecto causal para quienes responden al instrumento:
| Perfil | Tipo |
|---|---|
| Asiste solo si le toca el sorteo | complier |
| Siempre busca plaza, con o sin sorteo | always-taker |
| No asiste nunca, aunque le toque | never-taker |
El estimador de Wald identifica el LATE (Local Average Treatment Effect): el efecto del curso exclusivamente para los compliers. Son los únicos cuyo tratamiento varía con $Z_i$, y por tanto los únicos sobre los que el instrumento tiene información.
En el Bloque I formalizamos el modelo. Ahora lo traducimos a Python. El código que sigue simula exactamente el DGP que describiste en las ecuaciones — cada línea corresponde a una pieza del modelo teórico.
Lee el código completo antes de responder las preguntas. No hace falta que lo ejecutes todavía: el objetivo es que seas capaz de localizar en el código cada elemento del modelo y entender por qué está escrito como está.
El código está organizado en cuatro bloques funcionales. Lee cada uno con atención y localiza los elementos del modelo teórico del Bloque I antes de pasar a las preguntas.
1 · Parámetros del modelo
import numpy as np import pandas as pd import matplotlib.pyplot as plt import statsmodels.formula.api as smf from linearmodels.iv import IV2SLS # ── Todos los parámetros del DGP en un único diccionario ────────────────────── params = { "n" : 3000, # tamaño muestral "seed" : 123, # semilla para reproducibilidad "alpha" : 0.0, # intercepto del modelo de resultado "beta" : 0.10, # efecto causal verdadero del curso sobre el reempleo "gamma" : 1.0, # fuerza de la endogeneidad: cuánto influye U en D "delta_Z" : 1.2, # relevancia del instrumento: cuánto influye Z en D "delta_U" : 1.5, # efecto de U sobre Y (lo que OLS no puede ver) }
2 · Función que genera los datos
def simular_dgp(p): """ Genera los microdatos del experimento de formación. Recibe el diccionario de parámetros p y devuelve un DataFrame. Modelo: U_i ~ N(0,1) motivación latente (inobservable) Z_i ~ Bernoulli(0.5) resultado del sorteo (instrumento exógeno) D_i = 1 si γ·U + δ_Z·Z + ε > 0 (asistencia efectiva al curso) Y_i = 1 si α + β·D + δ_U·U + η > 0 (reempleo en 6 meses) """ np.random.seed(p["seed"]) # ── Variables latentes ─────────────────────────────────────────────────────── U = np.random.normal(0, 1, p["n"]) # inobservable: motivación Z = np.random.binomial(1, 0.5, p["n"]) # instrumento: sorteo aleatorio eps = np.random.normal(0, 1, p["n"]) # ruido en la elección de D eta = np.random.normal(0, 1, p["n"]) # ruido en el resultado Y # ── Tratamiento efectivo D (primera etapa) ─────────────────────────────────── # D=1 si el índice latente supera el umbral 0 D = (p["gamma"] * U + p["delta_Z"] * Z + eps > 0).astype(int) # ── Resultado Y (modelo estructural) ───────────────────────────────────────── Y_star = p["alpha"] + p["beta"] * D + p["delta_U"] * U + eta Y = (Y_star > 0).astype(int) # ── Clasificación de individuos por tipo ───────────────────────────────────── # Para clasificar necesitamos el tratamiento potencial bajo Z=0 y Z=1 D_z0 = (p["gamma"] * U + p["delta_Z"] * 0 + eps > 0).astype(int) D_z1 = (p["gamma"] * U + p["delta_Z"] * 1 + eps > 0).astype(int) tipos = [] for d0, d1 in zip(D_z0, D_z1): if d0 == 0 and d1 == 1: tipos.append("complier") elif d0 == 1 and d1 == 1: tipos.append("always_taker") elif d0 == 0 and d1 == 0: tipos.append("never_taker") else: tipos.append("defier") # raro: D=1→0 return pd.DataFrame({ "reempleo" : Y, "curso" : D, "sorteo" : Z, "U" : U, "tipo" : tipos, })
3 · Función que calcula los estimadores
def estimar_modelos(df, beta_verdadero): """ Calcula OLS, ITT y el estimador de Wald sobre el DataFrame df. Devuelve un diccionario con los resultados de cada estimador. """ resultados = {} # ── OLS: ignora la endogeneidad ────────────────────────────────────────────── ols = smf.ols("reempleo ~ curso", data=df).fit() resultados["OLS"] = { "estimacion" : round(ols.params["curso"], 4), "p_valor" : round(ols.pvalues["curso"], 4), "sesgo" : round(ols.params["curso"] - beta_verdadero, 4), } # ── ITT: efecto de la asignación (no del tratamiento efectivo) ─────────────── medias_Z = df.groupby("sorteo")["reempleo"].mean() itt = medias_Z[1] - medias_Z[0] resultados["ITT"] = { "estimacion" : round(itt, 4), "p_valor" : None, "sesgo" : round(itt - beta_verdadero, 4), } # ── Primera etapa: cuánto mueve el sorteo la asistencia ────────────────────── medias_D = df.groupby("sorteo")["curso"].mean() pi1 = medias_D[1] - medias_D[0] # ── Wald: ITT / primera etapa ──────────────────────────────────────────────── wald = itt / pi1 resultados["Wald"] = { "estimacion" : round(wald, 4), "p_valor" : None, "sesgo" : round(wald - beta_verdadero, 4), } # ── IV / 2SLS: verificación algebraica de Wald ─────────────────────────────── iv = IV2SLS.from_formula( "reempleo ~ 1 + [curso ~ sorteo]", data=df ).fit(cov_type="robust") resultados["IV"] = { "estimacion" : round(float(iv.params["curso"]), 4), "p_valor" : round(float(iv.pvalues["curso"]), 4), "sesgo" : round(float(iv.params["curso"]) - beta_verdadero, 4), } return resultados, round(pi1, 4)
4 · Función que genera los gráficos
def graficar_resultados(df, resultados, beta_verdadero): """ Genera dos gráficos: (1) Distribución de U por grupo Z: verifica la exogeneidad del instrumento (2) Comparación de estimadores vs. valor verdadero β """ fig, axes = plt.subplots(1, 2, figsize=(12, 4)) fig.patch.set_facecolor("#0f1117") # ── Gráfico 1: distribución de U por grupo Z ───────────────────────────────── ax1 = axes[0] ax1.set_facecolor("#0f1117") colores = {0: "#6b8cff", 1: "#4ecdc4"} etiquetas = {0: "No sorteado (Z=0)", 1: "Sorteado (Z=1)"} for grupo in [0, 1]: datos_grupo = df[df["sorteo"] == grupo]["U"] datos_grupo.plot.kde( ax=ax1, label=etiquetas[grupo], color=colores[grupo], linewidth=2 ) ax1.axvline(0, color="white", linestyle="--", alpha=0.2, linewidth=1) ax1.set_title("Distribución de U por grupo Z\n(debe ser idéntica si Z es exógeno)", color="#dde2f0", fontsize=10) ax1.tick_params(colors="#7a82a0") ax1.legend(fontsize=9, labelcolor="#dde2f0", facecolor="#1e2235", edgecolor="#2a2f45") for spine in ax1.spines.values(): spine.set_edgecolor("#2a2f45") # ── Gráfico 2: comparación de estimadores ──────────────────────────────────── ax2 = axes[1] ax2.set_facecolor("#0f1117") nombres = list(resultados.keys()) estimaciones = [resultados[k]["estimacion"] for k in nombres] colores_est = [] for k in nombres: sesgo = abs(resultados[k]["sesgo"]) if sesgo < 0.02: colores_est.append("#4ecdc4") # verde: cerca del verdadero elif sesgo < 0.08: colores_est.append("#ffd166") # amarillo: sesgo moderado else: colores_est.append("#ff7a6b") # rojo: sesgo grande barras = ax2.barh(nombres, estimaciones, color=colores_est, height=0.5, edgecolor="#2a2f45") ax2.axvline(beta_verdadero, color="white", linestyle="--", linewidth=1.5, label=f"β verdadero = {beta_verdadero}") for barra, val in zip(barras, estimaciones): ax2.text(val + 0.005, barra.get_y() + barra.get_height() / 2, f"{val:.3f}", va="center", color="#dde2f0", fontsize=9) ax2.set_title("Estimadores vs. valor verdadero β", color="#dde2f0", fontsize=10) ax2.tick_params(colors="#7a82a0") ax2.legend(fontsize=9, labelcolor="#dde2f0", facecolor="#1e2235", edgecolor="#2a2f45") for spine in ax2.spines.values(): spine.set_edgecolor("#2a2f45") plt.tight_layout() plt.savefig("resultados_formacion.png", dpi=150, facecolor="#0f1117", bbox_inches="tight") plt.show()
5 · Ejecución del análisis completo
# ── Generar datos ──────────────────────────────────────────────────────────── df = simular_dgp(params) # ── Resumen de tipos de individuos ─────────────────────────────────────────── conteo_tipos = df["tipo"].value_counts() for tipo, n in conteo_tipos.items(): pct = n / params["n"] * 100 print(f"{tipo:15s}: {n:5d} ({pct:.1f}%)") # ── Estimar y mostrar resultados ───────────────────────────────────────────── resultados, pi1 = estimar_modelos(df, params["beta"]) print(f"\nPrimera etapa π₁ = {pi1}") print(f"{'Estimador':10} {'Estimación':>12} {'Sesgo':>10}") print("-" * 36) for nombre, res in resultados.items(): print(f"{nombre:10} {res['estimacion']:>12.4f} {res['sesgo']:>10.4f}") # ── Gráficos ───────────────────────────────────────────────────────────────── graficar_resultados(df, resultados, params["beta"])
Responde las siguientes preguntas sin ejecutar el código todavía. Localiza en el código la línea o bloque que justifica cada respuesta.
params centraliza todos los
parámetros del DGP. ¿Qué ventaja tiene este diseño frente
a dispersar los valores numéricos por el código?
¿Qué parámetro controla la fuerza de la endogeneidad?
¿Y la relevancia del instrumento?
simular_dgp() la línea que implementa
la condición $\Cov(Z_i, U_i) = 0$.
¿Cómo garantiza el código que el sorteo es exógeno?
for d0, d1 in zip(D_z0, D_z1) clasifica
a cada individuo en su tipo. ¿Qué condición del condicional
corresponde a un complier? ¿Y a un never-taker?
¿Por qué se calculan D_z0 y D_z1
antes del bucle?
graficar_resultados() hay un condicional
que asigna un color a cada barra según el sesgo.
¿Qué color recibirá la barra de OLS con los parámetros
actuales? ¿Y la de Wald? Razona tu respuesta
sin ejecutar el código.
for nombre, res in resultados.items()
imprime la tabla de resultados. ¿Qué estructura de datos
es resultados? ¿Qué contiene cada valor
asociado a una clave como "OLS"?
¿Cómo accederías al p-valor del estimador IV?
Centralizar los parámetros en un diccionario tiene dos ventajas clave:
params; no hay
que buscar el valor disperso por el código. Esto es
especialmente útil en el Bloque III.
params["gamma"] sabes
exactamente qué parámetro se está usando.
params["gamma"]
— controla cuánto influye $U_i$ en $D_i$. params["delta_Z"]
— controla cuánto influye $Z_i$ en $D_i$.
La línea que garantiza la exogeneidad del instrumento es:
Z = np.random.binomial(1, 0.5, p["n"])
Z se genera de forma completamente independiente
de U: son dos llamadas separadas al generador
de números aleatorios que no comparten ningún argumento.
No hay ningún término de U en la construcción
de Z, lo que implementa directamente
$Z_i \perp U_i$.
Las condiciones del condicional son:
d0 == 0 and d1 == 1 —
no asiste si $Z=0$, asiste si $Z=1$
d0 == 0 and d1 == 0 —
no asiste en ningún caso
D_z0 y D_z1 se calculan antes
del bucle porque son arrays — se generan de una vez para
todos los individuos usando operaciones vectorizadas de numpy,
que son mucho más eficientes que calcularlos dentro del bucle
uno a uno.
Con los parámetros actuales (beta=0.10,
gamma=1.0):
resultados es un
diccionario de diccionarios:
las claves son los nombres de los estimadores
("OLS", "ITT",
"Wald", "IV") y cada valor
es a su vez un diccionario con tres claves:
"estimacion", "p_valor"
y "sesgo".
Para acceder al p-valor del estimador IV:
resultados["IV"]["p_valor"]
El código del Bloque II está diseñado para ser modificado. Todos los
parámetros del DGP están en el diccionario params —
basta con cambiar un valor, volver a ejecutar las tres funciones
y observar qué ocurre con los estimadores.
En cada ejercicio te pedimos que primero predijas el resultado usando la teoría del Bloque I, y luego lo verifiques ejecutando el código. La predicción es la parte más importante: si aciertas, es que has entendido el modelo.
Ponemos gamma = 0: la motivación latente $U_i$
ya no influye en la decisión de asistir al curso. La participación
sigue siendo voluntaria en apariencia, pero en realidad ya no
está correlacionada con $U_i$.
gamma = 0,
¿qué le ocurre a $\Cov(D_i, U_i)$? ¿Sigue siendo OLS sesgado?
¿Qué diferencia esperas ver entre OLS y Wald?
params y ejecuta el análisis completo:
params_exp1 = params.copy() params_exp1["gamma"] = 0.0 df1 = simular_dgp(params_exp1) res1, pi1_1 = estimar_modelos(df1, params_exp1["beta"]) graficar_resultados(df1, res1, params_exp1["beta"])¿Coincide el resultado con tu predicción? ¿Qué le ocurre al gráfico de densidades de $U$ por grupo $Z$?
gamma = 0, ¿tiene sentido usar IV en lugar de OLS?
¿Qué coste tiene usar IV cuando OLS ya sería válido?
Con gamma = 0, la ecuación de selección es:
D = (0 * U + delta_Z * Z + eps > 0)
$U_i$ desaparece de la decisión de asistir. Por tanto $\Cov(D_i, U_i) = 0$ y OLS recupera la consistencia. Esperamos que OLS y Wald produzcan estimaciones cercanas entre sí y ambas próximas a $\beta = 0.10$.
Al ejecutar, las barras de OLS y Wald quedan ambas
en verde — sesgo cercano a cero en los dos casos.
El gráfico de densidades muestra las distribuciones
de $U$ en $Z=0$ y $Z=1$ perfectamente superpuestas,
lo que confirma que $Z \perp U$ — aunque con
gamma=0 esto ya se cumplía también para $D$.
No tiene sentido usar IV si OLS es consistente. El coste de usar IV innecesariamente es una pérdida de eficiencia: IV tiene mayor varianza que OLS porque solo usa la variación en $D_i$ inducida por $Z_i$, que es una fracción de la variación total. Los errores estándar de IV son mayores y los intervalos de confianza más anchos.
Reducimos delta_Z de 1.2 a 0.1: el sorteo
apenas mueve la asistencia al curso. El instrumento sigue
siendo exógeno — el sorteo es aleatorio — pero casi nadie
cambia de comportamiento al ser sorteado.
delta_Z
es muy pequeño? ¿Cómo afecta eso al estimador de Wald?
Razona usando la fórmula
$\hat{\beta}^{Wald} = ITT / \hat{\pi}_1$.
params_exp2 = params.copy() params_exp2["delta_Z"] = 0.1 df2 = simular_dgp(params_exp2) res2, pi1_2 = estimar_modelos(df2, params_exp2["beta"]) graficar_resultados(df2, res2, params_exp2["beta"]) print(f"Primera etapa π₁ = {pi1_2}")¿Cuánto vale $\hat{\pi}_1$? ¿Es el estadístico F de la primera etapa mayor o menor que 10?
Con delta_Z muy pequeño, el sorteo apenas
cambia la probabilidad de asistir: $\hat{\pi}_1 \approx 0$.
El estimador de Wald divide el ITT (pequeño pero distinto
de cero) entre $\hat{\pi}_1$ (casi cero):
Dividir por un número cercano a cero amplifica enormemente cualquier variación muestral en el numerador. El estimador se vuelve muy inestable.
Con delta_Z = 0.1, la primera etapa
será muy débil: $\hat{\pi}_1 \approx 0.04$–$0.06$
y el estadístico F caerá muy por debajo de 10
(posiblemente por debajo de 2).
Los errores estándar de IV serán enormes —
el intervalo de confianza del estimador cubrirá
un rango amplísimo de valores.
El analista confunde consistencia con eficiencia. Con instrumento débil, IV sigue siendo consistente — en el límite con $n \to \infty$ converge a $\beta$ — pero en muestras finitas tiene dos problemas graves:
Diseña tú el experimento. Queremos saber cómo cambian los
estimadores cuando el efecto causal verdadero es mucho mayor:
beta = 0.40. Mantenemos el resto de parámetros
iguales al caso base.
beta?
¿Y el sesgo relativo (sesgo / valor verdadero)?
Razona usando la fórmula del sesgo.
params_exp3
como nombre del diccionario modificado.
La fórmula del sesgo de OLS es:
El sesgo depende de gamma y delta_U,
no de beta. Por tanto el sesgo absoluto
no cambia al cambiar beta. Sin embargo, el sesgo
relativo (sesgo / $\beta$) sí disminuye:
con $\beta = 0.40$, un sesgo absoluto de ~0.20 representa
el 50% del efecto verdadero, mientras que con $\beta = 0.10$
representaba el 200%.
params_exp3 = params.copy() params_exp3["beta"] = 0.40 df3 = simular_dgp(params_exp3) res3, pi1_3 = estimar_modelos(df3, params_exp3["beta"]) graficar_resultados(df3, res3, params_exp3["beta"])
La barra de OLS seguirá siendo roja:
el sesgo absoluto no cambia con beta,
sigue siendo ~0.20 > 0.08. La barra de Wald
seguirá siendo verde: el estimador
IV sigue siendo consistente independientemente del
valor verdadero de $\beta$.
Incluso con un efecto grande y ambos estimadores apuntando en la misma dirección, la magnitud exacta importa en política pública. Algunos ejemplos: