Sesión 16 (B) — Condicionar en \(\boldsymbol X\): la misma fórmula, dos experimentos

Qué es condicionar en \(\boldsymbol X\) ya está definido con precisión: la lección 9 presentó \(E[Z\mid\boldsymbol X]\) como la proyección ortogonal de \(Z\) sobre las funciones de \(\boldsymbol X\), y la lección 11 trabajó con esa definición. Esta práctica busca otra cosa: la consecuencia observable de condicionar, es decir, qué distingue en los números la varianza condicional evaluada en la muestra observada de la varianza incondicional.

En la práctica A (actividad 3 (html)) mantuvimos el regresor rooms fijo a lo largo de las 2000 réplicas, y solo dejamos cambiar la perturbación. Lo justificamos así: \(Var[\hat\beta_2\mid\boldsymbol X]\) es una función de \(\boldsymbol X\), y el experimento evalúa esa función en el \(\omega\) observado, donde \(\boldsymbol X\) vale rooms. Condicionar en \(\boldsymbol X\) es otra cosa: proyectar sobre las funciones de \(\boldsymbol X\) (lección 9); de esa proyección sale la función, y el ordenador solo la evalúa. En la imagen de la matriz de copias de la lección 11 (apéndice, sección 3 (html)), fijar el regresor es quedarse con sus \(n\) valores en la columna que ocurrió y recorrer, bloque de filas a bloque de filas, solo copias nuevas de la perturbación.

Esta práctica hace lo contrario. Vamos a repetir el mismo experimento de Montecarlo dos veces:

• Experimento A: el regresor se fija una vez y no cambia; solo cambia la perturbación.
• Experimento B: en cada réplica se sortea también un regresor nuevo.

Si la distinción fuera irrelevante, ambos experimentos darían lo mismo. Veremos que no.

Objetivos

1. Comprobar que la fórmula \(Var[\hat\beta_2\mid\boldsymbol X]=\sigma^2/\sum_i(X_i-\overline{X})^2\) describe el experimento A, no el B.
2. Ver qué consecuencia numérica tiene condicionar en \(\boldsymbol X\), operación que las lecciones 9 y 11 ya definieron con precisión.
3. Ver que la insesgadez, en cambio, se cumple en ambos experimentos, y entender por qué las dos propiedades se comportan de forma distinta.
4. Apreciar por qué el resultado que reportan los programas estadísticos es siempre el condicional.

Advertencia de diseño. Para que la diferencia se vea con claridad, trabajaremos con una muestra deliberadamente pequeña: doce viviendas. Cuanto mayor es la muestra, menos importa la distinción que queremos ilustrar, y con las 506 observaciones de la práctica A costaría apreciarla. En la Actividad 5 repetiremos el experimento con una muestra grande para comprobar eso.

Comandos nuevos de esta práctica. Todo el aparato es el de la práctica A: open, ols --quiet y sus accesores, loop --progressive, store y la apertura del fichero de réplicas con open "@workdir/…". Lo nuevo:

smpl 12 --random

elige al azar doce observaciones del conjunto de datos y restringe la muestra a ellas; a partir de ahí $nobs vale 12 y las funciones sobre series actúan solo sobre esas doce. Con set seed delante, siempre salen las mismas.

matrix pob_rooms = {rooms}

copia la serie en una matriz columna (práctica D de la sesión 14, actividad 2 (html)); la matriz no se ve afectada por smpl, y por eso conserva las 506 viviendas cuando la muestra se restringe a doce.

muniform(r, c)

matriz de \(r\times c\) sorteos uniformes entre 0 y 1.

ceil(x)

redondeo hacia arriba al entero; ceil(N * muniform(12,1)) da doce índices enteros entre 1 y \(N\).

pob_rooms[idx]

extrae de la matriz las filas cuyos índices están en idx; si un índice se repite, la misma vivienda sale dos veces (sorteo con reemplazamiento).

min(x), max(x)

mínimo y máximo de una serie.

Bloques que empiezan con open

cuando la salida de una actividad mezcla escalares teóricos y series de réplicas, los printf van en bloques distintos: el que calcula sobre las réplicas va precedido de open "@workdir/…", que abre su fichero, y desde esa línea los escalares anteriores ya no existen. Por eso los valores teóricos se imprimen siempre antes.

