Sesión 16 (A) — El experimento de Montecarlo: ¿de verdad está centrado \(\hat\beta_2\)?

La lección 11 demostró dos cosas sobre \(\hat\beta_2\): que es insesgado (\(\mathrm{E}(\hat\beta_2)=\beta_2\)) y que su varianza condicional vale exactamente \[ Var[\hat\beta_2\mid\boldsymbol X]=\frac{\sigma^2}{\sum_i(X_i-\overline{X})^2}. \] Son dos afirmaciones sobre todas las muestras que podríamos haber observado, y por eso resultan difíciles de creer con una sola muestra delante: por definición, de ese conjunto de muestras posibles solo vemos una.

Esta práctica hace trampa —una trampa legítima y muy instructiva—: fabrica ese conjunto de muestras posibles. Inventamos un mundo artificial en el que conocemos \(\beta_1\), \(\beta_2\) y \(\sigma^2\) porque los hemos elegido nosotros; generamos cientos de muestras de ese mundo; estimamos \(\hat\beta_2\) en cada una; y comparamos lo que sale con lo que la lección 11 predice. Es lo que se conoce como un experimento de Montecarlo.

Objetivos

1. Construir un mundo artificial con parámetros conocidos, anclado a un ajuste real sobre hprice2.
2. Generar una muestra simulada a mano y estimarla, para ver con claridad qué es lo que se va a repetir después.
3. Aprender la mecánica completa del comando loop de Gretl: cómo repetir un experimento, cómo guardar resultados de cada repetición y cómo analizarlos después.
4. Comprobar que la media de las estimaciones se acerca a \(\beta_2\) (insesgadez) y que su varianza se acerca a la fórmula de la lección 11.
5. Ver el histograma de \(\hat\beta_2\) y empezar a preguntarse por su forma —una pregunta que esta práctica deliberadamente deja sin responder.

Nota sobre el loop. Ya usamos este comando una vez, en la práctica de verificación geométrica que siguió a la lección 5, pero allí lo presentamos como una ``caja negra'': servía para repetir una comprobación muchas veces, sin explicar su mecánica. Hoy toca abrir la caja.

Para simular necesitamos un mundo del que conozcamos la verdad. Podríamos inventar los números de la nada, pero es más instructivo tomarlos prestados de un ajuste real: así el mundo artificial se parece al que ya hemos estado estudiando.

Partimos de hprice2 y de la regresión price~rooms que ya estimamos en la práctica que siguió a la lección 8. Sus resultados dejan de ser ``estimaciones'' para convertirse, por decreto, en los parámetros verdaderos de nuestro mundo inventado.

o bien teclee en línea de comandos:

open hprice2.gdt --quiet

ols price const rooms --quiet
scalar b1_true = $coeff(const)          # beta_1 "verdadero"
scalar b2_true = $coeff(rooms)          # beta_2 "verdadero"
scalar sigma   = sqrt($ess/$df)         # sigma "verdadera" (raiz de la cuasivarianza)
scalar n       = $nobs
scalar Sxx     = sum((rooms - mean(rooms))^2)
scalar var_teo = sigma^2 / Sxx          # formula de la leccion 11

printf "n            = %d\n", n
printf "beta_1       = %.4f\n", b1_true
printf "beta_2       = %.4f\n", b2_true
printf "sigma        = %.4f\n", sigma
printf "Sxx          = %.4f\n", Sxx
printf "Var[b2 | X]  = sigma^2/Sxx = %.4f\n", var_teo
printf "desv. tipica = %.4f\n", sqrt(var_teo)

n            = 506
beta_1       = -34796,2008
beta_2       = 9119,5484
sigma        = 6620,8652
Sxx          = 249,2872
Var[b2 | X]  = sigma^2/Sxx = 175844,7988
desv. tipica = 419,3385

Antes de repetir nada mil veces, conviene hacer una vez, despacio, lo que después se va a repetir.

Una ``muestra de nuestro mundo'' se construye en dos pasos: primero la parte sistemática \(\beta_1\boldsymbol 1+\beta_2\boldsymbol x\), que es siempre la misma; después una perturbación aleatoria nueva, que es lo único que cambia de una muestra a otra. Sumadas dan un regresando simulado.

