Sesión 14 (D) — La matriz de regresores en la práctica: verificación con datos simulados

Cerramos la sesión de regresión múltiple volviendo, por primera vez con Gretl, sobre la notación matricial de la lección 10: la matriz de regresores \(\boldsymbol{\mathsf{X}}\), las ecuaciones normales \(\boldsymbol{\mathsf{X}}^\top\boldsymbol{\mathsf{X}}\boldsymbol{\mathop{\widehat{\beta}}}=\boldsymbol{\mathsf{X}}^\top\boldsymbol{y}\), y la condición de ortogonalidad \(\boldsymbol{\mathsf{X}}^\top\boldsymbol{\mathop{\widehat{e}}}=\boldsymbol{0}\). Trabajamos con datos simulados porque ahora no queremos interpretar ningún fenómeno económico; queremos verificar que el álgebra de la lección 10 se cumple exactamente como se demostró: primero en un caso general, y después en el caso especial de regresores mutuamente ortogonales (``familia ortogonal''), donde cada coeficiente se puede calcular de forma aislada.

Objetivo

1. Construir explícitamente, con el comando matrix de Gretl, la matriz de regresores \(\boldsymbol{\mathsf{X}}\) de un modelo con dos regresores no constantes.
2. Calcular \(\boldsymbol{\mathsf{X}}^\top\boldsymbol{\mathsf{X}}\) y \(\boldsymbol{\mathop{\widehat{\beta}}}=(\boldsymbol{\mathsf{X}}^\top\boldsymbol{\mathsf{X}})^{-1}\boldsymbol{\mathsf{X}}^\top\boldsymbol{y}\), comparando el resultado con el que da ols directamente.
3. Verificar \(\boldsymbol{\mathsf{X}}^\top\boldsymbol{\mathop{\widehat{e}}}=\boldsymbol{0}\), en sus dos lecturas (matriz por vector y vector por matriz).
4. Construir un caso de familia ortogonal (regresores mutuamente ortogonales) y comprobar que cada \(\hat\beta_j\) se puede calcular de forma aislada, sin resolver ningún sistema.
5. Verificar, en ese mismo caso, que \(R^2=\rho_{\boldsymbol{z}_2\boldsymbol{y}}^2+\rho_{\boldsymbol{z}_3\boldsymbol{y}}^2\).
6. Ver qué hace Gretl ante la situación contraria, la colinealidad exacta: dos regresores que son la misma variable en dos unidades.

Comandos nuevos de esta práctica. Los datos se generan como en la sesión 7 (nulldata, set seed, uniform, normal, series, scalar, printf) y los modelos se estiman con ols y sus accesores (sesión 11). Lo nuevo es el lenguaje matricial de Gretl:

list Xlist = const x2 x3

una lista de series con nombre propio.

matrix X = { Xlist }

la matriz cuyas columnas son las series de la lista, en ese orden. Aquí es la matriz de regresores \(\boldsymbol{\mathsf X}\) de la lección 10, con \(n\) filas y \(k\) columnas.

X'X

el apóstrofo transpone: X'X es \(\boldsymbol{\mathsf X}^\top\boldsymbol{\mathsf X}\) y X'Y es \(\boldsymbol{\mathsf X}^\top\boldsymbol y\). El asterisco multiplica matrices.

invpd(A)

la inversa de una matriz simétrica definida positiva, como lo es \(\boldsymbol{\mathsf X}^\top\boldsymbol{\mathsf X}\) cuando los regresores son linealmente independientes. det(A) es el determinante.

print M

muestra una matriz con sus filas y columnas.

$coeff sin argumento

el vector columna con todos los coeficientes del último modelo. Con argumento, $coeff(x2), es un solo número (sesión 11).

var(x)

la varianza de una serie con divisor \(n-1\), como sd.

Generamos \(n=40\) observaciones con dos regresores no constantes, \(x_2\) y \(x_3\), deliberadamente correlacionados entre sí (para no obtener, sin advertirlo, el caso especial de la Actividad 4).

Cuarenta observaciones y semilla fija. x2 es uniforme entre \(0\) y \(10\); x3 se construye a partir de x2 más ruido normal, para que ambos estén correlacionados; u es la perturbación. Los tres escalares son los coeficientes verdaderos, y la línea de y aplica el modelo fila a fila. El printf final imprime la correlación entre los dos regresores con la función corr (sesión 7).