Reutilizamos el mundo artificial de la práctica A (mismos \(\beta_1\), \(\beta_2\) y \(\sigma\)), pero nos quedamos solo con doce viviendas, elegidas al azar. Además, guardamos las 506 viviendas en una matriz antes de restringir la muestra: son la población de la que el experimento B sorteará viviendas nuevas.

Dos bloques. El primero repite el mundo artificial de la práctica A y guarda los 506 valores de rooms en la matriz pob_rooms, la población del regresor, y calcula su varianza var_rooms con el divisor del curso. El segundo imprime la media y la varianza de esa población, restringe la muestra a doce viviendas con smpl 12 --random (después de set seed 160, para que sean siempre las mismas doce) y calcula sobre ellas Sxx_obs, la suma de cuadrados en desviaciones, y Var_b2_en_X_observada, la varianza que predice la lección 11. El «Sxx esperado» que también se imprime es \((n-1)\) veces esa varianza: el valor esperado de Sxx para doce viviendas sorteadas de la población.

en línea de comandos:

open hprice2.gdt --quiet

ols price const rooms --quiet
scalar beta1_verdadero = $coeff(const)
scalar beta2_verdadero = $coeff(rooms)
scalar sigma   = sqrt($ess/$df)

matrix pob_rooms = {rooms}              # las 506 viviendas: la "poblacion" del regresor
scalar N         = $nobs
scalar var_rooms = sum((rooms-mean(rooms))^2)/N

printf "Poblacion del regresor (las %d viviendas de hprice2): media = %.4f   varianza = %.4f\n", N, mean(rooms), var_rooms

set seed 160
smpl 12 --random                        # doce viviendas elegidas AL AZAR
scalar Sxx_obs  = sum((rooms-mean(rooms))^2)
scalar Var_b2_en_X_observada = sigma^2 / Sxx_obs

printf "\nSubmuestra observada (n = %d):\n", $nobs
printf "  media(rooms) = %.4f    varianza(rooms) = %.4f\n", mean(rooms), sum((rooms-mean(rooms))^2)/$nobs
printf "  Sxx de estas %d viviendas (suma de cuadrados en desviaciones) = %.4f\n", $nobs, Sxx_obs
printf "  Sxx esperado para %d viviendas sorteadas de la poblacion     = %.4f\n", $nobs, ($nobs-1)*var_rooms
printf "\nVar[b2 | X] en estas %d viviendas (formula de la leccion 11) = %.1f   (error estandar %.1f)\n", $nobs, Var_b2_en_X_observada, sqrt(Var_b2_en_X_observada)

Poblacion del regresor (las 506 viviendas de hprice2): media = 6,2841   varianza = 0,4927

Submuestra observada (n = 12):
  media(rooms) = 6,3917    varianza(rooms) = 0,2722
  Sxx de estas 12 viviendas (suma de cuadrados en desviaciones) = 3,2670
  Sxx esperado para 12 viviendas sorteadas de la poblacion     = 5,4193

Var[b2 | X] en estas 12 viviendas (formula de la leccion 11) = 13417907,5   (error estandar 3663,0)

Recuerde que Sxx es la abreviatura que introdujimos en la práctica A para \(\sum_i(x_i-\mu_{\boldsymbol x})^2\): el cuadrado de la norma euclídea del vector en desviaciones del regresor, o equivalentemente \(n\,\sigma_{\boldsymbol x}^2\) (evaluada sobre la muestra aleatoria, \(n\,S^2=\sum_i(X_i-\overline X)^2\), lección 11). No es notación del curso, solo una comodidad para los guiones.

Dos detalles del código merecen comentario. El primero es smpl 12 --random: elegimos las doce viviendas al azar, en lugar de tomar las doce primeras. Es importante, porque hprice2 viene ordenado por distritos censales y las primeras viviendas se parecen mucho entre sí; tomarlas directamente nos daría una submuestra atípica por un motivo ajeno a lo que queremos estudiar. El segundo es set seed, que fija qué doce viviendas salen, para que el resultado sea reproducible.

