Sesión 18 (C) — ¿Cubre de verdad el 95 por ciento? La cobertura del intervalo y la frecuencia de detección

La lección 12 distinguió lo que el intervalo de confianza dice y lo que no dice: dice que el suceso ``el intervalo contiene a \(\beta_2\)'' tiene probabilidad \(0{,}95\), con los extremos del intervalo como variables aleatorias; no dice que \(\beta_2\) esté en este intervalo con probabilidad \(0{,}95\). Una probabilidad se manifiesta como frecuencia al repetir, y esta práctica repite.

Reutilizamos el mundo artificial de la práctica B (actividad 2 (html): el panel disperso de la lección 12: doce valores del regresor de \(22{,}5\) a \(77{,}5\), \(\beta_1=10\), \(\beta_2=0{,}5\) y perturbaciones normales con \(\sigma=1{,}1\)) y la mecánica del loop de la práctica A de la sesión 16 (actividad 3 (html)). En cada réplica sorteamos doce perturbaciones nuevas, estimamos, construimos el intervalo y anotamos si contiene al \(0{,}5\) verdadero. Después contamos.

Los datos son simulados porque solo conociendo la verdad se puede saber si el intervalo la contiene.

Objetivos

1. Comprobar que el intervalo al \(95\%\) contiene al verdadero \(\beta_2\) en (aproximadamente) el \(95\%\) de las réplicas, y ver, réplica a réplica, cuáles fallan.
2. Medir el precio de estimar \(\sigma\): qué cobertura se obtiene si se usa \(1{,}96\) en lugar del valor crítico de la \(t_{10}\).
3. Contrastar una hipótesis falsa por poco con \(n=12\) y con \(n=120\): el mismo desvío, medido con un regresor \(\sqrt{10}\) veces más largo.
4. (Apéndice) Ver qué ocurre con la cobertura cuando las perturbaciones no son normales.

Comandos nuevos de esta práctica. La mecánica de loop R --progressive --quiet con store, y la reapertura del fichero de réplicas con open, son las de la práctica A de la sesión 16; critical, $coeff y $stderr, las de las prácticas A y B de hoy; randgen(X, 1), la de la práctica C de la sesión 16. Lo nuevo:

quantile(x, p)

el cuantil de orden p de la serie: el valor que deja por debajo una fracción p de los datos.

%

resto de la división entera: (index - 1) % 12 toma los valores \(0,1,\ldots,11\) y vuelve a empezar.

rech + expresión=

suma la expresión al escalar y guarda el resultado en él; es la forma abreviada de rech = rech + expresión. Como aquí la expresión vale \(1\) o \(0\), el escalar cuenta las réplicas que cumplen la condición.

smpl 1 12

restringe la muestra a las observaciones 1 a 12; smpl full la devuelve entera (sesión 14).

gnuplot ... --time-series --with-lines

dibuja las series frente al número de observación y con líneas en lugar de puntos.

Construimos el mundo, fijamos el valor crítico de la \(t\) con \(n-k=10\) grados de libertad, y repetimos el experimento \(4\,000\) veces guardando, en cada réplica, la estimación, su error estándar, los dos extremos del intervalo, un indicador de si contiene al verdadero \(\beta_2\) y el estadístico \(t\) para la hipótesis cierta \(H_0\!:\beta_2=0{,}5\).

Tres bloques. El primero construye el mundo (doce observaciones, el regresor disperso de la práctica B, los parámetros y el número de réplicas R) y guarda los dos valores críticos: el de la \(t_{10}\) y el \(1{,}96\) de la normal. Sigue el bucle: en cada vuelta sortea doce perturbaciones, forma y, estima con --quiet y guarda en escalares la estimación, su error estándar, los dos extremos del intervalo, el indicador dentro (vale \(1\) si el verdadero \(\beta_2\) está entre los dos extremos; && une las dos desigualdades), el estadístico \(t\) de la hipótesis cierta y el indicador análogo construido con \(1{,}96\); store añade al fichero de réplicas una fila con esos siete valores. El segundo bloque abre ese fichero, que tiene una observación por réplica, y vuelve a definir tc y zc. El tercero imprime el número de réplicas y dos medias: la de dentro, que es la fracción de intervalos que contienen a \(\beta_2\), y la de la condición \(|t|>t_c\), que es la fracción de rechazos.