en línea de comandos:

nulldata 40 --preserve
set seed 97531

series x2 = uniform(0,10)
series x3 = 0.6*x2 + normal(0,2)
series u  = normal(0,3)

scalar b1 = 5
scalar b2 = 2
scalar b3 = -1
series y = b1 + b2*x2 + b3*x3 + u

printf "Correlacion entre x2 y x3 = %.4f  (no son ortogonales en desviaciones)\n", corr(x2,x3)

Correlacion entre x2 y x3 = 0,5046  (no son ortogonales en desviaciones)

Recuerde la advertencia habitual sobre los vectores de diseño: u es el vector que nosotros usamos para generar los datos; no debe confundirse con el vector de residuos \(\boldsymbol{\mathop{\widehat{e}}}\) que calculará más adelante MCO a partir de los datos simulados.

Construimos explícitamente la matriz de regresores \(\boldsymbol{\mathsf{X}}=[\boldsymbol{1};\;\boldsymbol{x}_2;\;\boldsymbol{x}_3]\) (lección 10) con el comando matrix de Gretl, a partir de una lista de series.

Primero la lista Xlist con los tres regresores, la constante const incluida, y con ella la matriz X; lo mismo para el regresando, que queda en la matriz Y de una sola columna. Después los dos productos de las ecuaciones normales, X'X y X'Y, y print para verlos.

en línea de comandos:

list Xlist = const x2 x3
matrix X = { Xlist }
list ylist = y
matrix Y = { ylist }

matrix XtX = X'X
matrix XtY = X'Y

printf "Matriz X'X (3x3):\n"
print XtX
printf "\nVector X'y (3x1):\n"
print XtY

Matriz X'X (3x3):
XtX (3 x 3)

      40,000       193,05       114,98 
      193,05       1290,6       715,42 
      114,98       715,42       612,46 


Vector X'y (3x1):
XtY (3 x 1)

      472,54 
      2881,9 
      1421,0 

Mire esa matriz: los apuntes de ``Ecuaciones normales'' afirman que cada elemento de \(\boldsymbol{\mathsf{X}}^\top\boldsymbol{\mathsf{X}}\) es un producto escalar entre dos columnas de \(\boldsymbol{\mathsf{X}}\), y eso se puede comprobar aquí mismo, por simple inspección. El elemento \((1,1)\) es \(\langle\boldsymbol1,\boldsymbol1\rangle_e\), es decir, \(n\): debe valer \(40\), y vale \(40\). El elemento \((1,2)\) es \(\langle\boldsymbol1,\boldsymbol{x}_2\rangle_e=\sum_i x_{2i}=n\mu_{\boldsymbol{x}_2}\): divídalo entre \(40\) y obtendrá la media de x2. Lo mismo con el \((1,3)\) y la media de x3. La matriz es la tabla de todos los productos escalares entre regresores, y la primera fila y la primera columna son, sencillamente, sumas.

Ahora resolvemos las ecuaciones normales usando la inversa, tal como aparece en la mayoría de los manuales (lección 10), y comparamos el resultado con lo que da ols directamente.

invpd(XtX)*XtY calcula \((\boldsymbol{\mathsf X}^\top\boldsymbol{\mathsf X})^{-1}\boldsymbol{\mathsf X}^\top\boldsymbol y\). Después estimamos el mismo modelo con ols --quiet y recogemos sus tres coeficientes a la vez en una matriz con $coeff sin argumento. La resta de las dos matrices, componente a componente, debe dar tres ceros.

en línea de comandos:

matrix beta_matricial = invpd(XtX)*XtY
printf "beta calculado como (X'X)^(-1) X'y:\n"
print beta_matricial

ols y const x2 x3 --quiet
matrix beta_ols = $coeff
printf "\nbeta calculado por 'ols' de Gretl:\n"
print beta_ols

matrix diferencia = beta_matricial - beta_ols
printf "\nDiferencia (debe ser practicamente cero):\n"
print diferencia