Una cifra de la salida necesita explicación: el «Sxx esperado para 12 viviendas sorteadas de la población», \((n-1)\,\sigma^2_{\boldsymbol x}\), con \(\sigma^2_{\boldsymbol x}\) la varianza de rooms en las 506 viviendas (lección 4). Si cada una de las doce viviendas se sortea de las 506 con reemplazamiento, las doce son copias idénticas e independientes de una variable con esa varianza, y el valor esperado de \(\sum_i(X_i-\overline X)^2\) es \((n-1)\,\sigma^2_{\boldsymbol x}\), la regla de su curso de Estadística para la cuasivarianza. Es el patrón con el que vamos a comparar el Sxx de las doce viviendas que nos han tocado.

Es el experimento de la práctica A sin más cambio que el tamaño muestral: las doce viviendas en lugar de las 506. Cada vuelta sortea la perturbación, construye y, estima con ols --quiet y guarda el coeficiente b_fijo como una fila más de repFijo.gdt. Usamos R = 4000 réplicas, el doble que en la práctica A: la varianza empírica de \(R\) réplicas tiene su propio error de Montecarlo, del orden de \(\sqrt{2/R}\) en términos relativos (un \(3\%\) con 2000, un \(2\%\) con 4000), y queremos que sea pequeño frente a la diferencia entre los dos experimentos, que resultará ser mucho mayor.

en línea de comandos:

scalar R = 4000

set seed 160
loop R --progressive --quiet
    series u = normal(0, sigma)
    series y = beta1_verdadero + beta2_verdadero*rooms + u      # rooms NO cambia
    ols y const rooms --quiet
    scalar b_fijo = $coeff(rooms)
    store "@workdir/repFijo.gdt" b_fijo
endloop

printf "Hecho: %d replicas con el regresor fijo guardadas en repFijo.gdt\n", R

Hecho: 4000 replicas con el regresor fijo guardadas en repFijo.gdt

Todavía no miramos las estimaciones: las compararemos en la Actividad 4 junto a las del experimento B. Lo que tenemos ahora es un fichero con 4000 estimaciones de \(\beta_2\), todas obtenidas con las mismas doce viviendas y distintas perturbaciones.

Ahora decimos ``no solo podría haber sido otra la perturbación; también podrían haber sido otras las doce viviendas''. Para eso hace falta una población de la que sortear viviendas nuevas, y la tenemos: las 506 de hprice2, guardadas en pob_rooms. En cada réplica sorteamos doce de ellas con reemplazamiento: muniform(12,1) da doce números entre 0 y 1, ceil(N * ...) los convierte en doce índices entre 1 y 506, y pob_rooms[idx] extrae esas doce viviendas. Con reemplazamiento, cada sorteo es una copia idéntica e independiente de la distribución de rooms en hprice2: la muestra aleatoria de la lección 11, con esa distribución como población. Una misma vivienda puede salir dos veces, y eso es lo que hace independientes a los doce sorteos.

Antes de incluir el sorteo en el bucle, hacemos uno solo para ver qué cambia. El bloque sortea un xnew, calcula su Sxx y lo imprime junto al de las doce viviendas observadas.

en línea de comandos:

set seed 160
matrix idx  = ceil(N * muniform(12,1))      # doce indices entre 1 y 506
series xnew = pob_rooms[idx]                # doce viviendas de hprice2, con reemplazamiento
printf "Sxx de las doce viviendas observadas = %.4f\n", Sxx_obs
printf "Sxx de doce viviendas sorteadas       = %.4f\n", sum((xnew-mean(xnew))^2)

Sxx de las doce viviendas observadas = 3,2670
Sxx de doce viviendas sorteadas       = 4,5741

El Sxx de las doce viviendas sorteadas no coincide con el de las observadas, y en otro sorteo saldría otro. Eso es lo que el experimento B añade: en cada réplica, además de la perturbación, cambia el regresor y con él su Sxx. En el código del bucle las líneas nuevas son las dos primeras; después y se construye con xnew, la regresión se estima sobre xnew, y store guarda dos escalares por réplica en repVar.gdt: el coeficiente b_var y Sxx_i, la suma de cuadrados en desviaciones del regresor de esa réplica. En el experimento A ese número era una constante, \(3{,}27\); aquí es una variable aleatoria más, y en ella está la diferencia entre los dos experimentos.