en línea de comandos:

nulldata 12
setobs 1 1 --cross-section
series x  = 22.5 + 5*(index - 1)              # el regresor disperso de la leccion 12
scalar b1 = 10
scalar b2 = 0.5
scalar sig = 1.1
scalar R  = 4000
scalar tc = critical(t, 10, 0.025)            # valor critico de la t_10 al 5%
scalar zc = 1.96                              # el de la normal, para la Actividad 2

set seed 20261113
loop R --progressive --quiet
    series u  = normal(0, sig)
    series y  = b1 + b2*x + u
    ols y const x --quiet
    scalar b2hat  = $coeff(x)
    scalar ee     = $stderr(x)
    scalar cL     = b2hat - tc*ee                # extremos del IC al 95%
    scalar cR     = b2hat + tc*ee
    scalar dentro = (cL <= b2 && b2 <= cR)       # 1 si el IC contiene al verdadero beta_2
    scalar tstat  = (b2hat - b2)/ee              # t para la H0 CIERTA beta_2 = 0,5
    scalar dentroZ = (b2hat - zc*ee <= b2 && b2 <= b2hat + zc*ee)   # con 1,96 en vez de t_c
    store "@workdir/replicasIC.gdt" b2hat ee cL cR dentro tstat dentroZ
endloop

open "@workdir/replicasIC.gdt" --quiet
scalar tc = critical(t, 10, 0.025)            # al abrir otro fichero Gretl olvida los escalares
scalar zc = 1.96

printf "Replicas                                   : %d\n", $nobs
printf "Intervalos que contienen al verdadero b2   : %.4f   (objetivo: 0,95)\n", mean(dentro)
printf "Rechazos de la H0 cierta al 5%%             : %.4f   (objetivo: 0,05)\n", mean(abs(tstat) > tc)

Replicas                                   : 4000
Intervalos que contienen al verdadero b2   : 0,9520   (objetivo: 0,95)
Rechazos de la H0 cierta al 5%             : 0,0480   (objetivo: 0,05)

Atención a un detalle práctico ya conocido de la práctica A de la sesión 16 (actividad 4 (html)): al abrir el fichero de réplicas, Gretl descarta los escalares que teníamos definidos; por eso tc y zc se vuelven a definir después de abrirlo.

La lección 12 explicó por qué el estadístico sigue una \(t_{n-k}\) y no una normal: en el denominador va \(\mathfrak s\), que fluctúa de muestra en muestra, y eso ensancha las colas (la demostración está en el apéndice, sección 9 (html)). Con \(n-k=10\) la diferencia entre \(t_c=2{,}228\) y \(1{,}96\) es del \(14\%\). Ya hemos guardado, en cada réplica, si el intervalo construido con \(1{,}96\) contenía al verdadero \(\beta_2\); y tenemos los \(4\,000\) estadísticos \(t\) bajo la hipótesis cierta, cuya distribución podemos mirar directamente.

Seis printf. Los dos primeros comparan la cobertura con tc y con \(1{,}96\) (medias de dentro y de dentroZ). El tercero es la cabecera de los cuantiles; el cuarto y el quinto imprimen los cuantiles empíricos del \(2{,}5\%\) y del \(97{,}5\%\) de los \(4\,000\) estadísticos, calculados con quantile, junto a los valores que predicen la \(t_{10}\) y la normal. El sexto calcula la fracción de \(|t|\) mayores que \(1{,}96\). freq tstat dibuja el histograma de los estadísticos.

en línea de comandos:

printf "Cobertura del intervalo con t_c = %.3f : %.4f\n", tc, mean(dentro)
printf "Cobertura del intervalo con 1,96       : %.4f\n", mean(dentroZ)
printf "\nCuantiles empiricos del estadistico t bajo la H0 cierta:\n"
printf "   2,5%%  -> %.3f      (la t_10 predice %.3f; la normal, %.3f)\n", quantile(tstat, 0.025), -tc, -zc
printf "  97,5%%  -> %.3f      (la t_10 predice  %.3f; la normal,  %.3f)\n", quantile(tstat, 0.975), tc, zc
printf "Fraccion de |t| > 1,96                 : %.4f   (con 1,96 se rechazaria la H0 cierta mas del 5%%)\n", mean(abs(tstat) > zc)
freq tstat