o bien teclee en línea de comandos:

set seed 20261030                       # para que el resultado sea reproducible

series u_sim = normal(0, sigma)                    # perturbacion nueva
series y_sim = b1_true + b2_true*rooms + u_sim     # muestra simulada

ols y_sim const rooms
printf "\nbeta_2 verdadero  = %.4f\n", b2_true
printf "beta_2 estimado   = %.4f\n", $coeff(rooms)
printf "diferencia        = %.4f\n", $coeff(rooms) - b2_true

Modelo 1: MCO, usando las observaciones 1-506
Variable dependiente: y_sim

             coeficiente   Desv. típica   Estadístico t   valor p 
  ----------------------------------------------------------------
  const       -34761,1       2711,25         -12,82       9,11e-33 ***
  rooms         9140,64       428,784         21,32       2,34e-72 ***

Media de la vble. dep.  22679,15   D.T. de la vble. dep.   9326,634
Suma de cuad. residuos  2,31e+10   D.T. de la regresión    6769,993
R-cuadrado              0,474145   R-cuadrado corregido    0,473102
F(1, 504)               454,4397   Valor p (de F)          2,34e-72
Log-verosimilitud      -5180,030   Criterio de Akaike      10364,06
Criterio de Schwarz     10372,51   Crit. de Hannan-Quinn   10367,38


beta_2 verdadero  = 9119,5484
beta_2 estimado   = 9140,6352
diferencia        = 21,0868

Aquí está la novedad técnica de la sesión. El comando loop de Gretl repite un bloque de instrucciones un número dado de veces. Su forma general es

loop NUMERO_DE_REPETICIONES --progressive --quiet
    ... instrucciones que se repiten ...
endloop

Dos opciones merecen comentario:

• --quiet evita que Gretl imprima la salida completa de cada una de las repeticiones (imagine 1000 tablas de regresión seguidas).
• --progressive es la que hace el trabajo interesante: le dice a Gretl que va a acumular resultados a lo largo de las repeticiones, en lugar de olvidarlos al terminar cada una. Es lo que permite usar store dentro del bucle.

Y store es la pieza que convierte el experimento en datos analizables:

store "@workdir/fichero.gdt" escalar1 escalar2

Dentro de un bucle --progressive, store no sobrescribe el fichero en cada vuelta: va añadiendo una fila por repetición. Al terminar el bucle tenemos un conjunto de datos nuevo en el que cada observación es una repetición del experimento y cada variable es un resultado que hemos decidido guardar. Después basta con abrirlo, fuera del bucle, y tratarlo como cualquier otro conjunto de datos.

o bien teclee en línea de comandos:

scalar R = 2000                         # numero de replicas

set seed 20261030
loop R --progressive --quiet
    series u    = normal(0, sigma)
    series ysim = b1_true + b2_true*rooms + u
    ols ysim const rooms --quiet
    scalar b1hat = $coeff(const)
    scalar b2hat = $coeff(rooms)
    store "@workdir/replicas.gdt" b1hat b2hat
endloop

printf "Hecho: %d replicas guardadas en replicas.gdt\n", R

Hecho: 2000 replicas guardadas en replicas.gdt

Merece la pena detenerse un segundo en qué cambia y qué no dentro del bucle. Cambia una sola cosa: la perturbación u, que se sortea de nuevo en cada vuelta. No cambian ni rooms, ni b1_true, ni b2_true, ni sigma. Estamos generando, literalmente, las muestras que ``podrían haber salido'' del mismo mundo con el mismo regresor.

Abrimos ahora el fichero de réplicas y lo tratamos como lo que es: un conjunto de datos con \(R\) observaciones, donde cada observación es una estimación que podríamos haber obtenido.

Atención a un detalle práctico: al abrir un conjunto de datos nuevo, Gretl descarta los escalares que teníamos definidos. Por eso imprimimos los valores teóricos antes de abrir el fichero de réplicas, y los empíricos después.