en línea de comandos:

set seed 160
loop R --progressive --quiet
    matrix idx  = ceil(N * muniform(12,1))      # doce indices entre 1 y 506
    series xnew = pob_rooms[idx]                # doce viviendas NUEVAS, de hprice2
    series u    = normal(0, sigma)
    series y    = beta1_verdadero + beta2_verdadero*xnew + u
    ols y const xnew --quiet
    scalar b_var = $coeff(xnew)
    scalar Sxx_i = sum((xnew-mean(xnew))^2)     # el Sxx de ESTA replica
    store "@workdir/repVar.gdt" b_var Sxx_i
endloop

printf "Hecho: %d replicas con regresor nuevo guardadas en repVar.gdt\n", R

Hecho: 4000 replicas con regresor nuevo guardadas en repVar.gdt

Ya tenemos los dos ficheros de réplicas. Ahora comparamos tres cosas: lo que predice la lección 11 para nuestras doce viviendas, lo que ha salido con el regresor fijo y lo que ha salido con el regresor sorteado. Tres bloques, cada uno con su salida. El primero imprime la predicción con los escalares que todavía tenemos. El segundo va precedido de open "@workdir/repFijo.gdt", que abre el fichero de réplicas (y borra los escalares, por eso la predicción se imprime antes), y calcula la media y la varianza, con divisor \(n\), de b_fijo. El tercero va precedido de open "@workdir/repVar.gdt", hace lo mismo con b_var y añade, con mean, min y max, cuánto ha variado Sxx_i de una réplica a otra.

en línea de comandos:

printf "=== LO QUE PREDICE LA LECCION 11 (condicional en la X observada) ===\n"
printf "beta_2 verdadero  = %.2f\n", beta2_verdadero
printf "Var[b2 | X]       = %.1f\n\n", Var_b2_en_X_observada

=== LO QUE PREDICE LA LECCION 11 (condicional en la X observada) ===
beta_2 verdadero  = 9119,55
Var[b2 | X]       = 13417907,5

open "@workdir/repFijo.gdt" --quiet

printf "=== EXPERIMENTO A: regresor FIJO ===\n"
printf "media de b2hat    = %.2f\n", mean(b_fijo)
printf "varianza de b2hat = %.1f\n\n", sum((b_fijo-mean(b_fijo))^2)/$nobs

=== EXPERIMENTO A: regresor FIJO ===
media de b2hat    = 9061,70
varianza de b2hat = 13168607,8

open "@workdir/repVar.gdt" --quiet

printf "=== EXPERIMENTO B: regresor REMUESTREADO en cada replica ===\n"
printf "media de b2hat    = %.2f\n", mean(b_var)
printf "varianza de b2hat = %.1f\n", sum((b_var-mean(b_var))^2)/$nobs
printf "\nSxx a lo largo de las replicas del experimento B:\n"
printf "  medio = %.4f    minimo = %.4f    maximo = %.4f\n", mean(Sxx_i), min(Sxx_i), max(Sxx_i)

=== EXPERIMENTO B: regresor REMUESTREADO en cada replica ===
media de b2hat    = 9107,43
varianza de b2hat = 11510016,3

Sxx a lo largo de las replicas del experimento B:
  medio = 5,4394    minimo = 0,4215    maximo = 25,5699

Repetimos el experimento con cuatrocientas viviendas en lugar de doce, dejando todo lo demás igual. La pregunta es si la diferencia entre los dos experimentos sobrevive cuando la muestra es grande, y en qué sentido: la Actividad 4 ha dicho que el signo de la diferencia depende de si nuestra muestra concreta está más o menos dispersa que lo esperado, y con cuatrocientas viviendas esa comparación puede salir al revés.