beta calculado como (X'X)^(-1) X'y:
beta_matricial (3 x 1)

      4,4307 
      2,1139 
    -0,98078 


beta calculado por 'ols' de Gretl:
beta_ols (3 x 1)

      4,4307 
      2,1139 
    -0,98078 


Diferencia (debe ser practicamente cero):
diferencia (3 x 1)

  2,6645e-15 
  4,4409e-16 
 -1,3323e-15 

Recuerde la advertencia de la lección 10: esta fórmula con inversa EXISTE, y es lo que Gretl calcula internamente al ejecutar ols; pero el camino conceptual del curso nunca ha sido, ni será, invertir una matriz: es siempre la condición de ortogonalidad \(\boldsymbol{\mathsf{X}}^\top\boldsymbol{\mathop{\widehat{e}}}=\boldsymbol{0}\), que verificamos en la siguiente actividad.

Recuperamos el vector de residuos de la última regresión y verificamos que es ortogonal a las tres columnas de \(\boldsymbol{\mathsf{X}}\) simultáneamente, en sus dos lecturas equivalentes (lección 10, ``Transposición'').

Los residuos del último modelo, $uhat, pasan a una serie, de la serie a una lista y de la lista a la matriz E de una columna, el mismo camino que seguimos con Y. X'E es el producto de la transpuesta de \(\boldsymbol{\mathsf X}\) por el vector de residuos, tres números en columna; E'X es el mismo producto en el otro orden, una fila.

en línea de comandos:

series ehat = $uhat
list elist = ehat
matrix E = { elist }

matrix XtE = X'E
printf "X'e (columna 3x1, debe ser (aprox.) el vector nulo):\n"
print XtE

matrix eX = E'X
printf "\ne*X (Gretl lo presenta como fila 1x3; mismos tres numeros):\n"
print eX

X'e (columna 3x1, debe ser (aprox.) el vector nulo):
XtE (3 x 1)

 -1,5543e-13 
  5,7732e-13 
 -4,0945e-13 


e*X (Gretl lo presenta como fila 1x3; mismos tres numeros):
eX (1 x 3)

 -1,5543e-13   5,7732e-13  -4,0945e-13 

Observe que XtE y eX contienen los mismos tres números: es la identidad \(\boldsymbol{\mathop{\widehat{e}}}\boldsymbol{\mathsf{X}}=\boldsymbol{\mathsf{X}}^\top\boldsymbol{\mathop{\widehat{e}}}\) de la lección 10, comprobada aquí numéricamente.

Una advertencia sobre lo que ve en pantalla, para que no deshaga lo aprendido en la lección 10. Gretl presenta uno de los dos resultados ``en columna'' y el otro ``en fila'', y el guion ha tenido que escribir E'X, transponiendo lo que el curso llama un vector, porque el lenguaje matricial de Gretl no ofrece otra cosa: para Gretl un vector es una matriz de una sola fila o de una sola columna, y por tanto tiene orientación. Para nosotros no la tiene: un vector es una lista de coordenadas, y solo las matrices se transponen (lección 10, ``Transposición''; y la pregunta 2 del repaso de esa lección, que quizá acaba de responder). De hecho, el sentido de aquella transparencia es precisamente este: \(\boldsymbol{\mathop{\widehat e}}\boldsymbol{\mathsf{X}}\) no exige transponer \(\boldsymbol{\mathop{\widehat e}}\); la transposición recae sobre la matriz, no sobre el vector. Lo que ve en pantalla es una convención de Gretl, no un ingrediente del curso.

Vamos a construir ahora el caso especial de la transparencia ``Familia ortogonal'' de la lección 10. Y lo construiremos cambiando una sola cosa, para que la comparación sea interpretable.

Si cambiáramos varias cosas a la vez y el resultado cambiara, no sabríamos a cuál atribuirlo. Así que dejamos todo exactamente como en la Actividad 1 (las mismas \(n=40\) observaciones, los mismos \(\beta_1=5\), \(\beta_2=2\), \(\beta_3=-1\), y la misma perturbación \(\boldsymbol u\), sin volver a generarla) y cambiamos una sola cosa por decisión nuestra: la correlación entre los dos regresores, que pasará de \(0{,}50\) a cero exacto (con ella cambian, necesariamente, el regresor \(\boldsymbol x_3\) y por tanto el regresando; pero \(n\), los \(\beta\) y la perturbación son los mismos).

Para conseguirlo no hace falta ningún truco nuevo, solo la proyección ortogonal ya conocida: centramos \(\boldsymbol{x}_2\) para que tenga media nula, y a \(\boldsymbol{x}_3\) le quitamos su proyección sobre la constante y sobre \(\boldsymbol{z}_2\), es decir, nos quedamos con sus residuos. Por construcción (ecuaciones normales, lección 10), lo que queda es exactamente ortogonal a \(\boldsymbol1\) y a \(\boldsymbol{z}_2\).

z2 es x2 centrado. Para z3, estimamos con --quiet la regresión de x3 sobre la constante y z2 y nos quedamos con sus residuos, $uhat, que por las ecuaciones normales son ortogonales a los dos regresores. yo se construye con los mismos coeficientes y la misma perturbación u de la Actividad 1. Los cuatro printf comprueban medias, producto escalar y correlación con diez decimales.

en línea de comandos:

series z2 = x2 - mean(x2)        # centrado: media nula
ols x3 const z2 --quiet
series z3 = $uhat                # residuos: ortogonales a 1 y a z2 por construccion

series yo = b1 + b2*z2 + b3*z3 + u   # MISMOS beta, MISMA perturbacion u

printf "media(z2)    = %.10f   (debe ser 0)\n", mean(z2)
printf "media(z3)    = %.10f   (debe ser 0)\n", mean(z3)
printf "suma(z2*z3)  = %.10f   (debe ser 0: ortogonalidad mutua)\n", sum(z2*z3)
printf "corr(z2,z3)  = %.10f   (era 0,5046 entre x2 y x3)\n", corr(z2,z3)

media(z2)    = 0,0000000000   (debe ser 0)
media(z3)    = -0,0000000000   (debe ser 0)
suma(z2*z3)  = 0,0000000000   (debe ser 0: ortogonalidad mutua)
corr(z2,z3)  = 0,0000000000   (era 0,5046 entre x2 y x3)

Observe que la ortogonalidad no es aproximada sino exacta: las cifras son cero hasta el décimo decimal. Es consecuencia de las ecuaciones normales: el residuo de una regresión es siempre, por construcción, ortogonal a todos los regresores de esa regresión. Otra forma habitual de conseguir ortogonalidad exacta es un diseño factorial con valores \(\pm1\) (lo reencontraremos, con otro propósito, al estudiar variables dummy), pero exigiría cambiar los datos, y aquí queremos justamente conservarlos.

Recordando la transparencia ``Familia ortogonal'' de la lección 10: si los regresores son mutuamente ortogonales, cada \(\hat\beta_j\) puede calcularse de forma aislada, sin resolver ningún sistema, con la misma fórmula de Cauchy–Schwarz de la lección 3; y además \(R^2\) se descompone como suma de correlaciones al cuadrado.

Vamos a comprobarlo, pero haciendo los dos cálculos en los dos escenarios: el general de la Actividad 1 y el ortogonal que acabamos de construir. Porque comprobar que la fórmula funciona donde debe funcionar solo es la mitad del experimento; la otra mitad es ver qué ocurre donde no debe.

Dos bloques idénticos, uno por escenario. En cada uno se estima el modelo conjunto con --quiet y se imprimen tres filas con dos columnas: a la izquierda, el coeficiente del modelo conjunto ($coeff(x2)) y el \(R^2\) ($rsq); a la derecha, el cálculo aislado que valdría si los regresores fueran ortogonales: cov(x,y)/var(x), que es la pendiente de la regresión simple (lección 7; las dos funciones dividen por \(n-1\) y el divisor se cancela en el cociente), y la suma de las correlaciones al cuadrado.

en línea de comandos:

printf "================ CASO GENERAL (Actividad 1) ================\n"
printf "correlacion entre los regresores = %.4f\n\n", corr(x2,x3)
ols y const x2 x3 --quiet
printf "beta2:  conjunto = %10.6f     aislado (cov/var) = %10.6f\n", $coeff(x2), cov(x2,y)/var(x2)
printf "beta3:  conjunto = %10.6f     aislado (cov/var) = %10.6f\n", $coeff(x3), cov(x3,y)/var(x3)
printf "R2    : conjunto = %10.6f     suma de rho^2     = %10.6f\n", $rsq, corr(x2,y)^2 + corr(x3,y)^2

printf "\n================ CASO ORTOGONAL (Actividad 4) ==============\n"
printf "correlacion entre los regresores = %.4f\n\n", corr(z2,z3)
ols yo const z2 z3 --quiet
printf "beta2:  conjunto = %10.6f     aislado (cov/var) = %10.6f\n", $coeff(z2), cov(z2,yo)/var(z2)
printf "beta3:  conjunto = %10.6f     aislado (cov/var) = %10.6f\n", $coeff(z3), cov(z3,yo)/var(z3)
printf "R2    : conjunto = %10.6f     suma de rho^2     = %10.6f\n", $rsq, corr(z2,yo)^2 + corr(z3,yo)^2

================ CASO GENERAL (Actividad 1) ================
correlacion entre los regresores = 0,5046

beta2:  conjunto =   2,113860     aislado (cov/var) =   1,675269
beta3:  conjunto =  -0,980779     aislado (cov/var) =   0,222597
R2    : conjunto =   0,771932     suma de rho^2     =   0,651820

================ CASO ORTOGONAL (Actividad 4) ==============
correlacion entre los regresores = 0,0000

beta2:  conjunto =   2,122456     aislado (cov/var) =   2,122456
beta3:  conjunto =  -0,980779     aislado (cov/var) =  -0,980779
R2    : conjunto =   0,835808     suma de rho^2     =   0,835808

La lección 10 puso un ejemplo de regresores linealmente dependientes: la temperatura de un mismo día en grados Celsius y en Fahrenheit, \(F=32+1{,}8\,C\). Cada una varía; juntas no aportan dos direcciones, sino una. Fabriquemos ese caso con los datos de la Actividad 1, tomando x2 como si fuera una temperatura en Celsius.

C es una copia de x2 y F su transformación lineal exacta; corr debe dar \(1\). Con la lista Xcol construimos la matriz Xc y calculamos el determinante de Xc'Xc con det (%g imprime el número en la notación más corta, con exponente si hace falta). El último ols intenta estimar con los dos regresores, esta vez sin --quiet, para leer el mensaje de Gretl.

en línea de comandos:

series C = x2                     # "temperatura" en Celsius
series F = 32 + 1.8*C             # la misma temperatura, en Fahrenheit
printf "corr(C, F) = %.10f\n", corr(C, F)

list Xcol = const C F
matrix Xc = { Xcol }
printf "determinante de X'X con C y F = %g   (una matriz invertible tiene determinante distinto de cero)\n", det(Xc'Xc)

ols y const C F

corr(C, F) = 1,0000000000
determinante de X'X con C y F = 0   (una matriz invertible tiene determinante distinto de cero)

Modelo 1: MCO, usando las observaciones 1-40
Variable dependiente: y
Omitidas debido a colinealidad exacta: F

             coeficiente   Desv. típica   Estadístico t   valor p 
  ----------------------------------------------------------------
  const        3,72826       1,15050          3,241       0,0025   ***
  C            1,67527       0,202541         8,271       5,04e-10 ***

Media de la vble. dep.  11,81341   D.T. de la vble. dep.   6,338798
Suma de cuad. residuos  559,5827   D.T. de la regresión    3,837429
R-cuadrado              0,642903   R-cuadrado corregido    0,633506
F(1, 38)                68,41377   Valor p (de F)          5,04e-10
Log-verosimilitud      -109,5238   Criterio de Akaike      223,0476
Criterio de Schwarz     226,4253   Crit. de Hannan-Quinn   224,2688

El apéndice ``Estructura geométrica de la probabilidad y de la inferencia en el modelo de regresión'' (html) (sección 5) escribe las sumas de cuadrados de la lección 8 como sumas de cuadrados de coordenadas: en una base ortonormal de \(\mathcal L(\boldsymbol 1)^\perp\) cuyos \(k-1\) primeros vectores generen el subespacio de los regresores en desviaciones (el apéndice (html) lo llama \(\mathcal E\)) y cuyos \(n-k\) restantes generen el del residuo (\(\mathcal R\)), \(\mathrm{SEC}\) es la suma de los cuadrados de las \(k-1\) primeras coordenadas del vector de datos en desviaciones, \(\mathrm{SRC}\) la de las otras \(n-k\), y el \(F\) de significación conjunta es el cociente de las dos medias cuadráticas. Vamos a construir esa base con los datos de la Actividad 1 y a comprobarlo con números.

Primero la parte de \(\mathcal E\), ortogonalizando como en la Actividad 4: b_1 es \(\boldsymbol x_2\) en desviaciones dividido por su longitud, y b_2 es la parte de \(\boldsymbol x_3\) en desviaciones ortogonal a b_1, también normalizada (es la parte propia de \(\boldsymbol x_3\), lección 14, dividida por su longitud). Después la parte de \(\mathcal R\), con una regla fija: M es la matriz que proyecta sobre \(\mathcal L(\boldsymbol 1)^\perp\) (la identidad menos la matriz de unos dividida por \(n\)), PR la que proyecta sobre \(\mathcal R\) (se le quitan las dos direcciones de \(\mathcal E\)), y eigensym devuelve sus autovalores en lam y sus autovectores en V; los autovectores de autovalor \(1\) son una base ortonormal de \(\mathcal R\), y selifc se queda con esas columnas. Dos printf comprueban que la base entera es ortonormal y de media cero (las identidades (2) del apéndice (html)). Las coordenadas se obtienen multiplicando la traspuesta de cada bloque por el vector de datos en desviaciones v, y el resto son comparaciones: las dos sumas de cuadrados de coordenadas frente a \(\mathrm{SEC}\) y \(\mathrm{SRC}\) ($ess es la \(\mathrm{SRC}\) de Gretl), su suma frente a \(\mathrm{STC}\), el cociente de medias cuadráticas frente a $Fstat, y la coordenada sobre b_2 dividida por la longitud de la parte propia de \(\boldsymbol x_3\) frente a \(\hat\beta_3\).

en línea de comandos:

ols y const x2 x3 --quiet
scalar n = $nobs
scalar k = 3

series d2  = x2 - mean(x2)
series b_1 = d2 / sqrt(sum(d2^2))
series d3  = x3 - mean(x3)
series e3  = d3 - sum(d3*b_1)*b_1
series b_2 = e3 / sqrt(sum(e3^2))
matrix BE  = {b_1} ~ {b_2}

matrix M   = I(n) - ones(n,n)/n
matrix PR  = M - BE*BE'
matrix V
matrix lam = eigensym(PR, &V)
matrix BR  = selifc(V, (lam .> 0.5)')
printf "vectores de la base: %d en E, %d en R, total %d = n-1\n", cols(BE), cols(BR), cols(BE)+cols(BR)
printf "ortonormalidad: max|B'B - I| = %g ; media cero: max|B'1| = %g\n", max(abs(vec((BE~BR)'(BE~BR) - I(n-1)))), max(abs((BE~BR)'ones(n,1)))

matrix v  = {y} - mean(y)
matrix cE = BE'v
matrix cR = BR'v
printf "\ncoordenadas en E: c1 = %.4f, c2 = %.4f\n", cE[1], cE[2]
printf "primeras coordenadas en R: %.4f, %.4f, %.4f, ...\n", cR[1], cR[2], cR[3]

scalar STC = sum((y - mean(y))^2)
printf "\nsuma de cuadrados de las coordenadas de E = %10.4f   SEC = %10.4f\n", cE'cE, STC - $ess
printf "suma de cuadrados de las coordenadas de R = %10.4f   SRC = %10.4f\n", cR'cR, $ess
printf "suma de todas                             = %10.4f   STC = %10.4f\n", cE'cE + cR'cR, STC
printf "\n(SEC/(k-1)) / (SRC/(n-k)) = %.4f    F de Gretl = %.4f\n", (cE'cE/(k-1)) / (cR'cR/(n-k)), $Fstat
printf "c2 / ||parte propia de x3|| = %.6f    beta3 de Gretl = %.6f\n", cE[2]/sqrt(sum(e3^2)), $coeff(x3)

vectores de la base: 2 en E, 37 en R, total 39 = n-1
ortonormalidad: max|B'B - I| = 1,17961e-15 ; media cero: max|B'1| = 9,55833e-16

coordenadas en E: c1 = 31,7404, c2 = -14,2194
primeras coordenadas en R: -4,8481, 1,3516, -3,9413, ...

suma de cuadrados de las coordenadas de E =  1209,6443   SEC =  1209,6443
suma de cuadrados de las coordenadas de R =   357,3900   SRC =   357,3900
suma de todas                             =  1567,0343   STC =  1567,0343

(SEC/(k-1)) / (SRC/(n-k)) = 62,6162    F de Gretl = 62,6162
c2 / ||parte propia de x3|| = -0,980779    beta3 de Gretl = -0,980779

Hoy hemos verificado con Gretl el álgebra matricial completa de la lección 10: (i) que \(\boldsymbol{\mathop{\widehat{\beta}}}=(\boldsymbol{\mathsf{X}}^\top\boldsymbol{\mathsf{X}})^{-1}\boldsymbol{\mathsf{X}}^\top\boldsymbol{y}\) coincide con lo que calcula ols, y que \(\boldsymbol{\mathsf{X}}^\top\boldsymbol{\mathop{\widehat{e}}}=\boldsymbol{0}\) (que es lo mismo que \(\boldsymbol{\mathop{\widehat{e}}}\boldsymbol{\mathsf{X}}=\boldsymbol{0}\)); (ii) que con regresores mutuamente ortogonales cada coeficiente puede aislarse sin resolver ningún sistema, y \(R^2\) se descompone como suma de correlaciones al cuadrado; (iii) que en una base ortonormal adaptada a los regresores, \(\mathrm{SEC}\) y \(\mathrm{SRC}\) son sumas de cuadrados de coordenadas y el \(F\) es el cociente de sus medias por dimensión.

Y lo hemos verificado de una manera concreta que conviene retener, porque vale para cualquier experimento y no solo para este: comparando dos escenarios que se diferencian en una sola cosa. Si el caso ortogonal hubiera tenido además otro tamaño muestral, otros coeficientes y otro nivel de ruido, habríamos visto cambiar los números sin poder atribuir el cambio a la ortogonalidad. Al dejar fijo todo lo demás, incluida la perturbación, la única explicación posible de la diferencia entre los dos bloques de la Actividad 5 es la que queríamos estudiar.

1. En la Actividad 2, ¿por qué la diferencia entre beta_matricial y beta_ols no es cero, sino un número muy pequeño?
2. ¿Por qué necesitamos que x2 y x3 estén correlacionados en la Actividad 1 para que el ejercicio tenga sentido pedagógico? ¿Qué habría pasado si, por descuido, hubieran resultado casi ortogonales?
3. En la Actividad 4, ¿por qué es importante que \(\boldsymbol{z}_2\) y \(\boldsymbol{z}_3\) tengan media exactamente cero, y no solo estén incorrelados?
4. En la Actividad 5, el cálculo ``aislado'' \(\hat\beta_2\) (cov(x2,y)/var(x2)) sobre los datos de la Actividad 1 (x2, x3, correlacionados) no coincide con el coeficiente de x2 en el modelo conjunto ols y const x2 x3, y el de x3 ni siquiera tiene el mismo signo. ¿Por qué?
5. ¿Diría que el uso de la inversa \((\boldsymbol{\mathsf{X}}^\top\boldsymbol{\mathsf{X}})^{-1}\) en esta práctica contradice la advertencia de la lección 10 de que el curso ``nunca la usa como vía de razonamiento''? ¿Qué diferencia hay entre usarla para verificar una fórmula ya conocida y usarla como método para resolver MCO?
6. En la Actividad 6, ¿por qué el ajuste \(\boldsymbol{\mathop{\widehat y}}\) es único aunque los coeficientes no lo sean? ¿Qué subespacio generan \(\boldsymbol 1\), \(\boldsymbol C\) y \(\boldsymbol F\)?

Para profundizar:

• Véase la lección 10 (lección 10) para la notación matricial completa, el operador selector y la transparencia ``Familia ortogonal''.
• Cottrell, A. y Lucchetti, R. (2023). Gretl User's Guide. Sección sobre construcción y operaciones con matrices a partir de listas de series.

1. ¿Por qué la diferencia no es cero? Por errores de redondeo en aritmética de coma flotante, la misma razón por la que \(\sum\hat e_i\) no era cero en la práctica A de la sesión 11. Invertir una matriz \(3\times3\) y multiplicarla por otro vector implica varias operaciones aritméticas encadenadas, cada una con su pequeñísimo error de precisión; el resultado final coincide con el de ols hasta un margen del orden de \(10^{-10}\) o menor, que es indistinguible de cero a efectos prácticos.
2. ¿Por qué necesitamos correlación entre x2 y x3? Porque el objetivo pedagógico de la Actividad 1 es mostrar el caso general de la regresión múltiple, en el que resolver el sistema completo (o, equivalentemente, usar la matriz inversa) es imprescindible. Si x2 y x3 hubieran resultado casi ortogonales por accidente, el ejercicio se habría parecido demasiado al caso especial de la Actividad 4, y no habría quedado clara la diferencia entre ``hay que resolver el sistema'' (caso general) y ``cada coeficiente se calcula aislado'' (familia ortogonal): precisamente el contraste que se explota, con toda intención, entre las Actividades 1-3 y las Actividades 4-5.
3. ¿Por qué la media exactamente nula, no solo la incorrelación? Porque la fórmula de la transparencia ``Familia ortogonal'' (lección 10) exige que \(\boldsymbol{1},\boldsymbol{z}_2,\boldsymbol{z}_3\) sean mutuamente ortogonales: no basta con que \(\boldsymbol{z}_2\perp\boldsymbol{z}_3\), hace falta además que cada uno de ellos sea ortogonal a la constante \(\boldsymbol{1}\), es decir, que tenga media cero. Sin esa condición, \(\hat\beta_1\) (el coeficiente de la constante) no sería, en general, la media de \(\boldsymbol{y}\), y la familia dejaría de ser ortogonal: la fórmula de la familia ortogonal no se podría aplicar a los tres coeficientes a la vez. Para \(\hat\beta_2\) y \(\hat\beta_3\) bastaría con la incorrelación, porque la constante centra los regresores por su cuenta (lección 10: centrar no cambia los coeficientes de los regresores no constantes); la media exactamente cero es lo que hace que también \(\hat\beta_1\) sea una proyección aislada.
4. ¿Coincide el cálculo aislado con los datos de la Actividad 1? No, y ya no hace falta creerlo: la Actividad 5 lo calcula. La fórmula aislada \(\hat\beta_2=\sigma_{\boldsymbol{x}_2\boldsymbol{y}}/\sigma^2_{\boldsymbol{x}_2}\) solo reproduce el coeficiente del modelo conjunto cuando los regresores son mutuamente ortogonales. Con x2 y x3 correlacionados, el cálculo aislado no tiene en cuenta la información compartida entre ambos y asigna a uno parte de lo que corresponde al otro: \(\hat\beta_2\) pasa de \(2{,}11\) a \(1{,}68\), y \(\hat\beta_3\), y esto es lo llamativo, pasa de \(-0{,}98\) a \(+0{,}22\), cambiando de signo. Es el mismo fenómeno que se observó en la práctica de Ramanathan con sqft y bedrms, aquí con un mundo del que conocemos los coeficientes verdaderos y podemos por tanto ver cuál de los dos cálculos acierta: el modelo conjunto (\(-0{,}98\) frente a \(\beta_3=-1\)).
5. ¿Contradice esta práctica la advertencia de la lección 10 sobre la inversa? No, y la distinción es importante. La advertencia de la lección 10 es que el curso nunca usa la inversa como vía de razonamiento para justificar por qué MCO funciona, ni como herramienta habitual de cálculo: el argumento conceptual sigue siendo, siempre, la condición de ortogonalidad \(\boldsymbol{\mathsf{X}}^\top\boldsymbol{\mathop{\widehat{e}}}=\boldsymbol{0}\). Aquí hemos usado la inversa una única vez, con un propósito muy distinto: verificar numéricamente que una fórmula ya demostrada por otro camino (las ecuaciones normales, obtenidas de la ortogonalidad) da, en efecto, el mismo resultado que ols. Usarla como comprobación puntual de una fórmula ya justificada no es lo mismo que apoyarse en ella para explicar por qué MCO funciona.
6. Ajuste único, coeficientes no. \(\boldsymbol F=32\,\boldsymbol 1+1{,}8\,\boldsymbol C\) es combinación lineal de \(\boldsymbol 1\) y \(\boldsymbol C\), así que \(\mathcal L(\boldsymbol 1,\boldsymbol C,\boldsymbol F)=\mathcal L(\boldsymbol 1,\boldsymbol C)\): un plano, no un subespacio de dimensión tres. La proyección de \(\boldsymbol y\) sobre ese plano es única (es el punto más cercano), pero se puede escribir de infinitas maneras como combinación de tres generadores que solo aportan dos direcciones. Gretl elige una de esas maneras eliminando F; cualquier otra daría el mismo \(\boldsymbol{\mathop{\widehat y}}\), los mismos residuos y el mismo \(R^2\).