o bien teclee en línea de comandos:

printf "=== LO QUE PREDICE LA LECCION 11 ===\n"
printf "beta_2 verdadero          = %.4f\n", b2_true
printf "Var[b2 | X] = sigma^2/Sxx = %.4f\n", var_teo
printf "desviacion tipica         = %.4f\n\n", sqrt(var_teo)

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

printf "=== LO QUE SALE DE %d REPLICAS SIMULADAS ===\n", $nobs
printf "media de b2hat            = %.4f\n", mean(b2hat)
printf "varianza de b2hat         = %.4f\n", sum((b2hat-mean(b2hat))^2)/$nobs
printf "desviacion tipica         = %.4f\n", sd(b2hat)

=== LO QUE PREDICE LA LECCION 11 ===
beta_2 verdadero          = 9119,5484
Var[b2 | X] = sigma^2/Sxx = 175844,7988
desviacion tipica         = 419,3385

=== LO QUE SALE DE 2000 REPLICAS SIMULADAS ===
media de b2hat            = 9115,1726
varianza de b2hat         = 173769,1343
desviacion tipica         = 416,9605

Tenemos 2000 estimaciones. Veamos cómo se reparten.

o bien teclee en línea de comandos:

freq b2hat --normal --plot=display

La opción --normal superpone al histograma la densidad de una normal con la misma media y la misma varianza que los datos, además de contrastar formalmente si la distribución es normal.

Material complementario. Este apéndice va más allá de lo que da tiempo a hacer en clase; queda aquí para quien quiera completarlo por su cuenta, con el mismo estatus que el apéndice matricial de la práctica sobre datos simulados que siguió a la lección 10.

La última transparencia de la lección 11 enunció, sin demostrarla, la generalización matricial de lo que sí demostramos: \[ Var[\boldsymbol{\mathop{\widehat{\beta}}}\mid\boldsymbol X]=\sigma^2(\boldsymbol X^\top\boldsymbol X)^{-1}. \] Podemos verificarla numéricamente con el mismo experimento, guardando también \(\hat\beta_1\) y comparando la matriz de varianzas y covarianzas empírica de las réplicas con la que predice la fórmula.

o bien teclee en línea de comandos:

open hprice2.gdt --quiet                  # recuperamos los regresores
ols price const rooms --quiet
scalar sigma = sqrt($ess/$df)