Dos cambios respecto a las Actividades 2 y 3 merecen explicación. El primero: los dos experimentos van en un solo bucle, que en cada vuelta estima primero con rooms fijo (coeficiente bg_fijo) y después con un regresor nuevo xnew y una perturbación nueva u2 (coeficiente bg_var); así las dos series de réplicas tienen el mismo número de vueltas y la misma semilla, y el cociente entre sus varianzas no depende de eso. El segundo: \(40\,000\) réplicas en lugar de \(4\,000\). Con cuatrocientas viviendas el efecto que buscamos es pequeño, en torno a un \(4{,}7\%\), y con \(4\,000\) réplicas el error de Montecarlo de un cociente de varianzas supera el \(3\%\): repitiendo el experimento con distintas semillas, el cociente oscilaba entre \(0{,}99\) y \(1{,}07\), y con alguna semilla el efecto desaparecía. Con \(40\,000\) se estabiliza por encima de uno. El coste son unos segundos más. Un experimento simulado tiene su propia precisión, y conviene asegurarse de que su error es menor que el efecto que se quiere ver.

Cuatro bloques. El primero reabre hprice2 y recupera los escalares (se perdieron al abrir los ficheros de réplicas), toma cuatrocientas viviendas con smpl 400 --random y ejecuta el bucle; store guarda los dos coeficientes por réplica en repGrande.gdt. El segundo imprime las cifras teóricas, entre ellas el cociente entre el Sxx observado y el esperado, que es lo que cabe esperar del cociente entre las dos varianzas. El tercero, open "@workdir/repGrande.gdt", abre el fichero de réplicas, y el cuarto imprime las dos varianzas empíricas y su cociente.

en línea de comandos:

open hprice2.gdt --quiet
ols price const rooms --quiet
scalar beta1_verdadero = $coeff(const)
scalar beta2_verdadero = $coeff(rooms)
scalar sigma   = sqrt($ess/$df)
matrix pob_rooms = {rooms}
scalar N         = $nobs
scalar var_rooms = sum((rooms-mean(rooms))^2)/N

set seed 160
smpl 400 --random
scalar Sxx_obs  = sum((rooms-mean(rooms))^2)
scalar Var_b2_en_X_observada = sigma^2 / Sxx_obs

set seed 160
loop 40000 --progressive --quiet      # diez veces mas replicas que en la Actividad 2: ver el texto
    series u = normal(0, sigma)
    series y = beta1_verdadero + beta2_verdadero*rooms + u
    ols y const rooms --quiet
    scalar bg_fijo = $coeff(rooms)
    matrix idx  = ceil(N * muniform(400,1))
    series xnew = pob_rooms[idx]
    series u2   = normal(0, sigma)
    series y2   = beta1_verdadero + beta2_verdadero*xnew + u2
    ols y2 const xnew --quiet
    scalar bg_var = $coeff(xnew)
    store "@workdir/repGrande.gdt" bg_fijo bg_var
endloop

printf "Con n = 400:  Sxx observado = %.2f   (esperado: %.2f)\n", Sxx_obs, ($nobs-1)*var_rooms
printf "             Var[b2|X] en las 400 viviendas observadas = %.2f\n", Var_b2_en_X_observada
printf "             cociente Sxx observado/esperado = %.3f  <- lo que cabe esperar de var_B/var_A\n", Sxx_obs/(($nobs-1)*var_rooms)

Con n = 400:  Sxx observado = 205,76   (esperado: 196,57)
             Var[b2|X] en las 400 viviendas observadas = 213046,34
             cociente Sxx observado/esperado = 1,047  <- lo que cabe esperar de var_B/var_A

open "@workdir/repGrande.gdt" --quiet

printf "  regresor FIJO         : var = %.2f\n", sum((bg_fijo-mean(bg_fijo))^2)/$nobs
printf "  regresor REMUESTREADO : var = %.2f\n", sum((bg_var-mean(bg_var))^2)/$nobs
printf "  cociente observado    : %.3f\n", (sum((bg_var-mean(bg_var))^2))/(sum((bg_fijo-mean(bg_fijo))^2))

  regresor FIJO         : var = 211595,99
  regresor REMUESTREADO : var = 224930,68
  cociente observado    : 1,063