Cobertura del intervalo con t_c = 2,228 : 0,9520
Cobertura del intervalo con 1,96       : 0,9188

Cuantiles empiricos del estadistico t bajo la H0 cierta:
   2,5%  -> -2,176      (la t_10 predice -2,228; la normal, -1,960)
  97,5%  -> 2,236      (la t_10 predice  2,228; la normal,  1,960)
Fraccion de |t| > 1,96                 : 0,0812   (con 1,96 se rechazaria la H0 cierta mas del 5%)

Hasta aquí hemos contrastado una hipótesis cierta. Ahora contrastamos una falsa por poco: \(H_0\!:\beta_2=0{,}47\), cuando el verdadero es \(0{,}5\). ¿Con qué frecuencia la detecta el contraste? Lo hacemos dos veces: con las doce observaciones de siempre, y con ciento veinte: el mismo patrón de regresor repetido diez veces, el mismo \(\sigma\), la misma hipótesis. Lo único que cambia es el tamaño de la muestra.

Tres bloques. El primero crea un conjunto de 120 observaciones (--preserve conserva los escalares) con el patrón del regresor repetido diez veces gracias al resto % 12, y vuelve a fijar los parámetros y el valor falso c_falsa. El segundo restringe la muestra a las doce primeras, calcula el valor crítico con los grados de libertad $nobs - 2, fija la semilla, pone a cero dos contadores y repite R veces: sortear, estimar, y sumar a rech un \(1\) si el estadístico de la hipótesis falsa supera el valor crítico y a cov un \(1\) si el intervalo contiene al verdadero \(\beta_2\). Este bucle no lleva --progressive ni store: solo acumula en los dos escalares, y el printf final los divide por R. El tercero repite lo mismo con smpl full, es decir, con las 120 observaciones.

en línea de comandos:

nulldata 120 --preserve
setobs 1 1 --cross-section
series x = 22.5 + 5*((index - 1) % 12)       # el mismo patron, diez veces
scalar b1 = 10                               # los parametros del mundo, de nuevo
scalar b2 = 0.5
scalar sig = 1.1
scalar R = 4000
scalar c_falsa = 0.47

smpl 1 12                                    # una sola vuelta del patron: n = 12
scalar tc = critical(t, $nobs - 2, 0.025)
set seed 20261113
scalar rech = 0
scalar cov = 0
loop R --quiet
    series u = normal(0, sig)
    series y = b1 + b2*x + u
    ols y const x --quiet
    rech += (abs(($coeff(x) - c_falsa)/$stderr(x)) > tc)
    cov  += ($coeff(x) - tc*$stderr(x) <= b2 && b2 <= $coeff(x) + tc*$stderr(x))
endloop
printf "n = %3d:  t_c = %.3f   rechazos de la H0 falsa (0,47) = %.3f   cobertura del IC = %.3f\n", $nobs, tc, rech/R, cov/R

n =  12:  t_c = 2,228   rechazos de la H0 falsa (0,47) = 0,319   cobertura del IC = 0,952

smpl full                                    # las diez vueltas: n = 120
scalar tc = critical(t, $nobs - 2, 0.025)
set seed 20261113
scalar rech = 0
scalar cov = 0
loop R --quiet
    series u = normal(0, sig)
    series y = b1 + b2*x + u
    ols y const x --quiet
    rech += (abs(($coeff(x) - c_falsa)/$stderr(x)) > tc)
    cov  += ($coeff(x) - tc*$stderr(x) <= b2 && b2 <= $coeff(x) + tc*$stderr(x))
endloop
printf "n = %3d:  t_c = %.3f   rechazos de la H0 falsa (0,47) = %.3f   cobertura del IC = %.3f\n", $nobs, tc, rech/R, cov/R