matrix X      = {const, rooms}
matrix invXTX = invpd(X'X)
matrix Vteo   = sigma^2 * invXTX
printf "Matriz sigma^2 (X'X)^{-1} predicha por la leccion 11:\n"
print Vteo

matrix B    = {b1hat}~{b2hat}
matrix Vemp = mcov(B)
printf "\nMatriz de varianzas y covarianzas EMPIRICA de las replicas:\n"
print Vemp

Matriz sigma^2 (X'X)^{-1} predicha por la leccion 11:
Vteo (2 x 2)

  7,0306e+06  -1,1050e+06 
 -1,1050e+06   1,7584e+05 


Matriz de varianzas y covarianzas EMPIRICA de las replicas:
Vemp (2 x 2)

  6,9249e+06  -1,0905e+06 
 -1,0905e+06   1,7386e+05 

Compare las dos matrices elemento a elemento. El elemento \((2,2)\) es la varianza de \(\hat\beta_2\) que ya habíamos comprobado en la Actividad 4; el \((1,1)\) es la de \(\hat\beta_1\); y los de fuera de la diagonal son la covarianza entre ambos coeficientes —que, como señalaba la lección 11, es negativa: si la recta se inclina más, tiene que bajar el punto de corte para seguir pasando cerca de la nube de puntos.

1. En la Actividad 2 obtuvimos una estimación distinta del valor verdadero. ¿Significa eso que MCO ha fallado en esa muestra? ¿Qué habría tenido que ocurrir para que coincidieran?
2. ¿Por qué mantenemos rooms fijo dentro del bucle, en lugar de generarlo de nuevo en cada repetición? ¿Con qué expresión de la lección 11 tiene eso que ver?
3. La media de las 2000 estimaciones no coincide con \(\beta_2\). ¿Contradice eso la insesgadez? ¿Qué esperaría que ocurriese con esa diferencia si hiciéramos 200 000 réplicas en lugar de 2000?
4. Suponga que repetimos todo el experimento con un \(\sigma\) diez veces mayor. ¿Qué cambiaría en la media de las estimaciones? ¿Y en su dispersión?
5. Un compañero sugiere que, como el histograma de la Actividad 5 parece una campana, ya podemos calcular la probabilidad de que \(\hat\beta_2\) se aleje más de 500 de su valor verdadero. ¿Qué objeción le haría?

Para profundizar:

• Lección 11 (lección 11): insesgadez y varianza de \(\hat\beta_2\).
• Wooldridge, J. M. (2020). Introductory Econometrics, cap. 2, Teoremas 2.1 y 2.2.
• Manual de Gretl, capítulo sobre loop constructs (en particular, los bucles --progressive).

1. ¿Ha fallado MCO? No. MCO ha hecho exactamente lo que sabe hacer: encontrar la combinación lineal de los regresores más próxima al regresando de esa muestra. Lo que ocurre es que el regresando de esa muestra incluye una perturbación concreta que no conocemos y que MCO no puede separar de la parte sistemática. Según la identidad de la lección 11, \(\hat\beta_2=\beta_2\,\mathit 1+\sum_iW_iU_i\): para que coincidieran tendría que anularse \(\sum_iW_iU_i\), es decir, la perturbación sorteada tendría que resultar ortogonal al regresor centrado. Con una perturbación generada al azar eso no ocurre nunca (salvo casualidad numérica); en la práctica C lo forzaremos a propósito para ver qué pasa.
2. ¿Por qué el regresor se queda fijo? Porque lo que la lección 11 calculó es la varianza condicional \(Var[\hat\beta_2\mid\boldsymbol X]\): la dispersión de las estimaciones para un conjunto dado de valores del regresor. Mantener rooms fijo en todas las réplicas es, literalmente, la versión experimental de condicionar en \(\boldsymbol X\): fijamos la matriz de regresores y dejamos variar solo lo que el modelo considera aleatorio, la perturbación. La práctica B explora qué ocurre si dejamos de condicionar.
3. La media no coincide. No lo contradice. La insesgadez es una afirmación sobre la esperanza, es decir, sobre el promedio en el límite de infinitas muestras posibles. Con \(R\) réplicas obtenemos un promedio muestral de estimaciones, que es él mismo una variable aleatoria con varianza aproximadamente \(Var[\hat\beta_2\mid\boldsymbol X]/R\). Al multiplicar \(R\) por cien, esa varianza se divide por cien y la desviación típica por diez: cabe esperar que la diferencia con \(\beta_2\) se reduzca aproximadamente a la décima parte.
4. Un \(\sigma\) diez veces mayor. La media de las estimaciones seguiría centrada en \(\beta_2\): la insesgadez no depende de \(\sigma\) en absoluto (la demostración de la lección 11 solo usa \(E[U_i\mid\boldsymbol X]=\mathit 0\)). La dispersión, en cambio, se multiplicaría por diez la desviación típica —y por cien la varianza—, porque \(\sigma^2\) aparece como factor en el numerador de la fórmula. Dicho de otro modo: más ruido no descentra el estimador, solo lo vuelve menos preciso.
5. La objeción al compañero. Que la forma de campana del histograma no la hemos descubierto, sino impuesto: hemos generado las perturbaciones con normal(0, sigma) porque así lo decidimos al construir el mundo artificial. Nada de lo demostrado en la lección 11 garantiza esa forma; allí se demostraron la media y la varianza, y ambas se obtuvieron sin suponer nada sobre la distribución de \(U\). Para calcular probabilidades hace falta un supuesto adicional sobre esa distribución, y ese supuesto —con sus consecuencias— es justamente el asunto de la lección 12. La práctica C de hoy muestra, además, que con perturbaciones muy poco normales la media y la varianza siguen cumpliéndose igual de bien.