1. ¿Qué pregunta responde el experimento A y qué pregunta, distinta, responde el experimento B? Formule ambas con palabras, sin fórmulas.
2. ¿Por qué la insesgadez se cumple en los dos experimentos, mientras que la fórmula de la varianza solo describe el primero? (Pista: ¿en qué paso de la demostración de la lección 11 interviene \(\boldsymbol X\) de forma esencial?)
3. Cuando Gretl le muestra el error estándar de un coeficiente, ¿está respondiendo a la pregunta del experimento A o a la del B? ¿Por qué cree que es así?
4. Imagine que un investigador diseña el experimento y puede elegir los valores del regresor (por ejemplo, las dosis en un ensayo). A la vista de la fórmula de la lección 11, ¿qué le convendría hacer con esas dosis?
5. En la Actividad 5, con \(n=400\), las dos varianzas casi coinciden. ¿Diría entonces que la distinción entre la varianza condicional y la incondicional es irrelevante en la práctica? Argumente a favor y en contra.

Para profundizar:

• Lección 11 (lección 11), sección ``Varianza'', en particular el párrafo final sobre la varianza incondicional y la descomposición de la varianza (apéndice, sección 1.3).
• Lección 9 (lección 9) para la esperanza condicional como proyección.
• Wooldridge, J. M. (2020). Introductory Econometrics, cap. 2, sección 2.5.

1. Las dos preguntas. El experimento A pregunta: ``dadas estas doce viviendas, con estos valores concretos de rooms, ¿cuánto variaría mi estimación si el azar hubiera repartido otras perturbaciones?''. El experimento B pregunta algo más amplio: ``¿cuánto variaría mi estimación entre todos los estudios posibles de doce viviendas, considerando que también podrían haber sido otras las viviendas?''. La primera toma la muestra de regresores como un dato del problema; la segunda la considera parte de lo que podría haber sido distinto.
2. Por qué la insesgadez se mantiene. Porque la demostración de la insesgadez pasa por \(E[\hat\beta_2\mid\boldsymbol X]=\beta_2\,\mathit 1\) y después aplica la ley de las esperanzas iteradas: el resultado condicional no depende de \(\boldsymbol X\), es la misma constante sea cual sea la matriz de regresores, y por tanto su promedio sobre todas las \(\boldsymbol X\) posibles vuelve a ser esa constante. Con la varianza ocurre lo contrario: \(Var[\hat\beta_2\mid\boldsymbol X]=\sigma^2/\mathrm{Sxx}\) sí depende de \(\boldsymbol X\), así que promediar sobre las \(\boldsymbol X\) posibles da un número distinto del que corresponde a una \(\boldsymbol X\) concreta.
3. Lo que reporta Gretl. Responde a la pregunta del experimento A: el error estándar se calcula con la matriz de regresores efectivamente observada, sustituyendo \(\sigma^2\) por su estimación \(\mathfrak s^2\). Tiene sentido que sea así por dos motivos. Uno práctico: no conocemos la distribución de la que salieron los regresores, así que no podríamos promediar sobre ella aunque quisiéramos. Y otro conceptual: la pregunta que le interesa a quien analiza unos datos es cuánta confianza merece su estimación, con su muestra, no cuánta merecería en promedio un estudio hipotético que nunca hizo.
4. Si se pueden elegir los valores del regresor. La fórmula dice que la varianza es \(\sigma^2\) dividido por \(\sum_i(X_i-\overline{X})^2\). Conviene, por tanto, hacer ese denominador lo más grande posible: elegir dosis muy separadas entre sí, en los extremos del rango admisible, en lugar de concentrarlas en torno a un valor central. Es el mismo principio que hemos visto en sentido contrario en esta práctica: nuestras doce viviendas tenían valores de rooms poco dispersos, y eso aumentaba la varianza de la estimación. (Esto supone que el modelo lineal es válido en todo ese rango: separar mucho las dosis mejora la precisión, pero solo si la relación sigue siendo lineal allí.)
5. ¿Es irrelevante con \(n\) grande? A favor de que lo es: numéricamente las dos varianzas se acercan mucho (la diferencia pasa de un \(13\%\) con doce observaciones a un \(6\%\) con cuatrocientas), así que a efectos de calcular un error estándar la distinción cambia poco las cifras. En contra: primero, ``\(n\) grande'' no siempre se tiene, y con muestras pequeñas, que abundan en economía aplicada, la diferencia es la que hemos visto. Segundo, la distinción no es solo numérica, sino conceptual: aclara qué población de muestras estamos imaginando cuando decimos que un estimador ``varía''. Confundir ambas preguntas lleva a interpretar mal qué significa un error estándar.