n = 120:  t_c = 1,980   rechazos de la H0 falsa (0,47) = 0,999   cobertura del IC = 0,946

Material complementario, con el mismo estatus que el apéndice matricial (html) de la práctica A de la sesión 16: va más allá de lo que da tiempo en clase y queda aquí para quien quiera completarlo.

La lección 12 fue explícita: la normalidad es el único supuesto nuevo, y su cometido es dar forma a la distribución de \(\hat\beta_2\). La práctica C de la sesión 16 (actividad 1 (html)) mostró que insesgadez y varianza no la necesitan; la pregunta natural es si la cobertura del intervalo, que sí depende de la forma, la necesita. Repetimos la Actividad 1 con perturbaciones \(\chi^2\) de un grado de libertad, centradas y reescaladas a la misma \(\sigma\): una distribución muy asimétrica, la misma con la que la práctica C de la sesión 16 obtuvo (actividad 2 (html)), con \(n=12\), un histograma que el contraste de normalidad rechazaba.

Dos bloques con la misma estructura que los dos últimos bloques de la Actividad 3, cambiando la perturbación: randgen(X, 1) sortea una \(\chi^2\) con un grado de libertad, que tiene esperanza \(1\) y varianza \(2\); restar \(1\) y dividir por \(\sqrt2\) la centra y le da desviación típica \(1\), y multiplicar por sig le da la de nuestro mundo. Los contadores son ahora la cobertura y los rechazos de la hipótesis cierta.

en línea de comandos:

smpl 1 12
scalar tc = critical(t, $nobs - 2, 0.025)
set seed 20261113
scalar cov = 0
scalar rech0 = 0
loop R --quiet
    series u = (randgen(X, 1) - 1)/sqrt(2) * sig     # chi2(1) centrada, con desviacion tipica sig
    series y = b1 + b2*x + u
    ols y const x --quiet
    cov   += ($coeff(x) - tc*$stderr(x) <= b2 && b2 <= $coeff(x) + tc*$stderr(x))
    rech0 += (abs(($coeff(x) - b2)/$stderr(x)) > tc)
endloop
printf "chi2(1) centrada, n = %3d:  cobertura del IC = %.4f   rechazos de la H0 cierta = %.4f\n", $nobs, cov/R, rech0/R

chi2(1) centrada, n =  12:  cobertura del IC = 0,9620   rechazos de la H0 cierta = 0,0380

smpl full
scalar tc = critical(t, $nobs - 2, 0.025)
set seed 20261113
scalar cov = 0
scalar rech0 = 0
loop R --quiet
    series u = (randgen(X, 1) - 1)/sqrt(2) * sig
    series y = b1 + b2*x + u
    ols y const x --quiet
    cov   += ($coeff(x) - tc*$stderr(x) <= b2 && b2 <= $coeff(x) + tc*$stderr(x))
    rech0 += (abs(($coeff(x) - b2)/$stderr(x)) > tc)
endloop
printf "chi2(1) centrada, n = %3d:  cobertura del IC = %.4f   rechazos de la H0 cierta = %.4f\n", $nobs, cov/R, rech0/R

chi2(1) centrada, n = 120:  cobertura del IC = 0,9478   rechazos de la H0 cierta = 0,0523

1. En la Actividad 1, la frecuencia de cobertura no ha salido exactamente \(0{,}95\). ¿Qué esperaría que ocurriese con la diferencia si el experimento se repitiera con \(400\,000\) réplicas en lugar de \(4\,000\)? ¿Y si se repitiera con otras \(4\,000\), cambiando la semilla?
2. Escoja uno de los intervalos de la figura que no contiene a \(0{,}5\). Si ese fuera el único que tuviera, como ocurre siempre con datos reales, ¿qué diría de él? ¿Habría alguna forma de saber que es de los que fallan?
3. En la Actividad 2, ¿por qué usar \(1{,}96\) con \(n-k=10\) hace que el contraste rechace la hipótesis cierta más del \(5\%\) de las veces, y no menos?
4. En la Actividad 3, la cobertura es la misma con doce y con ciento veinte observaciones, pero la detección no. ¿Qué propiedad del intervalo se conserva al cambiar \(n\), y cuál cambia? Relaciónelo con la anchura de los intervalos de la práctica A.
5. Tras el Apéndice: si en un estudio con muestra pequeña las perturbaciones fueran claramente asimétricas, ¿en qué sentido habría que tomar con cautela un valor \(p\) de \(0{,}04\)? ¿Y uno de \(0{,}0001\)?

• Lección 12 (lección 12), secciones ``El intervalo de confianza'' y ``Por qué el valor crítico depende de \(n-k\)''.
• Wooldridge, J. M. (2020). Introductory Econometrics, cap. 4, sección 4.3 (intervalos de confianza) y cap. 5, sección 5.2 (validez aproximada de los contrastes en muestras grandes sin normalidad).
• Manual de Gretl, capítulo sobre loop constructs; funciones quantile y randgen.

1. Más réplicas, otra semilla. La frecuencia de cobertura observada es la media de \(R\) indicadores \(0/1\) independientes, cada uno con probabilidad \(0{,}95\) de valer \(1\); su desviación típica es \(\sqrt{0{,}95\cdot0{,}05/R}\), unos \(0{,}0034\) con \(R=4\,000\) y diez veces menos con \(R=400\,000\). Cabe esperar, por tanto, que la distancia a \(0{,}95\) se reduzca a la décima parte. Con otra semilla y las mismas \(4\,000\) réplicas, la frecuencia sería otra, distinta de la obtenida en el mismo orden de magnitud (unas milésimas): la propia frecuencia empírica tiene su dispersión.
2. Un intervalo que falla, visto en solitario. Nada lo distingue de los que aciertan: tiene el mismo aspecto, la misma anchura aproximada, se ha construido con la misma receta. Con datos reales no hay forma de saber si el intervalo que tenemos delante es uno de los del \(5\%\). Eso es lo que significa que la garantía sea sobre el procedimiento y no sobre el intervalo concreto: sabemos cuántas veces acierta el método, no si ha acertado esta vez.
3. Por qué \(1{,}96\) rechaza de más. Porque \(1{,}96<2{,}228\): la región crítica ``\(\vert t\vert>1{,}96\)'' es más grande que ``\(\vert t\vert>2{,}228\)'', y el estadístico, que sigue una \(t_{10}\), con colas más gruesas que la normal, cae en ella con más frecuencia que el \(5\%\). Usar el valor crítico de la normal es actuar como si \(\sigma\) fuese conocida; como no lo es y se estima, el estadístico se dispersa más de lo que la normal supone, y el contraste, calibrado para una normal, se vuelve demasiado laxo.
4. Lo que se conserva y lo que cambia con \(n\). Se conserva la cobertura: sea cual sea \(n\), el intervalo al \(95\%\) contiene al parámetro el \(95\%\) de las veces, porque su valor crítico se ajusta a los grados de libertad. Cambia la anchura: con más observaciones el error estándar es menor (el regresor en desviaciones es más largo) y el intervalo, más estrecho; por eso un valor falso por poco queda fuera de él con más frecuencia. Es lo que se veía en la práctica A al pasar de doce viviendas a quinientas: la misma garantía del \(95\%\), con un intervalo diez veces más estrecho y, en consecuencia, mucho más informativo.
5. Valores \(p\) con muestra pequeña y perturbaciones asimétricas. El Apéndice muestra que, en ese caso, el nivel real del contraste no coincide con el nominal: aquí el contraste resultó conservador, pero con otras distribuciones podría ser lo contrario, y en cualquier caso el \(0{,}05\) del programa deja de ser una cifra fiable. Un valor \(p\) de \(0{,}04\) está tan cerca del umbral que esa descalibración puede cambiar el veredicto: hay que tomarlo como ``en el límite'', no como ``significativo''. Uno de \(0{,}0001\) está tan lejos del umbral que ninguna descalibración razonable lo convierte en no significativo: la decisión de rechazar es robusta aunque la cifra exacta no lo sea. La regla práctica es que la cautela hace falta cerca del umbral, y que con muestras grandes el problema se desvanece.